CN103942381A - 用于飞机铝合金结构性能预测的状态近场动力学方法 - Google Patents

Abstract

一种力学数值模拟技术领域的用于预测飞机铝合金结构性能的状态近场动力学方法，首先将铝合金结构离散成一系列空间域内的物质点；其次根据近场动力学状态线性化理论推导出节点刚度矩阵，并在此基础上改进Kilic针对近场动力学键理论提出动力松弛法，求解计算所需参量；然后将铝合金结构总的位移载荷分为若干增量步，对每个增量步采用改进的动力松弛法迭代计算；结合相对和绝对准则改进不平衡力准则，选用该准则判定系统的平衡与否，若平衡则进入下一增量步，直至加载加速，否则进入下一轮迭代，最终得到铝合金结构的性能数据。通过本发明能够实现运用近场动力学状态理论的常规材料模型求解复杂平衡路径准静态问题。

Description

用于飞机铝合金结构性能预测的状态近场动力学方法

技术领域

本发明涉及的是一种力学数值模拟技术领域的方法，具体是一种基于近场动力学状态理论常规材料模型的准静态数值模拟方法。

背景技术

近场动力学(Peridynamics，简称PD)是近期发展起来的一种多尺度力学方法，能有效解决宏、微观力学领域中的不连续问题。它的优势源于理论本身的特点：它将物体离散成一系列空间域内的物质点，基于非局部作用思想建立物质点间的作用关系，通过空间积分方程描述物质点的运动，可采用统一的模型和求解体系描述原子尺度到宏观尺度的力学行为，不连续现象自然产生，突破了连续性假设和空间微分方程在不连续问题上出现的求解瓶颈。PD方法兴起于2000年，目前存在两个理论分支——近场动力学键理论和近场动力学状态理论，其中：近场动力学状态理论突破了近场动力学键理论对传统材料建模的局限，更具先进性。按物质点间作用力方向的不同，可将近场动力学状态理论中的材料模型分为两类：常规材料和非常规材料，根据Seleson P,Parks M L.“On the role of the influence function in theperidynamic theory”([J].Journal for Multiscale Computational Engineering,2011,9(6):689–706.)，其中：常规材料模型与近场动力学键理论中的材料模型在一定条件下具有一致性。近场动力学状态理论提出于2007年，目前关于常规材料模型准静态问题模拟的数值方法研究很少，Mitchell J A.在“A nonlocal,ordinary,state‐based plasticity model for peridynamics”([R].Sandia National Laboratory Report,SAND2011‐3166,Albuquerque,New Mexico,2011.)中只提到了一种类似于混合法的数值求解思想，但没有给出具体实施方法及应用实例，且该方法对一些特殊问题的求解会失效，如极值型失稳问题。Kilic B,Madenci E.在“An adaptive dynamicrelaxtion method for quasi‐static simulation using the peridynamic theory”([J].Theoretical andApplied Fracture Mechanics.2010,53:194‐204.)中提出的动力松弛法能解决此类问题，但该方法仅限于近场动力学键理论。

发明内容

本发明针对现有技术存在的上述不足，提出一种用于飞机铝合金结构性能预测的状态近场动力学方法，通过本发明能够实现运用近场动力学状态理论的常规材料模型求解复杂平衡路径准静态问题。

本发明是通过以下技术方案实现的，本发明首先将结构离散成一系列空间域内的物质点；其次根据近场动力学状态线性化理论推导出节点刚度矩阵，并在此基础上改进Kilic针对近场动力学键理论提出动力松弛法，求解计算所需参量；然后将结构总的位移载荷分为若干增量步，对每个增量步采用改进的动力松弛法迭代计算，增量步开始时根据节点的空间位置选取初始位移场，使结构的空间构型接近新的平衡状态以提高算法的收敛速度；结合相对和绝对准则改进不平衡力准则，选用该准则判定系统的平衡与否，若平衡则进入下一增量步，直至加载结束，否则进入下一轮迭代；最终得到铝合金结构的性能数据。

本发明具体包括以下步骤：

步骤一、将结构离散成一系列空间域内的物质点：将物体离散成均匀的正方体晶格，节点(积分计算点)置于正方体晶格的中心。

步骤一中，当采用薄板铝合金结构时，将其离散成板状晶格以减少节点数量提高计算效率。

所述的薄板铝合金结构中的板状晶格在其他方向的尺寸一致，以使节点键族分布均匀。

所述的正方体晶格在各方向的晶格数量不少于2，以避免节点刚度矩阵和某方向相关的元素都为0导致的计算失败。

步骤二、根据近场动力学状态线性化理论推导出节点刚度矩阵，具体步骤包括：

设结构离散后，物质点x附近有N_x个节点q_i(1≤i≤N_x)，满足|q_i-x|<2δ，根据PD状态线性化理论，导出节点x的增量平衡方程为：

