何歡，宋大鵬，張晨凱, 2，陳國平, 2

1.南京航空航天大學(xué) 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室,南京 210016 2. 南京航空航天大學(xué) 振動(dòng)工程研究所,南京 210016

機(jī)翼作為飛機(jī)的重要承力部件，對飛機(jī)的飛行性能有著重要的影響[1]。在結(jié)構(gòu)設(shè)計(jì)與強(qiáng)度分析階段，對機(jī)翼進(jìn)行有限元建模及分析是不可或缺的一步。目前機(jī)翼最為常見的建模方法是將機(jī)翼簡化為“桿板組合結(jié)構(gòu)”[2]，但該建模方法工作量較大，用于靜力學(xué)分析尚可，動(dòng)力學(xué)分析時(shí)就會(huì)顯得計(jì)算量過大。在機(jī)翼結(jié)構(gòu)初級設(shè)計(jì)以及結(jié)構(gòu)優(yōu)化階段，機(jī)翼各構(gòu)件的參數(shù)會(huì)不斷調(diào)整變化，采用“桿板組合結(jié)構(gòu)”建模方法就會(huì)顯得較為繁雜，效率不高[3]。為了避免分析過于復(fù)雜，特別是在動(dòng)力學(xué)分析中，只需要能夠反映結(jié)構(gòu)總體性能的簡化模型即可，因而需要尋求一種更為簡化的模型。

建立等效有限元模型是解決上述問題的一種途徑，如等效梁模型[4-6]和等效板模型[7-9]。然而，等效有限元模型的建立需要結(jié)合工程經(jīng)驗(yàn)，而且精度也受到一定的限制。此外，同時(shí)考慮靜力與動(dòng)力學(xué)分析的等效有限元模型缺乏全面的研究[10]。因此為了簡單、有效地進(jìn)行結(jié)構(gòu)特性分析，建模時(shí)盡可能減少結(jié)構(gòu)自由度，建立更加簡潔的模型以及分析方法就顯得很有必要。考慮到機(jī)翼以及機(jī)翼各主要部件(蒙皮、長桁、翼梁等)通常都是細(xì)長結(jié)構(gòu)的特點(diǎn)，如果能夠建立一種梁模型，使得機(jī)翼中各細(xì)長部件均采用一維梁單元進(jìn)行建模，那么模型將十分簡潔，結(jié)構(gòu)自由度也將大大降低，進(jìn)而給機(jī)翼的有限元建模與分析帶來便利。

目前最為常見的梁模型有Euler-Bernoulli梁模型和Timoshenko梁模型。Euler-Bernoulli梁忽略了橫向剪力和橫向正應(yīng)變的影響，適用于剪切變形影響較小的細(xì)長梁[11]，Timoshenko梁假設(shè)剪應(yīng)變在給定橫截面上為常值，而切應(yīng)力卻不是常數(shù)，與彈性力學(xué)不符。對于高跨比較大的深梁而言，其切應(yīng)力沿截面變化也不再滿足一階剪切變形理論，故Timoshenko梁理論也不再適用[12-16]。在一些工程分析中，如機(jī)翼等在復(fù)雜載荷作用下，高階剪切變形、翹曲變形較大，而上述梁理論均無法反映截面面內(nèi)變形，因此存在計(jì)算精度不足的問題。在一些商用軟件中如需考慮高階剪切變形的影響，或者處理復(fù)雜截面形狀結(jié)構(gòu)時(shí)，則推薦使用實(shí)體單元，而實(shí)體單元建模又會(huì)帶來結(jié)構(gòu)自由度大幅增加、效率偏低的問題。

為此，本文引入若干個(gè)截面特征點(diǎn)，以截面特征點(diǎn)位移作為未知量，引入拉格朗日插值函數(shù)對梁軸向坐標(biāo)上的任意截面上的特征點(diǎn)的位移進(jìn)行插值，描述梁截面的復(fù)雜變形模式，建立一種可用于任意截面梁的有限元模型，結(jié)合虛功原理推導(dǎo)了截面插值梁單元的剛度矩陣和質(zhì)量矩陣。最后，利用截面插值梁模型對機(jī)翼各部件進(jìn)行建模，并展開典型工況下的靜力學(xué)分析以及動(dòng)力學(xué)分析(模態(tài)分析、顫振分析、動(dòng)響應(yīng)分析)，通過與Nastran實(shí)體單元計(jì)算結(jié)果對比驗(yàn)證該模型的可靠性。

