基于非线性有限元降阶方法的网格加筋筒壳结构屈曲承载特性研究

殷毓基，梁 珂，孙 秦

(西北工业大学航空学院, 西安 710072)

0 引言

作为航空航天工程的典型结构构型，薄壁结构具有质量小、承载效率高及可设计性强等技术优势，其科学问题的复杂性及工业应用的安全可靠性历来备受学术及工业界的关注，被认为是未来航空航天结构低成本技术发展的主要方向之一。由于飞行器发动机的强大推力，薄壁结构在服役过程中会受到沿轴向的各种复杂载荷作用，极易发生屈曲失稳。薄壁结构屈曲时不仅会发生较大的变形，而且会降低结构的离面刚度，呈现出明显的几何非线性效应，增大工程应用风险，因此结构的非线性屈曲响应分析是轻质薄壁结构设计的重要技术环节，是确保飞机、火箭以及人造卫星等大型航空航天器服役安全的关键[1]。

为计及几何非线性效应，当前计算结构屈曲响应的常规方法是基于Newton增量迭代技术的非线性有限元法。潘光等[2]采用非线性有限元方法分析了复合材料圆柱壳的屈曲性能。该方法需要对有限元离散后的大规模非线性方程组进行迭代求解，存在求解效率低、计算规模大等问题;该计算量问题在需要开展结构重分析的结构缺陷敏感度分析中显得尤为突出。与上述基于有限元全模型的分析相对应，旨在缩减有限元模型中自由度数目的降阶方法近些年引起不少学者的关注，其中最为活跃的当属基于Koiter[3]摄动理论提出的一系列Koiter降阶方法[4-7]。该方法在屈曲分支点处采用屈曲模态对结构的初始后屈曲平衡路径进行摄动展开，从而将原先大规模的非线性有限元平衡方程组缩减为关于若干摄动参数的小规模非线性方程组，即降阶模型，进而求解。近两年，梁珂等[8-9]人将Koiter摄动理论与Newton法的增量迭代思想相结合，并经有限元实现后提出了一种全新的非线性有限元降阶算法。如图1所示，该方法改进传统的Koiter摄动理论，使其能够在结构非线性平衡路径上的任意平衡状态点处使用，并在每个载荷步中采用非线性预测+迭代修正策略，实现对结构静力屈曲非线性平衡路径的高效、高精度跟踪。

本文采用这种新型的非线性有限元降阶方法对航天结构工程中常用的网格加筋筒壳结构进行分析计算，重点研究其在轴向压缩载荷作用下的稳定性承载特性，以及结构缺陷因素对结构屈曲载荷的影响。

图1 非线性有限元降阶方法的基本思想Fig.1 Basic idea of the nonlinear finite element reduced-order method

1 非线性有限元降阶方法

非线性有限元降阶方法的基本算法步骤如下:

首先，在结构的已知平衡状态点处建立有限元全模型规模的非线性平衡方程以及相应的降阶模型。定义结构的3个位移场变量u0、u、q和3个载荷系数变量λ0、Δλ、λ。它们之间满足以下的两个关系式

q=u0+u

(1)

λ=λ0+Δλ

(2)

式(1)(2)中，(u0,λ0)为结构平衡路径上已知平衡状态点,(q,λ)为已知平衡点附近的一个未知平衡状态点,(u,Δλ)为这两个平衡状态点之间的增量。在结构的已知平衡状态点(从未变形状态点处开始)上建立离散化的非线性有限元平衡方程

K(q)q=λfext

(3)

式中，q为结构的变形位移场向量,K(q)为与位移场相关的结构的非线性有限元刚度矩阵,fext为外载荷向量,λ为外载荷向量系数。当前传统方法是直接求解该非线性方程组获得结构的载荷位移(q-λ)曲线。该有限元全模型K的阶数较大，传统方法在增量迭代求解过程中需要反复分解大规模矩阵,导致计算量大、分析时间长。

结构的降阶模型是基于改进的Koiter摄动展开理论建立的。在结构的已知平衡点上求解结构的线性屈曲特征值问题

Ktw=zKgw

(4)

式中，Kt和Kg分别为结构的切线刚度矩阵和几何刚度矩阵；w为屈曲模态；z为特征值；提取结构的密集屈曲模态，个数为m，来构建降阶模型。将Kg与特征向量w的乘积定义为摄动载荷。