$\left[\begin{array}{cccc}C(x,x)& C(x,q_1)& ...& C(x,q_{N_x})\end{array}\right]\left[\begin{array}{c}u\left(x\right)\\ u\left(q_1\right)\\ .\\ .\\ .\\ u\left(q_{N_x}\right)\end{array}\right]+b\left(x\right)=0,$ 其中：物质点x的节点刚度矩阵为记n为任意单位矢量，C(x,q_i)n表示保持除q_i外的其他物质点不动，将物质点q_i移动位移矢量n对物质点x产生的内力密度，C(x,q_i)的计算方式如下：

其中p点满足

{p||p-x|<δ且|p-q|<δ}，[]<>为模量双状态，其定义可见文献Silling S A.“Linearizedtheory of peridynamic states”([J].Journal of Elasticity,2010,99:85‐111.)中的公式(32)及文献Mitchell J A.“A nonlocal,ordinary,state‐based plasticity model for peridynamics”([R].SandiaNational Laboratory Report,SAND2011‐3166,Albuquerque,New Mexico,2011.)的公式(26)。

步骤三、根据步骤二中得到的C(x,x)计算虚拟密度矩阵虚拟密度矩阵为对角阵，其对角元素计算如下：

其中：上标0表示在结构的初始构型中求解该参量。

步骤四、将结构总的位移载荷分为若干增量步，对每个增量步采用改进的动力松弛法迭代计算。

步骤五、增量步开始时，根据节点的空间位置，插值得到初始位移场。

在第n个增量步初，在结构位移约束处，取j方向位移固定的约束节点x_b1和j方向位移已知的加载节点x_b2，根据节点的空间位置进行插值，得到j方向位移自由的节点x在此增量步初该方向的初始位移增量： ${\mathrm{\&Delta;u}}_j^{n,0}\left(x\right)=\frac{{\mathrm{\&Delta;u}}^n{\left(x_{b2}\right)}_j}{{(x_{b2}-x_{b1})}_j}{(x-x_{b1})}_j$

步骤六、结合相对和绝对准则改进不平衡力准则，选用该准则判定系统的平衡与否，若平衡则进入下一增量步，直至加载结束，否则进入下一轮迭代。

改进后的不平衡力准则为： $\left\{\begin{array}{cccc}\underset{j}{\max }\left|\frac{L_{i,j}^{n,k}}{L_{i,j}^{n,0}}\right|\&le;{\mathrm{tol}}_{r,j}& & & k<n_b\\ \underset{j}{\max }\left|\frac{L_{i,j}^{n,k}}{L_{i,j}^{n,0}}\right|\&le;{\mathrm{tol}}_{r,j}& \mathrm{or}& \underset{j}{\max }\left|L_{i,j}^{n,k}\right|\&le;{\mathrm{tol}}_{a,j}& k\&GreaterEqual;n_b\end{array}\right.,$ 其中：表示在第n个增量步初节点x_i的节点在第j个自由度上的不平衡内力，表示第n个增量步在第k次迭代后编号为i的节点在第j个自由度上的不平衡内力。n_b表示平衡准则转换的临界次数，tol_r,j表示节点在第j个自由度上的相对平衡容差，tol_a,j表示节点在第j个自由度上的绝对平衡容差。

考虑到权函数对局部作用的影响，所述的节点在不同自由度上可采用不同的相对或绝对平衡容差。

步骤七、根据步骤六中迭代平衡后的系统状态，可输出节点变量在某增量步的场量云图及节点变量随加载变化的历史曲线。优选地，当在结构稳定性分析中，则输出铝合金结构的位移场云图及位移载荷曲线，得到铝合金结构的失稳模态及失稳载荷性能数据。

技术效果

与现有技术相比较，本发明技术效果包括：

1.本发明能求解复杂平衡路径下基于上述理论的准静态模拟问题，如极值型失稳问题。

2.本发明基于上述理论提出的薄板结构离散方法，能合理分布节点，有效提高计算效率。

3.本发明通过推导节点刚度矩阵改进的Kilic基于近场动力学键理论的动力松弛法，能有效适用于近场动力学状态理论的求解。

4.本发明提出的初始位移选择办法能使结构的空间构型接近新的平衡状态，提高算法的收敛速度。

5.本发明改进的迭代收敛准则，能改善绝对准则对求解精度的影响及相对准则对求解速度的影响，有效提高求解效率。

附图说明

图1为本发明算法流程示意图。

图2为模型初始几何缺陷云图。

图3为载荷随轴向位移变化的曲线(平板)。

图4为中心点挠度随载荷变化的曲线(平板)。

图5为平板初次屈曲模态。

图6为平板二次屈曲模态图。

图7为平板三次屈曲前模态图。

图8为平板三次屈曲后模态图。

图9为面内刚度随轴向位移变化的曲线(平板)。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1

如图1所示，本实施例包括以下步骤：

a)输入模型参数