1 構(gòu)造梁位移場

1.1 拉格朗日插值多項(xiàng)式介紹

二元拉格朗日插值多項(xiàng)式一個(gè)重要的特點(diǎn)是其在自身點(diǎn)處值等于1，而在其他點(diǎn)處值等于0，且所有插值函數(shù)之和恒為1[17]。根據(jù)插值點(diǎn)數(shù)選擇的不同，可以得到不同精度的插值函數(shù)，與有限元理論中的二維拉格朗日單元插值函數(shù)相同。三點(diǎn)插值(L3)、四點(diǎn)插值(L4)和九點(diǎn)插值(L9)的拉格朗日多項(xiàng)式分別對應(yīng)于傳統(tǒng)有限元中三節(jié)點(diǎn)三角形單元、四節(jié)點(diǎn)四邊形單元和九節(jié)點(diǎn)四邊形單元的形函數(shù)。圖1所示為局部坐標(biāo)系ξ-η平面上的標(biāo)準(zhǔn)拉格朗日單元。

對三角形單元可以直接利用廣義拉格朗日插值函數(shù)法(劃線法)構(gòu)造其插值函數(shù)，而構(gòu)造拉格朗日矩形單元插值函數(shù)的一個(gè)簡便而系統(tǒng)的方法是通過2個(gè)坐標(biāo)方向適當(dāng)方次的拉格朗日多項(xiàng)式的乘積得到[18]。圖1中3種單元對應(yīng)的拉格朗日多項(xiàng)式為[19]

L3:F1=1-ξ-η,F2=ξ,F3=η

(1)

k=1,2,3,4

(2)

L9:

(3)

式中：ξk和ηk為特征點(diǎn)k的坐標(biāo)。

圖1 拉格朗日單元Fig.1 Lagrange element

1.2 梁位移場

k=1,2,…,m

(4)

(5)

式中:uxk、uyk和uzk(k=1,2,3,4)為截面特征點(diǎn)3個(gè)方向上的位移分量。

圖2 截面插值梁示意圖Fig.2 Cross-section interpolation beam diagram

2 截面插值梁的有限元模型

2.1 梁單元位移函數(shù)

n節(jié)點(diǎn)梁單元截面特征點(diǎn)位移向量可表示為

(6)

(7)

(8)

2.2 應(yīng)變和應(yīng)力

將式(7)代入幾何方程和物理方程，可得單元應(yīng)變矩陣和應(yīng)力矩陣分別為

(9)

k=1,2,…,m;i=1,2,…,n

(10)

式中：L為微分算子矩陣，D彈性矩陣，考慮各向同性材料，其分別為

(11)

(12)

其中：

其中：E為彈性模量；υ為泊松比。

2.3 單元?jiǎng)偠染仃嚺c質(zhì)量矩陣

假定單元發(fā)生虛位移，根據(jù)式(7)和式(9)有

(13)

(14)

式中：l為截面特征點(diǎn)數(shù)；j為節(jié)點(diǎn)號(hào)。

若材料密度為ρ，根據(jù)達(dá)朗貝爾原理，單元內(nèi)單位體積分布慣性力為

(15)

在外力中考慮慣性力，則虛功原理可寫為

(16)

將式(10)、式(13)和式(14)代入整理可得式(16) 中各項(xiàng)分別為

(17)

(18)

(19)

式中：Pv和Ps分別為作用在單元上的體積力和面力；Pp為集中力。

式(16)進(jìn)一步整理可得

(20)

(21)

令

(22)

則式(21)可整理為

(23)

可得

(24)

式中：下標(biāo)組合lj代表著梁第j個(gè)節(jié)點(diǎn)所在平面的第l個(gè)插值點(diǎn)。

2.4 梁單元截面的離散與組集

對于復(fù)雜形狀截面梁單元，可以通過將橫截面細(xì)分為幾個(gè)拉格朗日單元的方法進(jìn)行有限元分析。由于位移函數(shù)中未知量為各插值點(diǎn)位移，所以單元疊加時(shí)只需疊加相應(yīng)插值點(diǎn)即可，如圖3所示。

圖3 單元截面疊加Fig.3 Superimposed cross-section element

3 算例1

考慮圖4所示矩形截面梁，梁長L=2 m，截面形狀為正方形，橫截面尺寸為a=0.1 m。彈性模量E=75 GPa，泊松比υ=0.33，材料密度為ρ=2 700 kg/m3。在梁自由端作用豎向集中力P=2 000 N，作用部位如圖4所示。