然后，将结构平衡方程(3)在已知平衡状态点处关于位移场u展开至3阶项，可得

U′(u)+U″(u,u)+U″″(u,u,u)=Fφ

(5)

式中，U为结构的应变能，U右上方小撇的个数表示对应变能求导的阶数；F为载荷向量矩阵，它的第1列由外载荷向量构成，其余各列由摄动载荷向量构成；φ为载荷系数向量。

接下来，将位移场u和载荷系数向量φ关于摄动参数向量a分别展开至3阶项，可得

u=uiai+uijaiaj+uijkaiajak

(6)

φ=Q(a)+S(a)+P(a)

(7)

式中，i、j、k= 1, 2, …,m+1；ui、uij、uijk分别为结构的1阶、2阶和3阶位移场；ai、aj、ak为摄动参数向量a的各分量；Q、S、P

分别为待求解的1次、2次和3次算子。

最后，将式(6)和式(7)分别代入式(5)的左右两端，令摄动参数a的各次幂的系数为0，即可获得式(7)中各算子的表达式。此时，结构的降阶模型，即式(8)即可建立

Q(a)+S(a)+P(a)=φ

(8)

该降阶模型本质上是一个(1+m)阶的非线性方程组，其阶数一般<10。

采用弧长法求解可得到载荷系数λ和摄动参数a的关系，再代入位移展开式(6)就可获得结构的非线性响应，即λ-u曲线。

以上阐述了单个摄动展开步中降阶模型的建立和求解过程，对于某些变形较大而非线性效应不是非常明显的结构，往往一个摄动展开步(在结构的未变形状态处进行一次摄动展开)就可以跟踪到结构的屈曲临界载荷，并获得结构初始后屈曲阶段的响应。但是降阶模型相比有限元全模型而言终归还是一个近似模型，在单个摄动展开步中，当结构变形增大时它的计算误差也将随之增大。为了让该算法实现对结构非线性平衡路径的自动跟踪，仍需要给出每个摄动展开步的有效范围，即确定下一个摄动展开步的起始点。具体过程如下：首先，在求解当前摄动展开步中降阶模型的同时，通过时刻计算结构的残余力来判断当位移达到多大时该降阶模型已经丧失求解精度。接下来，在修正步中将降阶模型已经偏离了的解拉回到结构真实的平衡路径上，从而获得真实解。最后，以该真实解作为下一个摄动展开步的已知起始点建立新的降阶模型并进行求解。这个过程可自动反复进行直到获得想要的结构响应为止。

2 轴压网格加筋筒壳结构的屈曲分析

网格加筋筒壳结构的具体几何尺寸见表1，几何模型如图2所示。该结构由筒壳和筋条两部分组成，筋条截面形状为矩形， 含有0°、+60°、-60°这3种方向的筋条各80根。筒壳和筋条的材料属性相同，即弹性模量为68246MPa，泊松比为0.3，加筋筒壳基于等效刚度模型，采用板壳单元进行模拟。加载端约束除了轴向以外的全部自由度，约束端固支。

表1 网格加筋筒壳结构尺寸参数(单位：mm)

图2 金属加筋筒几何模型Fig.2 Model of the metal grid-stiffened cylinder

2.1 特征值屈曲分析

由第1节非线性有限元降阶方法的分析步骤可知，首先需要获得结构在当前已知平衡状态点处的前m阶密集屈曲模态，然后用这些模态来构造摄动载荷，进而建立当前状态下的结构降阶模型。为此，在加载端均匀施加轴向压载荷，进行结构特征值屈曲分析，得到该结构的临界屈曲载荷为13795kN，前5阶屈曲载荷如表2所示。

表2 特征值分析得到的前5阶屈曲载荷值(单位：kN)

2.2 非线性屈曲分析