本实施例为平板的轴压失稳模拟，平板长300mm，宽200mm，厚2mm，杨氏模量为69GPa，泊松比为0.33，载荷形式为沿长度方向的压缩位移1.2mm，边界条件为四边简支。权函数形式为键长的倒数。

b)模型离散化

本实施例沿长度方向划分了30个晶格，沿宽度方向划分20个晶格，沿厚度方向划分2个晶格，方板晶格尺寸为10mm×10mm×1mm，节点置于方板晶格中心，节点的近场范围为14.2mm(约为方板晶格的对角线长度)。

c)初始化节点信息

本实施例为稳定性问题，考虑到初始几何缺陷在生产加工过程中不可避免，且对结构稳定性有一定影响，因此本实施例引入了初始几何缺陷，初始几何缺陷通过有限元软件ABAQUS的BUCKLE分析模块生成，有限元模型的节点与本实施例的节点一一对应，材料属性和边界设置同本实施例一致。引入本实施例的初始几何缺陷为有限元模型的一阶屈曲模态，如图2，且按厚度的十分之一进行缩放，即几何缺陷缩小为图2所示缺陷的0.2倍。

引入几何缺陷后，重新计算节点的初始位置，并根据近场范围，计算各节点的键族信息。

d)计算节点刚度

e)计算节点虚拟密度矩阵

f)将载荷分100个增量步进行加载，平衡标志设置为否，增量步计算器设置为1，平衡迭代计算器设置为0。

g)判断平衡迭代计数器是否为0，若为0则计算该增量步的初始位移场增量，更新节点位移初始位移场增量，否则根据上一轮迭代计算得到的位移场增量，更新节点位移。

h)由新的节点位移，根据文献Kilic B,Madenci E.在“An adaptive dynamic relaxtionmethod for quasi‐static simulation using the peridynamic theory”([J].Theoretical and AppliedFracture Mechanics.2010,53:194‐204.)中的式(81～83)、(85～87)、(112～113)更新节点的键族信息。

i)计算节点的不平衡力。

j)根据不平衡力准则更新平衡标志，其中：n_b设为1000，各方向的绝对容差因子设为0.001，相对容差因子设为0.01。

k)根据平衡标志判断，若为“是”，则执行步骤m，若为“否”，则执行步骤l。

l)按动力松弛法计算节点本轮迭代的位移增量，平衡迭代计算器增加1，随后执行步骤g，进入下一轮迭代。

m)平衡迭代次数设置为0，增量步计数器增加1。

n)判断增量步计算器是否指定步数，本实施例中，若增量步计数器为101，则计算结束，否则执行步骤g，进入下一增量步。

本实施例模拟得到平板简支加载边沿长度方向的外力(即轴向外力)随轴向位移的变化如图3，中心点挠度随轴向外力的变化如图4中菱形标记。

从图中可以看出：a点为平板的初次屈曲点，中心点挠度变化率在初次屈曲前后有较大变化，初次屈曲载荷为9.608kN，初次屈曲模态如图5所示；b点为平板的二次分支屈曲点，中心点挠度达到局部最大值，二次屈曲载荷为35.2kN，二次屈曲模态如图6所示；从c点到d点，平板发生了三次跳跃屈曲，三次屈曲载荷为39.6kN，c的模态如图7所示，d点的模态如图8所示。面内压缩刚度随轴向位移的变化如图9所示，可以看出初次屈曲和三次屈曲时，面内刚度都发生了突变。

技术效果评价

根据飞机设计手册总编委会《飞机设计手册第9册》([M].航空工业出版社,2001:391‐393.)中给出了轴压载荷作用下矩形平板弹性稳定性临界应力σ_cr的计算公式： 其中：E和μ_e为材料的弹性模量和泊松比，δ和b为平板的厚度和加载边宽度。K_c为压缩临界应力系数，与平板的边界条件及长宽比有关。

四边简支条件下，当平板的长宽比大于1时，K_c取4，将E=69×10³MPa，μ_e＝0.33，δ＝2mm，b=200mm代入σ_cr的计算公式中，得到临界应力σ_cr=25.47MPa，临界应力与加载边截面面积相乘后得到临界压力，临界压力为10.190kN，见图4的三角形标记，数值模拟得到的初次屈曲载荷与该值相近，相差约5.7%，与目前薄板屈曲的研究成果一致。且初次屈曲载荷的三倍以内模态形式不变，这与沈惠申.“板壳后屈曲行为”([M].上海科学技术出版社,2002:228‐295.)以及周祥玉,范祖尧.“受压薄板后屈曲性态的试验研究”([J].上海交通大学学报,1992,26(4):14‐18.)中记载的大多薄板轴压屈曲试验的研究成果一致。