采用截面插值梁單元(L9)對梁結(jié)構(gòu)進(jìn)行分析，單元節(jié)點(diǎn)數(shù)取4，同時(shí)利用實(shí)體單元進(jìn)行建模作為對比，實(shí)體單元建模時(shí)在每一個(gè)特征點(diǎn)處設(shè)節(jié)點(diǎn)以保證2種模型有著相同的自由度。

圖5給出了隨梁單元數(shù)變化時(shí)2種模型自由端撓度uz的計(jì)算結(jié)果。從圖中可以看出，雖然實(shí)體單元模型和插值梁單元模型都對截面面內(nèi)進(jìn)行了離散，但是當(dāng)軸向網(wǎng)格數(shù)較少時(shí)，實(shí)體單元模型軸向尺寸將明顯大于其他2個(gè)方向尺寸，單元?jiǎng)偠染仃嚨难趴杀刃辛惺竭h(yuǎn)小于1，使得單元?jiǎng)偠染仃囍斜碚鞑煌较虻膭偠认禂?shù)相差過大，降低了撓度計(jì)算精度。而本文方法只要軸向劃分的單元數(shù)足夠描述梁的撓曲線特征即可達(dá)到滿意的計(jì)算精度。從圖5中可以看出，當(dāng)截面插值梁單元模型的單元數(shù)為5時(shí)，撓度計(jì)算結(jié)果即可收斂，而實(shí)體單元模型則需要遠(yuǎn)多于梁單元模型自由度時(shí)才能達(dá)到相同的收斂結(jié)果，因此截面插值梁單元在處理細(xì)長結(jié)構(gòu)時(shí)相比實(shí)體單元優(yōu)勢明顯。

圖4 懸臂梁示意圖Fig.4 Diagram of cantilever

圖5 撓度隨軸向單元數(shù)變化曲線Fig.5 Curves of variation of deflection with element number along axis

4 算例2

考慮圖6所示一端固支工字截面梁，其中W=100 mm，H=100 mm，t1=8 mm，t2=5 mm，長度L=1 000 mm。彈性模量E=210 GPa，泊松比υ=0.29。在自由端作用如圖6所示的一組集中力，大小為2 000 N。

根據(jù)Euler-Bernoulli梁以及Timoshenko梁理論，這些力自相平衡，截面不會(huì)發(fā)生變形，而在實(shí)際中很容易想象出，截面會(huì)發(fā)生明顯變形。截面插值梁模型可以完整描述這種變形，自由端截面變形如圖7所示。圖中7L9和14L9分別表示采用了7個(gè)和14個(gè)截面插值梁單元，從結(jié)果中也可以看出，通過增加截面面內(nèi)的單元數(shù)可以有效提高模型的精度。

圖6 工字梁截面示意圖Fig.6 I-shaped beam profile diagram

圖7 截面變形Fig.7 Cross-section deformation

5 算例3

考慮圖8所示直機(jī)翼結(jié)構(gòu)，機(jī)翼展長l=6 m，翼型為NACA2415，弦長c=1 m，蒙皮厚度為3 mm， 翼梁厚度為5 mm，其余具體尺寸如圖8所示。為方便分析，機(jī)翼各部件采用相同鋁合金材料，彈性模量E=75 GPa，泊松比υ=0.33，材料密度為ρ=2 700 kg/m3。

5.1 靜力學(xué)分析

采用截面插值梁單元對機(jī)翼各部件進(jìn)行建模分析，并分別與經(jīng)典梁單元以及Nastran實(shí)體單元計(jì)算結(jié)果進(jìn)行對比分析。截面插值點(diǎn)如圖9所示，截面通過圖示方法離散與組集，機(jī)翼展向單元數(shù)取8。經(jīng)典梁單元考慮Euler-Bernolli梁以及Timoshenko梁2種梁單元，機(jī)翼梁模型的各形狀屬性數(shù)據(jù)借助Patran軟件得到。同時(shí)為得到機(jī)翼結(jié)構(gòu)復(fù)雜應(yīng)力分布情況，對機(jī)翼結(jié)構(gòu)全部采用實(shí)體單元建模(如圖10所示)作為驗(yàn)證標(biāo)準(zhǔn)。實(shí)體模型中共有125 160個(gè)單元，其中124 140個(gè)Hex8單元，1 020個(gè)Wedge6單元，總節(jié)點(diǎn)數(shù)為169 290。

