殷毓基，梁 珂，孫 秦

(西北工業(yè)大學(xué)航空學(xué)院, 西安 710072)

0 引言

作為航空航天工程的典型結(jié)構(gòu)構(gòu)型，薄壁結(jié)構(gòu)具有質(zhì)量小、承載效率高及可設(shè)計性強等技術(shù)優(yōu)勢，其科學(xué)問題的復(fù)雜性及工業(yè)應(yīng)用的安全可靠性歷來備受學(xué)術(shù)及工業(yè)界的關(guān)注，被認為是未來航空航天結(jié)構(gòu)低成本技術(shù)發(fā)展的主要方向之一。由于飛行器發(fā)動機的強大推力，薄壁結(jié)構(gòu)在服役過程中會受到沿軸向的各種復(fù)雜載荷作用，極易發(fā)生屈曲失穩(wěn)。薄壁結(jié)構(gòu)屈曲時不僅會發(fā)生較大的變形，而且會降低結(jié)構(gòu)的離面剛度，呈現(xiàn)出明顯的幾何非線性效應(yīng)，增大工程應(yīng)用風(fēng)險，因此結(jié)構(gòu)的非線性屈曲響應(yīng)分析是輕質(zhì)薄壁結(jié)構(gòu)設(shè)計的重要技術(shù)環(huán)節(jié)，是確保飛機、火箭以及人造衛(wèi)星等大型航空航天器服役安全的關(guān)鍵[1]。

為計及幾何非線性效應(yīng)，當(dāng)前計算結(jié)構(gòu)屈曲響應(yīng)的常規(guī)方法是基于Newton增量迭代技術(shù)的非線性有限元法。潘光等[2]采用非線性有限元方法分析了復(fù)合材料圓柱殼的屈曲性能。該方法需要對有限元離散后的大規(guī)模非線性方程組進行迭代求解，存在求解效率低、計算規(guī)模大等問題;該計算量問題在需要開展結(jié)構(gòu)重分析的結(jié)構(gòu)缺陷敏感度分析中顯得尤為突出。與上述基于有限元全模型的分析相對應(yīng)，旨在縮減有限元模型中自由度數(shù)目的降階方法近些年引起不少學(xué)者的關(guān)注，其中最為活躍的當(dāng)屬基于Koiter[3]攝動理論提出的一系列Koiter降階方法[4-7]。該方法在屈曲分支點處采用屈曲模態(tài)對結(jié)構(gòu)的初始后屈曲平衡路徑進行攝動展開，從而將原先大規(guī)模的非線性有限元平衡方程組縮減為關(guān)于若干攝動參數(shù)的小規(guī)模非線性方程組，即降階模型，進而求解。近兩年，梁珂等[8-9]人將Koiter攝動理論與Newton法的增量迭代思想相結(jié)合，并經(jīng)有限元實現(xiàn)后提出了一種全新的非線性有限元降階算法。如圖1所示，該方法改進傳統(tǒng)的Koiter攝動理論，使其能夠在結(jié)構(gòu)非線性平衡路徑上的任意平衡狀態(tài)點處使用，并在每個載荷步中采用非線性預(yù)測+迭代修正策略，實現(xiàn)對結(jié)構(gòu)靜力屈曲非線性平衡路徑的高效、高精度跟蹤。

本文采用這種新型的非線性有限元降階方法對航天結(jié)構(gòu)工程中常用的網(wǎng)格加筋筒殼結(jié)構(gòu)進行分析計算，重點研究其在軸向壓縮載荷作用下的穩(wěn)定性承載特性，以及結(jié)構(gòu)缺陷因素對結(jié)構(gòu)屈曲載荷的影響。

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

1 非線性有限元降階方法

非線性有限元降階方法的基本算法步驟如下:

首先，在結(jié)構(gòu)的已知平衡狀態(tài)點處建立有限元全模型規(guī)模的非線性平衡方程以及相應(yīng)的降階模型。定義結(jié)構(gòu)的3個位移場變量u0、u、q和3個載荷系數(shù)變量λ0、Δλ、λ。它們之間滿足以下的兩個關(guān)系式

q=u0+u

(1)

λ=λ0+Δλ