Claims (6)

1.一种用于预测飞机铝合金结构性能的状态近场动力学方法，其特征在于，首先将铝合金结构离散成一系列空间域内的物质点；其次根据近场动力学状态线性化理论推导出节点刚度矩阵，并在此基础上改进Kilic针对近场动力学键理论提出动力松弛法，求解计算所需参量；然后将铝合金结构总的位移载荷分为若干增量步，对每个增量步采用改进的动力松弛法迭代计算；结合相对和绝对准则改进不平衡力准则，选用该准则判定系统的平衡与否，若平衡则进入下一增量步，直至加载加速，否则进入下一轮迭代，最终得到铝合金结构的性能数据。

2.根据权利要求1所述的方法，其特征是，所述方法具体包括以下步骤：

步骤一、将物体离散成均匀的正方体晶格，节点即积分计算点置于正方体晶格的中心；

步骤二、物质点x附近有N_x个节点q_i(1≤i≤N_x)，满足|q_i-x|<2δ，根据PD状态线性化理论，导出节点x的增量平衡方程为：

$\left[\begin{array}{cccc}C(x,x)& C(x,q_1)& ...& C(x,q_{N_x})\end{array}\right]\left[\begin{array}{c}u\left(x\right)\\ u\left(q_1\right)\\ .\\ .\\ .\\ u\left(q_{N_x}\right)\end{array}\right]+b\left(x\right)=0,$ 其中：物质点x的节点刚度矩阵为记n为任意单位矢量，C(x,q_i)n表示保持除q_i外的其他物质点不动，将物质点q_i移动位移矢量n对物质点x产生的内力密度；

步骤三、根据步骤二中得到的C(x,x)计算虚拟密度矩阵虚拟密度矩阵为对角阵，其对角元素为：其中：上标0表示在结构的初始构型中求解该参量；

步骤四、将结构总的位移载荷分为若干增量步，对每个增量步采用改进的动力松弛法迭代计算；

步骤五、增量步开始时，根据节点的空间位置，插值得到初始位移场，具体为：在第n个增量步初，在结构位移约束处，取j方向位移固定的约束节点x_b1和j方向位移已知的加载节点x_b2，根据节点的空间位置进行插值，得到j方向位移自由的节点x在此增量步初该方向的初始位移增量： ${\mathrm{\&Delta;u}}_j^{n,0}\left(x\right)=\frac{{\mathrm{\&Delta;u}}^n{\left(x_{b2}\right)}_j}{{(x_{b2}-x_{b1})}_j}{(x-x_{b1})}_j;$

步骤六、结合相对和绝对准则改进不平衡力准则，选用该准则判定系统的平衡与否，若平衡则进入下一增量步，直至加载结束，否则进入下一轮迭代，具体为：改进后的不平衡力准则为： $\left\{\begin{array}{cccc}\underset{j}{\max }\left|\frac{L_{i,j}^{n,k}}{L_{i,j}^{n,0}}\right|\&le;{\mathrm{tol}}_{r,j}& & & k<n_b\\ \underset{j}{\max }\left|\frac{L_{i,j}^{n,k}}{L_{i,j}^{n,0}}\right|\&le;{\mathrm{tol}}_{r,j}& \mathrm{or}& \underset{j}{\max }\left|L_{i,j}^{n,k}\right|\&le;{\mathrm{tol}}_{a,j}& k\&GreaterEqual;n_b\end{array}\right.,$ 其中：表示在第n个增量步初节点x_i的节点在第j个自由度上的不平衡内力，表示第n个增量步在第k次迭代后编号为i的节点在第j个自由度上的不平衡内力，n_b表示平衡准则转换的临界次数，tol_r,j表示节点在第j个自由度上的相对平衡容差，tol_a,j表示节点在第j个自由度上的绝对平衡容差；

步骤七、根据步骤六中迭代平衡后的系统状态，输出节点变量在某增量步的场量云图及节点变量随加载变化的历史曲线。

3.根据权利要求1或2所述的方法，其特征是，当采用薄板铝合金结构时，将其离散成板状晶格以减少节点数量提高计算效率。

4.根据权利要求3所述的方法，其特征是，所述的薄板铝合金结构中的正方体晶格在其他方向的尺寸一致，以使节点键族分布均匀。

5.根据权利要求2所述的方法，其特征是，所述的正方体晶格在各方向的晶格数量不少于2，以避免节点刚度矩阵和某方向相关的元素都为0导致的计算失败。

6.根据权利要求1或2所述的方法，其特征是，所述的节点在不同自由度上可采用不同的相对或绝对平衡容差。