基于有限元法的轴对称结构极限下限分析

王 崴，王颜辉，孙 万，李睿智，巩少鹏

（山西能源学院，山西 晋中 030600）

0 引言

轴对称结构广泛应用于航空航天、化工、机械等工程领域中。通过确定各类轴对称结构的极限载荷，可以为工程设计和安全评定提供准确可靠的理论依据和参考[1-3]。传统的“许用应力法”对结构极限承载能力的评估往往趋于保守，而弹塑性增量法由于计算难度大、计算成本高等缺点使得其很难应用于实际工程中。相较于前两种方法，极限分析能够在考虑材料塑性性质的基础上得到真实反映材料安全裕度的重要参数，而且计算过程不涉及结构的加载历史，执行起来简单易行，是研究轴对称结构承载能力的一种直接有效的方式。

目前，将数值方法与数学规划理论相结合已成为分析极限分析问题的主流方法。这种方法将极限分析归纳为一个非线性数学规划问题，但维数障碍使得该问题的求解异常困难。众多学者为克服维数障碍问题进行了研究。李春光等[4]对四边形网格的平衡方程进行积分，减少了自变量和约束方程的数量；钱向东等[5]对路径跟踪内点法进行改进，在算法中引入子迭代程序，使大量计算可以在单元一级完成；秦方等[6]将下限分析列式转化为圆锥二次优化问题进行数值求解。

本文基于塑性极限下限分析理论，结合有限元离散技术，构造了求解轴对称极限下限问题的相应算法。通过减缩基技术将下限分析列式归结为一系列数学规划子问题，并在每个子问题中通过弹塑性增量分析得到的结果形成自平衡应力场搜索子空间。根据数学规划问题中变量少、约束条件多的特点，采用复合形法进行求解，并编制了相应的计算程序。最后通过对两个轴对称结构极限下限载荷系数的计算和对比，验证了本文方法的可行性和有效性。

1 数学规划格式

现考虑一个理想弹塑性轴对称结构V，为简便计算，取轴对称结构中的任一对称面S进行分析。极限下限定理可以表述为[7-9]：结构未达到极限状态的充要条件是，存在一个自平衡应力场ρij（x），与虚拟弹性应力场（x）相叠加后产生的总应力场处处不违反屈服条件f[·]。而极限下限载荷就是结构在未达到极限状态时所能承受的最大载荷，由此得到其数学规划格式为：

其中χ为待求的极限载荷系数，nj表示边界单位外法线向量，σs是Mises 屈服应力。表达式（2）为结构应满足的屈服条件，（3）和（4）分别表示自平衡应力场需要满足的应力边界条件和平衡方程。

对于一个连续体结构，为减少数学规划式（1）～（4）的未知数和约束条件的数量，必须对结构进行离散化处理。根据虚位移原理[10]，应力边界条件（3）和平衡方程（4）的等效积分弱形式表示为：

采用有限元法进行离散，假定将对称面S划分为N个单元的的组合，则每个单元中的应变列向量εe的计算a表达式为：

式中，Be为单元应变位移关系矩阵，u表示单元节点位移列向量，根据式（6）可以得到：

将单元节点应变向量εe组合成为总体节点应变向量ε，并代入式（5），由Gauss 积分法得到：

式中，U表示结构的总体节点位移向量，B为扩展之后与U同阶的应变位移关系矩阵，N是离散之后结构中的单元总数，G为每个单元内部Gauss 积分点的数目，I表示全部Gauss 积分点的集合，（ρe）i为单元e中第i个Gauss 积分点的自平衡应力，ri和分别为单元e中第i个Gauss 积分点的径向半径和Jacobi 矩阵的值。考虑到δU的任意性，式（8）进一步写为：

令xi表示离散之后第e个单元中第i个Gauss积分点的坐标值，则虚拟弹性应力场。综合上述分析过程，离散之后的轴对称极限下限分析数学规划格式为：

2 减缩基技术

求解非线性规划问题（12）的最大困难在于未知数和约束条件的数量较为庞大，因此本文采用减缩基技术[11]将自平衡应力场表示为一组自平衡应力基矢量的线性组合：