(2)

式(1)(2)中，(u0,λ0)為結(jié)構(gòu)平衡路徑上已知平衡狀態(tài)點,(q,λ)為已知平衡點附近的一個未知平衡狀態(tài)點,(u,Δλ)為這兩個平衡狀態(tài)點之間的增量。在結(jié)構(gòu)的已知平衡狀態(tài)點(從未變形狀態(tài)點處開始)上建立離散化的非線性有限元平衡方程

K(q)q=λfext

(3)

式中，q為結(jié)構(gòu)的變形位移場向量,K(q)為與位移場相關(guān)的結(jié)構(gòu)的非線性有限元剛度矩陣,fext為外載荷向量,λ為外載荷向量系數(shù)。當(dāng)前傳統(tǒng)方法是直接求解該非線性方程組獲得結(jié)構(gòu)的載荷位移(q-λ)曲線。該有限元全模型K的階數(shù)較大，傳統(tǒng)方法在增量迭代求解過程中需要反復(fù)分解大規(guī)模矩陣,導(dǎo)致計算量大、分析時間長。

結(jié)構(gòu)的降階模型是基于改進的Koiter攝動展開理論建立的。在結(jié)構(gòu)的已知平衡點上求解結(jié)構(gòu)的線性屈曲特征值問題

Ktw=zKgw

(4)

式中，Kt和Kg分別為結(jié)構(gòu)的切線剛度矩陣和幾何剛度矩陣；w為屈曲模態(tài)；z為特征值；提取結(jié)構(gòu)的密集屈曲模態(tài)，個數(shù)為m，來構(gòu)建降階模型。將Kg與特征向量w的乘積定義為攝動載荷。

然后，將結(jié)構(gòu)平衡方程(3)在已知平衡狀態(tài)點處關(guān)于位移場u展開至3階項，可得

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

(5)

式中，U為結(jié)構(gòu)的應(yīng)變能，U右上方小撇的個數(shù)表示對應(yīng)變能求導(dǎo)的階數(shù)；F為載荷向量矩陣，它的第1列由外載荷向量構(gòu)成，其余各列由攝動載荷向量構(gòu)成；φ為載荷系數(shù)向量。

接下來，將位移場u和載荷系數(shù)向量φ關(guān)于攝動參數(shù)向量a分別展開至3階項，可得

u=uiai+uijaiaj+uijkaiajak

(6)

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

(7)

式中，i、j、k= 1, 2, …,m+1；ui、uij、uijk分別為結(jié)構(gòu)的1階、2階和3階位移場；ai、aj、ak為攝動參數(shù)向量a的各分量；Q、S、P

分別為待求解的1次、2次和3次算子。

最后，將式(6)和式(7)分別代入式(5)的左右兩端，令攝動參數(shù)a的各次冪的系數(shù)為0，即可獲得式(7)中各算子的表達式。此時，結(jié)構(gòu)的降階模型，即式(8)即可建立

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

(8)

該降階模型本質(zhì)上是一個(1+m)階的非線性方程組，其階數(shù)一般<10。

采用弧長法求解可得到載荷系數(shù)λ和攝動參數(shù)a的關(guān)系，再代入位移展開式(6)就可獲得結(jié)構(gòu)的非線性響應(yīng)，即λ-u曲線。

以上闡述了單個攝動展開步中降階模型的建立和求解過程，對于某些變形較大而非線性效應(yīng)不是非常明顯的結(jié)構(gòu)，往往一個攝動展開步(在結(jié)構(gòu)的未變形狀態(tài)處進行一次攝動展開)就可以跟蹤到結(jié)構(gòu)的屈曲臨界載荷，并獲得結(jié)構(gòu)初始后屈曲階段的響應(yīng)。但是降階模型相比有限元全模型而言終歸還是一個近似模型，在單個攝動展開步中，當(dāng)結(jié)構(gòu)變形增大時它的計算誤差也將隨之增大。為了讓該算法實現(xiàn)對結(jié)構(gòu)非線性平衡路徑的自動跟蹤，仍需要給出每個攝動展開步的有效范圍，即確定下一個攝動展開步的起始點。具體過程如下：首先，在求解當(dāng)前攝動展開步中降階模型的同時，通過時刻計算結(jié)構(gòu)的殘余力來判斷當(dāng)位移達到多大時該降階模型已經(jīng)喪失求解精度。接下來，在修正步中將降階模型已經(jīng)偏離了的解拉回到結(jié)構(gòu)真實的平衡路徑上，從而獲得真實解。最后，以該真實解作為下一個攝動展開步的已知起始點建立新的降階模型并進行求解。這個過程可自動反復(fù)進行直到獲得想要的結(jié)構(gòu)響應(yīng)為止。