機(jī)翼坐標(biāo)系及本文計(jì)算輸出點(diǎn)如圖11所示。

圖8 機(jī)翼截面及尺寸Fig.8 Size and shape of wing cross-section

圖9 機(jī)翼截面插點(diǎn)值示意圖Fig.9 Interpolation points of wing cross-section diagram

圖10 Patran實(shí)體單元模型Fig.10 Solid element model of Patran

圖11 機(jī)翼坐標(biāo)系Fig.11 Coordinate system of wing

兼顧機(jī)翼實(shí)際受載情況(慣性載荷、局部質(zhì)量)以及簡化分析(集中載荷)的需要，針對表1所示工況條件展開分析，其中集中力及局部質(zhì)量點(diǎn)均位于點(diǎn)4、y=4 m處，g取9.8 m/s2，慣性力方向?yàn)閦方向。

本文采用73個(gè)L9單元模擬整片機(jī)翼。表2給出了2種工況條件下部分特征點(diǎn)的位移以及應(yīng)力值，并對截面插值梁模型、經(jīng)典梁模型計(jì)算結(jié)果與有限元軟件Nastran實(shí)體單元計(jì)算結(jié)果進(jìn)行了對比，其中EBBM表示Euler-Bernolli梁模型，TBM表示Timoshenko梁模型。表2中同時(shí)給出了各種模型自由度對比情況。

表1 工況表Table 1 Condition sheet

表2 部分位移及應(yīng)力結(jié)果Table 2 Partial displacement and stress results

圖12為機(jī)翼自由端截面變形圖。從表2及圖12可以看出，相比經(jīng)典梁單元，截面插值梁單元計(jì)算結(jié)果與實(shí)體單元計(jì)算結(jié)果更為接近，但差異很小。從自由度對比情況可以看出，截面插值梁模型自由度遠(yuǎn)小于實(shí)體模型自由度，因而有著更高的計(jì)算效率。

圖13為點(diǎn)1、點(diǎn)2、點(diǎn)3及點(diǎn)4所在軸線正應(yīng)力沿軸線y分布曲線圖，從圖中可以看出，截面插值梁單元計(jì)算結(jié)果精度整體高于經(jīng)典梁單元，尤其是圖13(d)中，在點(diǎn)4載荷作用部位正應(yīng)力變化趨勢突變，截面插值梁模型可以準(zhǔn)確反映出這種變化，與實(shí)體單元模型一致。

圖14為點(diǎn)5處剪應(yīng)力分布曲線，其中圖14(a)為沿軸線方向剪應(yīng)力分布圖，圖14(b)為截面y=3 m 處沿z方向剪應(yīng)力分布圖，從圖中可以看出，截面插值梁模型剪應(yīng)力計(jì)算結(jié)果與實(shí)體模型結(jié)果基本一致。圖15為點(diǎn)6處剪應(yīng)力分布圖，圖15(a)為沿軸線方向剪應(yīng)力分布圖，圖15(b)為截面y=3 m處沿z方向剪應(yīng)力分布。

綜合以上2種工況的靜力分析結(jié)果來看，雖然經(jīng)典梁單元可以得到部分點(diǎn)的位移分析結(jié)果，部分點(diǎn)的正應(yīng)力分析結(jié)果誤差也不大，但其無法得到結(jié)構(gòu)的剪應(yīng)力分布結(jié)果，而截面插值梁單元可以得到完整的靜力分析結(jié)果，且均有著堪比實(shí)體單元模型的計(jì)算精度，模型自由度也遠(yuǎn)小于實(shí)體模型，因而基于截面插值梁的機(jī)翼模型在靜力學(xué)分析中有著獨(dú)特的優(yōu)勢。

圖12 機(jī)翼自由端截面變形圖Fig.12 Deformation of wingtip corss-section

圖13 正應(yīng)力分布曲線(點(diǎn)1～4,工況1)Fig.13 Positive stress distribution curves (Points 1-4, Case 1)

圖14 剪應(yīng)力分布曲線(點(diǎn)5，工況1)Fig.14 Shear stress distribution curves (Point 5, Case 1)

5.2 動(dòng)力學(xué)分析

