殷毓基，梁 珂，孫 秦

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

0 引言

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

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

本文采用這種新型的非線性有限元降階方法對航天結(jié)構(gòu)工程中常用的網(wǎng)格加筋筒殼結(jié)構(gòu)進(jìn)行分析計(jì)算，重點(diǎn)研究其在軸向壓縮載荷作用下的穩(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)點(diǎn)處建立有限元全模型規(guī)模的非線性平衡方程以及相應(yīng)的降階模型。定義結(jié)構(gòu)的3個(gè)位移場變量u0、u、q和3個(gè)載荷系數(shù)變量λ0、Δλ、λ。它們之間滿足以下的兩個(gè)關(guān)系式

q=u0+u

(1)

λ=λ0+Δλ

(2)

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

K(q)q=λfext

(3)

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

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

Ktw=zKgw

(4)

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

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

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

(5)

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

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

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為攝動(dòng)參數(shù)向量a的各分量；Q、S、P

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

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

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

(8)

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

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

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

表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)點(diǎn)處的前m階密集屈曲模態(tài)，然后用這些模態(tài)來構(gòu)造攝動(dòng)載荷，進(jìn)而建立當(dāng)前狀態(tài)下的結(jié)構(gòu)降階模型。為此，在加載端均勻施加軸向壓載荷，進(jìn)行結(jié)構(gòu)特征值屈曲分析，得到該結(jié)構(gòu)的臨界屈曲載荷為13795kN，前5階屈曲載荷如表2所示。

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

2.2 非線性屈曲分析

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

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

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

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

3 結(jié)論

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

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

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