2 軸壓網(wǎng)格加筋筒殼結(jié)構(gòu)的屈曲分析

網(wǎng)格加筋筒殼結(jié)構(gòu)的具體幾何尺寸見表1，幾何模型如圖2所示。該結(jié)構(gòu)由筒殼和筋條兩部分組成，筋條截面形狀為矩形， 含有0°、+60°、-60°這3種方向的筋條各80根。筒殼和筋條的材料屬性相同，即彈性模量為68246MPa，泊松比為0.3，加筋筒殼基于等效剛度模型，采用板殼單元進行模擬。加載端約束除了軸向以外的全部自由度，約束端固支。

表1 網(wǎng)格加筋筒殼結(jié)構(gòu)尺寸參數(shù)(單位：mm)

圖2 金屬加筋筒幾何模型Fig.2 Model of the metal grid-stiffened cylinder

2.1 特征值屈曲分析

由第1節(jié)非線性有限元降階方法的分析步驟可知，首先需要獲得結(jié)構(gòu)在當(dāng)前已知平衡狀態(tài)點處的前m階密集屈曲模態(tài)，然后用這些模態(tài)來構(gòu)造攝動載荷，進而建立當(dāng)前狀態(tài)下的結(jié)構(gòu)降階模型。為此，在加載端均勻施加軸向壓載荷，進行結(jié)構(gòu)特征值屈曲分析，得到該結(jié)構(gòu)的臨界屈曲載荷為13795kN，前5階屈曲載荷如表2所示。

表2 特征值分析得到的前5階屈曲載荷值(單位：kN)

2.2 非線性屈曲分析

依照第1節(jié)中介紹的非線性有限元降階方法分析步驟，在未變形狀態(tài)點處建立降階模型時需基于該結(jié)構(gòu)的前5階密集屈曲模態(tài)，得到的降階模型總自由度數(shù)為6。通過求解降階模型可以獲得該載荷步的非線性預(yù)測解，當(dāng)預(yù)測解精度不足時，采用Newton迭代格式進行修正，隨后在新的平衡狀態(tài)點處重新建立降階模型，開始下一個載荷步求解，最終得到結(jié)構(gòu)壓縮位移隨加載變化曲線，如圖3所示。圖3中的兩條結(jié)構(gòu)響應(yīng)曲線分別為采用基于常規(guī)結(jié)構(gòu)非線性有限元方法的ABAQUS軟件和本文介紹的非線性有限元降階方法計算獲得。通過比較圖3中的曲線可知，采用兩種方法獲得的承載響應(yīng)曲線在極值點處吻合較好，提取曲線最高點處數(shù)據(jù)得到結(jié)構(gòu)的非線性屈曲載荷為12563.3kN，曲線的微小差異主要歸結(jié)于兩種方法所采用的板殼單元構(gòu)造原理不同。獲得圖3中的結(jié)構(gòu)響應(yīng)曲線，常規(guī)結(jié)構(gòu)非線性有限元方法需要計算17個載荷步，而非線性有限元降階方法僅需要計算5個載荷步。非線性有限元降階方法在每個載荷步中建立降階模型需要求解1個有限元模型規(guī)模的線性代數(shù)方程組，對預(yù)測解進行修正時需要求解2個線性代數(shù)方程組，因此每個載荷步需要求解3個線性代數(shù)方程組。通過控制求解參數(shù)常規(guī)結(jié)構(gòu)非線性有限元方法的每個載荷步也大約需要求解3個線性代數(shù)方程組，這兩種方法在每個載荷步中的計算量相當(dāng)，但非線性降階方法需要更少的載荷步，因而可見降階方法因具有更大的載荷步長而在非線性計算效率方面的優(yōu)勢明顯，計算量僅為常規(guī)方法的1/3。