對機(jī)翼結(jié)構(gòu)進(jìn)行模態(tài)分析，分別考慮無集中質(zhì)量和含有集中質(zhì)量2種工況(如表1所示)。表3給出了幾種不同機(jī)翼模型的前8階固有頻率。從表中可以看出，經(jīng)典梁單元只能得到低階彎曲模態(tài)，無法得到較高階的扭轉(zhuǎn)模態(tài)、彎扭耦合模態(tài)。圖16給出了截面插值梁模型與實(shí)體模型振型結(jié)果對比圖(不含集中質(zhì)量)，圖17為相應(yīng)的MAC柱狀圖，對角線上前6階模態(tài)置信度均大于0.95，前8階大于0.90，說明采用本文方法的模型模態(tài)計(jì)算結(jié)果與精細(xì)化的實(shí)體單元機(jī)翼模態(tài)計(jì)算結(jié)果吻合程度很好。

圖15 剪應(yīng)力分布曲線(點(diǎn)6，工況1)Fig.15 Shear stress distribution curves (Point 6, Case 1)

表3 固有頻率計(jì)算結(jié)果Table 3 Natural frequency calculation results

圖16 振型結(jié)果對比Fig.16 Comparison of mode shapes

在模態(tài)分析的基礎(chǔ)上，對該機(jī)翼展開顫振分析。為簡化分析，在空氣動(dòng)力方面采用片條理論，直接利用二維流的已知結(jié)果計(jì)算氣動(dòng)力。在顫振初步估算階段，對于大展弦比飛機(jī)，無集中質(zhì)量的情況下可取機(jī)翼的一階彎曲和扭轉(zhuǎn)振型作為位移函數(shù)將機(jī)翼運(yùn)動(dòng)簡化為2個(gè)自由度[20-21]。利用v-g法進(jìn)行顫振分析，相應(yīng)的v-g曲線和v-f曲線如圖18所示。

圖17 MAC柱狀圖Fig.17 MAC histogram

圖18 顫振計(jì)算結(jié)果Fig.18 Flutter calculation results

通過插值的方法得到g=0時(shí)對應(yīng)的顫振速度和相應(yīng)的顫振頻率，表4中給出了2種模型的顫振特性對比結(jié)果。從計(jì)算結(jié)果可以看出，截面插值梁模型在顫振初步估算階段可以得到較為滿意的顫振計(jì)算結(jié)果。

為了進(jìn)一步驗(yàn)證本文模型的有效性，對截面插值梁機(jī)翼模型展開時(shí)域響應(yīng)分析。在點(diǎn)4、y=4 m 處沿z方向施加幅值為3 000 N、頻率為20 Hz的正弦激勵(lì)。分別采用振型疊加法和Newmark 法進(jìn)行動(dòng)響應(yīng)計(jì)算，并與Nastran實(shí)體單元中的模態(tài)疊加和直接積分2種計(jì)算結(jié)果進(jìn)行對比。部分點(diǎn)位移響應(yīng)如圖19和圖20所示。

圖19為振型疊加法計(jì)算結(jié)果對比，選取的模態(tài)數(shù)為10，從圖中可以看出，2種模型計(jì)算結(jié)果基本一致，這也進(jìn)一步驗(yàn)證了前述模態(tài)分析結(jié)果的可靠性。圖20為直接積分法計(jì)算結(jié)果，從圖中也可以看出，2種模型計(jì)算結(jié)果一致，基于截面插值梁的機(jī)翼模型可以很好地用于機(jī)翼結(jié)構(gòu)的動(dòng)響應(yīng)分析。

表4 顫振計(jì)算結(jié)果Table 4 Flutter calculation results

圖19 位移響應(yīng)圖(振型疊加法)Fig.19 Displacement results (modal superposition method)

圖20 位移響應(yīng)圖(直接積分法)Fig.20 Displacement results (direct integration method)

6 結(jié) 論

1) 截面插值梁單元通過截面特征點(diǎn)位移插值得到截面位移場，可以得到截面面內(nèi)、面外完整截面變形。

2) 截面插值梁單元與經(jīng)典梁單元相比，不僅有著更高的計(jì)算精度，還可以有效處理復(fù)雜截面形狀細(xì)長結(jié)構(gòu)；與實(shí)體單元相比，同等網(wǎng)格尺寸、同樣誤差等級下，截面插值梁單元模型自由度遠(yuǎn)遠(yuǎn)小于實(shí)體模型。

3) 截面插值梁單元可以用于簡單機(jī)翼結(jié)構(gòu)建模，并能夠得到滿意的靜力學(xué)、動(dòng)力學(xué)分析結(jié)果，為機(jī)翼的建模與分析提供了一種一維梁簡化建模方法。