其中M1，M2，…，MK表示基矢量对应的待定系数。通过减缩基技术，整个极限问题的求解将在维数很低的自平衡应力场子空间中进行迭代。具体过程如下：首先根据4.1 节的计算得到结构的弹性极限载荷，之后在此基础上施加一个载荷增量△χ1，结合4.2 节的步骤得到一组自平衡应力场基矢量，并通过复合形法[12]对该数学规划问题进行求解，得到近似解χ1。由于此时规划式中的未知数只有K+1 个（即待求的极限载荷系数和基矢量对应的K个待定系数），因此计算规模得到很大程度的减小。χ1的精度一般来说并不高，还需进行第2 次数学规划问题的计算。以第1 次数学规划问题得到的平衡应力场为基础，施加第2 个载荷增量△χ2，同时将求解第1 次数学规划问题得到的自平衡应力场作为第K个基矢量，由此得到一组新的自平衡应力场基矢量，并再次通过复合形法计算得到近似解χ2。如此重复进行，直到满足收敛条件：

则计算终止，式中vol 表示误差容限。

3 平衡应力场构造

3.1 虚拟弹性应力场

假设在轴对称面S上分布着体力b= [br bz]T，面力边界Γt作用有给定的面力，位移边界Γu有已知位移，则通过加权余量法可得到轴对称弹性力学问题控制方程的等效弱形式为：

式中，r代表轴对称半径，σ表示应力向量：

根据有限元法插值方法，将求解域内的位移场函数代入式（15）中，并注意到位移变分δu具有任意性，由此进一步推导出轴对称结构线弹性问题的控制方程：

式中u代表求解域中所有节点的位移向量，K和f分别表示总体刚度矩阵和总体载荷向量。

3.2 自平衡应力场

约束条件式（10）表明，求解极限下限问题的关键在于找到一个能让载荷系数χ取最大值的自平衡应力场。为了让迭代得到良好的收敛效果，本文利用弹塑性增量分析中平衡迭代的结果构造自平衡应力基矢量，具体过程如下：

以第m次迭代为例，根据上一次迭代得到的载荷系数χ（m-1）和自平衡应力场ρm-1i，则各Gauss 积分点总应力为：

现在此基础上增加一个载荷增量△χm使结构进一步屈服。为提高计算效率，本文选用修正的Newton-Raphson 迭代法[11]进行计算。则第1 次和第p迭代的平衡表达式为：

其中F代表基准载荷作用时对应的节点载荷向量，De为轴对称问题的弹性矩阵。用（20）减去（19），得到：

式中ρp是一个与外载荷有关的自平衡应力场。将每一次迭代得到的ρp作为一个基矢量，就形成了自平衡应力场的搜索子空间。

4 算例验证

4.1 厚壁圆筒

一个理想弹塑性厚壁圆筒结构，受到均布内压P的作用，分别用r和R表示厚壁圆筒的内径和外径，其计算模型如图1 所示。根据塑性力学理论，厚壁圆筒极限载荷的解析解计算表达式可写为[13]：

图1 受内压厚壁圆筒

本文对R/r在2.0 ～4.0 的一组厚壁圆筒进行极限分析。对于不同的内外径之比，设置不同的单元数量进行计算。图2 所示为当R= 2r时有限元网格的划分情况。为验证本文算法的有效性，将采用本文算法得到的数值解与解析解进行对比，具体情况如图3所示。从图中可以看出，数值解和解析解形成的曲线吻合很好，说明了所提算法的有效性。

图2 R/r=2 时厚壁圆筒单元布置

图3 本文数值解与解析解的比较

4.2 厚壁球壳

对一个内径为r，外径为R的厚壁球壳进行分析，球壳内部受到均布载荷P的作用。图4 为R=1.3r时厚壁球壳的计算模型，则球壳结构的极限载荷解析解计算表达式为[13]：

图4 受内压厚壁球壳

为进一步验证所提方法的有效性，算例中计算了外径与内径之比为1.1～1.3 之间的一组厚壁球壳，计算结果如表1 所示。表1 结果显示，采用本文算法得到的数值解和理论计算得到的解析解之间的误差比较小，都在1%以内。为了探究基矢量数目对计算精度的影响，算例中以R/r= 1.3 时的厚壁球壳为研究对象，分别设置不同数目的基矢量进行计算，结果如表2 所示。表2 结果显示，不同基矢量数目得到的计算结果均稳定在0.52 附近，这充分说明了本文算法的合理性和稳定性。

表1 厚壁球壳极限下限载荷系数的计算结果

表2 R=1.3r 时厚壁球壳不同基矢量数目的计算结果

5 结语

本文基于塑性极限分析理论，将有限元法应用于求解轴对称结构的极限下限问题中，并建立了相应算法。通过将整个下限分析转化为一系列未知数较少的非线性数学规划子问题，解决了维数障碍问题。每个子问题中通过弹塑性增量法迭代的结果构造出相应的自平衡应力场，并采用复合形法直接进行求解。算例的计算结果表明，本文提出的算法具有数值稳定性好、计算精度高的优点。