一種考慮截面任意變形模式梁的有限元模型

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

1.南京航空航天大學 機械結構力學及控制國家重點實驗室,南京 210016 2. 南京航空航天大學 振動工程研究所,南京 210016

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

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

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

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

1 構造梁位移場

1.1 拉格朗日插值多項式介紹

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

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

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

(1)

k=1,2,3,4

(2)

L9:

(3)

式中：ξk和ηk為特征點k的坐標。

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

1.2 梁位移場

k=1,2,…,m

(4)

(5)

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

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

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

2.1 梁單元位移函數

n節點梁單元截面特征點位移向量可表示為

(6)

(7)

(8)

2.2 應變和應力

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

(9)

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

(10)

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

(11)

(12)

其中：

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

2.3 單元剛度矩陣與質量矩陣

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

(13)

(14)

式中：l為截面特征點數；j為節點號。

若材料密度為ρ，根據達朗貝爾原理，單元內單位體積分布慣性力為

(15)

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

(16)

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

(17)

(18)

(19)

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

式(16)進一步整理可得

(20)

(21)

令

(22)

則式(21)可整理為

(23)

可得

(24)

式中：下標組合lj代表著梁第j個節點所在平面的第l個插值點。

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

對于復雜形狀截面梁單元，可以通過將橫截面細分為幾個拉格朗日單元的方法進行有限元分析。由于位移函數中未知量為各插值點位移，所以單元疊加時只需疊加相應插值點即可，如圖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)對梁結構進行分析，單元節點數取4，同時利用實體單元進行建模作為對比，實體單元建模時在每一個特征點處設節點以保證2種模型有著相同的自由度。

圖5給出了隨梁單元數變化時2種模型自由端撓度uz的計算結果。從圖中可以看出，雖然實體單元模型和插值梁單元模型都對截面面內進行了離散，但是當軸向網格數較少時，實體單元模型軸向尺寸將明顯大于其他2個方向尺寸，單元剛度矩陣的雅可比行列式遠小于1，使得單元剛度矩陣中表征不同方向的剛度系數相差過大，降低了撓度計算精度。而本文方法只要軸向劃分的單元數足夠描述梁的撓曲線特征即可達到滿意的計算精度。從圖5中可以看出，當截面插值梁單元模型的單元數為5時，撓度計算結果即可收斂，而實體單元模型則需要遠多于梁單元模型自由度時才能達到相同的收斂結果，因此截面插值梁單元在處理細長結構時相比實體單元優勢明顯。

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

圖5 撓度隨軸向單元數變化曲線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。

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

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

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

5 算例3

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

5.1 靜力學分析

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

機翼坐標系及本文計算輸出點如圖11所示。

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

圖9 機翼截面插點值示意圖Fig.9 Interpolation points of wing cross-section diagram

圖10 Patran實體單元模型Fig.10 Solid element model of Patran

圖11 機翼坐標系Fig.11 Coordinate system of wing

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

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

表1 工況表Table 1 Condition sheet

表2 部分位移及應力結果Table 2 Partial displacement and stress results

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

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

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

綜合以上2種工況的靜力分析結果來看，雖然經典梁單元可以得到部分點的位移分析結果，部分點的正應力分析結果誤差也不大，但其無法得到結構的剪應力分布結果，而截面插值梁單元可以得到完整的靜力分析結果，且均有著堪比實體單元模型的計算精度，模型自由度也遠小于實體模型，因而基于截面插值梁的機翼模型在靜力學分析中有著獨特的優勢。

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

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

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

5.2 動力學分析

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

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

表3 固有頻率計算結果Table 3 Natural frequency calculation results

圖16 振型結果對比Fig.16 Comparison of mode shapes

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

圖17 MAC柱狀圖Fig.17 MAC histogram

圖18 顫振計算結果Fig.18 Flutter calculation results

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

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

圖19為振型疊加法計算結果對比，選取的模態數為10，從圖中可以看出，2種模型計算結果基本一致，這也進一步驗證了前述模態分析結果的可靠性。圖20為直接積分法計算結果，從圖中也可以看出，2種模型計算結果一致，基于截面插值梁的機翼模型可以很好地用于機翼結構的動響應分析。

表4 顫振計算結果Table 4 Flutter calculation results

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

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

6 結 論

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

2) 截面插值梁單元與經典梁單元相比，不僅有著更高的計算精度，還可以有效處理復雜截面形狀細長結構；與實體單元相比，同等網格尺寸、同樣誤差等級下，截面插值梁單元模型自由度遠遠小于實體模型。

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