依照第1节中介绍的非线性有限元降阶方法分析步骤，在未变形状态点处建立降阶模型时需基于该结构的前5阶密集屈曲模态，得到的降阶模型总自由度数为6。通过求解降阶模型可以获得该载荷步的非线性预测解，当预测解精度不足时，采用Newton迭代格式进行修正，随后在新的平衡状态点处重新建立降阶模型，开始下一个载荷步求解，最终得到结构压缩位移随加载变化曲线，如图3所示。图3中的两条结构响应曲线分别为采用基于常规结构非线性有限元方法的ABAQUS软件和本文介绍的非线性有限元降阶方法计算获得。通过比较图3中的曲线可知，采用两种方法获得的承载响应曲线在极值点处吻合较好，提取曲线最高点处数据得到结构的非线性屈曲载荷为12563.3kN，曲线的微小差异主要归结于两种方法所采用的板壳单元构造原理不同。获得图3中的结构响应曲线，常规结构非线性有限元方法需要计算17个载荷步，而非线性有限元降阶方法仅需要计算5个载荷步。非线性有限元降阶方法在每个载荷步中建立降阶模型需要求解1个有限元模型规模的线性代数方程组，对预测解进行修正时需要求解2个线性代数方程组，因此每个载荷步需要求解3个线性代数方程组。通过控制求解参数常规结构非线性有限元方法的每个载荷步也大约需要求解3个线性代数方程组，这两种方法在每个载荷步中的计算量相当，但非线性降阶方法需要更少的载荷步，因而可见降阶方法因具有更大的载荷步长而在非线性计算效率方面的优势明显，计算量仅为常规方法的1/3。

图3 结构非线性分析获得承载响应曲线Fig.3 Loading response curve obtained from the structural nonlinear analysis

2.3 结构几何形状缺陷分析

采用侧向小扰动载荷所形成的结构表面形变来模拟结构的初始几何形状缺陷。小扰动载荷施加的位置在筒壳的轴向中点处，并垂直壳表面指向圆筒中心。如图4所示，几何形状缺陷的大小可由施加扰动载荷的大小来控制。改变扰动载荷的大小并进行结构非线性分析，得到各个扰动载荷下结构的屈曲载荷，并以归一化载荷(载荷值/特征值线性屈曲载荷值)的形式来表示。

采用本文的非线性有限元降阶方法，针对14种不同的扰动载荷大小进行结构分析，获得图5中的屈曲载荷随扰动载荷变化曲线。由图5可以看出，结构屈曲载荷先是随扰动载荷的增大而降低，但扰动载荷达到一定程度后曲线逐渐放平，结构的屈曲载荷趋于一个定值，此时归一化载荷大约为0.75。需要注意的是，采用常规非线性有限元方法(如ABAQUS)，对于不同的扰动载荷工况需要重新计算大规模的结构非线性问题，因而图5的14种扰动载荷工况的总计算量为14T，这里T为单次结构非线性分析的计算量。而本文采用的非线性有限元降阶方法在不同扰动载荷下不需要再重新建立降阶模型，只用重新求解一遍小规模(7自由度)的降阶模型即可，而该部分的计算量是可以忽略的，因而对14种扰动载荷工况的总计算量仅相当于一次结构非线性分析，约为T/3(单次结构非线性分析计算量仅为常规方法的1/3)。由此可见，非线性有限元降阶方法在对含缺陷结构进行重分析时的计算效率优势更加明显。

图4 扰动载荷施加位置示意图Fig.4 The applied position of the disturbing load

图5 不同扰动载荷下的结构屈曲载荷Fig.5 Buckling loads with different disturbing load values

3 结论

本文采用结构非线性有限元降阶方法对金属网格加筋筒壳结构的稳定性承载特性进行了细致分析。具体工作总结如下：

1)首先对轴压载荷作用下的网格加筋筒壳结构进行了线性特征值屈曲分析，获得建立降阶模型所需要的密集屈曲模态，随即采用非线性有限元降阶方法计算了结构的非线性屈曲载荷，通过与常规非线性有限元方法比较，验证了降阶方法的分析精度，并展示了其高效的计算效率，即其对结构非线性屈曲分析的计算量仅为常规方法的1/3；

2)对包含几何形状缺陷的加筋筒壳进行了非线性屈曲分析，通过与常规非线性有限元方法比较，发现非线性有限元降阶方法在对含缺陷结构进行重分析时的计算效率优势更加明显。同时研究了几何形状缺陷对结构稳定性承载的影响规律。从结果曲线来看，承载力变化呈现出先随缺陷增大而下降，但缺陷达到一定程度后承载力变化又趋于变缓，即变化曲线出现平台现象。这个承载力变化的下限值对结构的安全可靠性设计有着十分重要的意义。