圖3 結(jié)構(gòu)非線性分析獲得承載響應(yīng)曲線Fig.3 Loading response curve obtained from the structural nonlinear analysis

2.3 結(jié)構(gòu)幾何形狀缺陷分析

采用側(cè)向小擾動載荷所形成的結(jié)構(gòu)表面形變來模擬結(jié)構(gòu)的初始幾何形狀缺陷。小擾動載荷施加的位置在筒殼的軸向中點處，并垂直殼表面指向圓筒中心。如圖4所示，幾何形狀缺陷的大小可由施加擾動載荷的大小來控制。改變擾動載荷的大小并進行結(jié)構(gòu)非線性分析，得到各個擾動載荷下結(jié)構(gòu)的屈曲載荷，并以歸一化載荷(載荷值/特征值線性屈曲載荷值)的形式來表示。

采用本文的非線性有限元降階方法，針對14種不同的擾動載荷大小進行結(jié)構(gòu)分析，獲得圖5中的屈曲載荷隨擾動載荷變化曲線。由圖5可以看出，結(jié)構(gòu)屈曲載荷先是隨擾動載荷的增大而降低，但擾動載荷達到一定程度后曲線逐漸放平，結(jié)構(gòu)的屈曲載荷趨于一個定值，此時歸一化載荷大約為0.75。需要注意的是，采用常規(guī)非線性有限元方法(如ABAQUS)，對于不同的擾動載荷工況需要重新計算大規(guī)模的結(jié)構(gòu)非線性問題，因而圖5的14種擾動載荷工況的總計算量為14T，這里T為單次結(jié)構(gòu)非線性分析的計算量。而本文采用的非線性有限元降階方法在不同擾動載荷下不需要再重新建立降階模型，只用重新求解一遍小規(guī)模(7自由度)的降階模型即可，而該部分的計算量是可以忽略的，因而對14種擾動載荷工況的總計算量僅相當(dāng)于一次結(jié)構(gòu)非線性分析，約為T/3(單次結(jié)構(gòu)非線性分析計算量僅為常規(guī)方法的1/3)。由此可見，非線性有限元降階方法在對含缺陷結(jié)構(gòu)進行重分析時的計算效率優(yōu)勢更加明顯。

圖4 擾動載荷施加位置示意圖Fig.4 The applied position of the disturbing load

圖5 不同擾動載荷下的結(jié)構(gòu)屈曲載荷Fig.5 Buckling loads with different disturbing load values

3 結(jié)論

本文采用結(jié)構(gòu)非線性有限元降階方法對金屬網(wǎng)格加筋筒殼結(jié)構(gòu)的穩(wěn)定性承載特性進行了細致分析。具體工作總結(jié)如下：

1)首先對軸壓載荷作用下的網(wǎng)格加筋筒殼結(jié)構(gòu)進行了線性特征值屈曲分析，獲得建立降階模型所需要的密集屈曲模態(tài)，隨即采用非線性有限元降階方法計算了結(jié)構(gòu)的非線性屈曲載荷，通過與常規(guī)非線性有限元方法比較，驗證了降階方法的分析精度，并展示了其高效的計算效率，即其對結(jié)構(gòu)非線性屈曲分析的計算量僅為常規(guī)方法的1/3；

2)對包含幾何形狀缺陷的加筋筒殼進行了非線性屈曲分析，通過與常規(guī)非線性有限元方法比較，發(fā)現(xiàn)非線性有限元降階方法在對含缺陷結(jié)構(gòu)進行重分析時的計算效率優(yōu)勢更加明顯。同時研究了幾何形狀缺陷對結(jié)構(gòu)穩(wěn)定性承載的影響規(guī)律。從結(jié)果曲線來看，承載力變化呈現(xiàn)出先隨缺陷增大而下降，但缺陷達到一定程度后承載力變化又趨于變緩，即變化曲線出現(xiàn)平臺現(xiàn)象。這個承載力變化的下限值對結(jié)構(gòu)的安全可靠性設(shè)計有著十分重要的意義。