一類考慮溫度效應的截面插值梁模型

潘成浩， 陳國平,2， 何 歡,2

(1. 南京航空航天大學 機械結(jié)構(gòu)力學及控制國家重點實驗室， 南京 210016；2. 南京航空航天大學 振動工程研究所， 南京 210016)

承受橫向載荷而受彎的構(gòu)件稱為梁。它是生產(chǎn)生活中常見的結(jié)構(gòu)之一，也是工程實際中廣泛采用的基本構(gòu)件之一。常見的梁結(jié)構(gòu)可以分為實體梁和薄壁梁兩大類：實體梁變形以彎曲為主，薄壁梁變形以彎曲和扭轉(zhuǎn)為主。傳統(tǒng)的梁理論本質(zhì)上是將三維的實際梁結(jié)構(gòu)基于一定的假設(shè)，通過對梁的軸線和截面分別進行一定的處理，最終轉(zhuǎn)化為一維問題，以此簡化計算。由此建立的梁理論有經(jīng)典梁理論(Euler-Bernoulli梁理論)、一階剪切變形理論(Timoshenko梁理論)[1]以及高階剪切變形理論(Reddy三階梁理論)等。這些梁理論只適用于一定的范圍，有一定的局限性。如Euler-Bernoulli梁理論適用于受剪切變形影響較小的細長梁的求解問題，但其忽略了橫向剪應力和橫向剪應變的影響，不適用于短梁、復合材料梁、夾層結(jié)構(gòu)以及梁的振動等問題。Timoshenko梁理論[2]考慮了梁的彎曲變形引起的轉(zhuǎn)動慣量和梁的剪切變形，適用于短深梁，但仍存在能量不自洽、剪切自鎖等問題。許多學者對于上述三種理論進行了深入研究：Newmark[3]、Levinson[4]、Li[5]、周倩南[6]、何歡等[7]、Wan等[8]針對這些局限性，考慮剪切滑移、橫截面翹曲，提出了多種新的梁模型。經(jīng)典梁理論由于采用了平截面假設(shè)，所以在截面變形的描述上有很大的不足，例如機翼在復雜載荷作用下翹曲變形較大，傳統(tǒng)的梁理論不能反映復雜的截面變形。此外，針對復雜截面而言，傳統(tǒng)的梁理論不能精確地描述截面的幾何形狀，由此對力學特性的描述就不夠精確。商業(yè)有限元軟件在處理復雜截面、復雜工況的梁結(jié)構(gòu)時推薦采用實體單元，但是這又帶來自由度增加計算規(guī)模增大的問題。

為此，本文提出一種基于截面插值的考慮溫度效應和截面形狀的梁的有限元單元模型。該模型摒棄之前梁理論的剛性軸假設(shè)和中性軸假設(shè)，將梁單元分為在截面上的截面特征點和沿軸向的單元節(jié)點，如此就可以由軸向單元節(jié)點的位移經(jīng)插值得到任意一個截面上截面特征點的位移，再由截面特征點的位移插值得到該截面內(nèi)任意一點的位移，這樣就得到了截面插值梁模型的位移場。經(jīng)過這種方法得到的位移場能夠很輕易地描述截面的形狀，同時在面對溫度梯度時，也能很好地反映不同溫度帶來的不同變形。再基于虛功原理得到單元剛度矩陣、溫度剛度矩陣和質(zhì)量矩陣，對結(jié)構(gòu)進行靜力分析和模態(tài)分析。

另外，熱模態(tài)分析是研究結(jié)構(gòu)的熱振動和熱顫振的基礎(chǔ)。分析熱模態(tài)關(guān)鍵是分析溫度變化對結(jié)構(gòu)的影響。而溫度變化對結(jié)構(gòu)的影響主要體現(xiàn)在兩方面，一是溫度變化影響材料的特性，例如在高溫情況下某些金屬材料的彈性模量會降低，影響結(jié)構(gòu)的整體剛度。二是溫度變化產(chǎn)生的熱應力會形成初應力剛度，影響結(jié)構(gòu)的剛度。國內(nèi)外學者在這方面進行了大量的研究。吳曉等[9-12]研究了材料特性隨溫度變化后對結(jié)構(gòu)振動的影響；Snyder等[13-15]則研究了熱應力對結(jié)構(gòu)振動的影響。Kehoe等[16-19]研究對比了這兩方面各自對結(jié)構(gòu)振動的影響。

針對提出的考慮溫度效應的截面插值梁模型，本文同樣研究了其熱模態(tài)的計算，并與Nastran實體單元計算結(jié)果比較，驗證模型的可靠性。

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

本文中梁位移場的構(gòu)造需要進行兩步插值：第一步在截面特征點之間進行插值，第二步要在軸向節(jié)點之間進行插值。本文采用的插值方法均為拉格朗日多項式插值。

對于第一步插值采用二元拉格朗日多項式進行插值。根據(jù)在自身處二元拉格朗日插值多項式取值為1，在其他點處多項式取值0，并且所有插值多項式的和恒為1這個特點，可以很輕易地構(gòu)造出插值函數(shù)。本文的插值函數(shù)也對應于傳統(tǒng)單元的插值函數(shù)：三節(jié)點插值多項式對應于三節(jié)點三角形單元的形函數(shù)，四節(jié)點則對應于四節(jié)點矩形單元，九節(jié)點則對應于九節(jié)點矩形單元。圖1為局部坐標系下的標準拉格朗日單元。

(a) L3

(b) L4

(c) L9圖1 局部坐標系下的標準拉格朗日單元Fig.1 Standard Lagrange element in local coordinate system

由此可以得到對應的拉格朗日插值多項式

F1=1-ξ-η

L3:F2=ξ

(1)

F3=η

(2)

(3)

式中，ξk，ηk為截面特征點k的坐標。

對于第二步的插值直接采用一維拉格朗日插值多項式，對于n個節(jié)點的一維單元，Ni(y)可采用n-1次拉格朗日插值多項式，即

(4)

式中：n是單元軸向的節(jié)點數(shù)；y1,y2,…,yn是n個節(jié)點的坐標。

利用拉格朗日插值函數(shù)可以得到梁單元的位移場

u=FkNiuki,k=1,2,…,m,i=1,2,…,n

(5)

式中：u為單元內(nèi)任意一點的位移向量；m為截面特征點數(shù)；n為軸向節(jié)點數(shù)；uki為軸向節(jié)點i所在截面內(nèi)第k個截面特征點的位移向量

uki=[ukix,ukiy,ukiz]T

(6)

用矩陣來表示就是

u=FN·u′

(7)

其中：

(8)

(9)

在之后的計算中還需要根據(jù)已知的節(jié)點溫度來插值得到單元中任意一點的溫度，采用的插值方法和位移場的構(gòu)建方法相同，故可以得到單元內(nèi)任意一點的溫度為

T=FkNiTki

(10)

式中，Tki為軸向節(jié)點i所在截面內(nèi)第k個截面特征點的溫度。同樣可以用矩陣表示

T=FN′·T′

(11)

2 有限元分析

2.1 應變和應力

將式(7)代入幾何方程，可以導出單元總應變與單元節(jié)點位移的關(guān)系，即

ε=Lu=L·FN·u′

(12)

式中，L為微分算子矩陣

(13)

溫度變化在單元內(nèi)產(chǎn)生的應變，即初應變?yōu)?/p>

(14)

式中，α為線膨脹系數(shù)。

將式(12)和式(14)代入熱彈性理論下的物理方程，可以導出單元應力與單元節(jié)點位移、單元節(jié)點溫度的關(guān)系，即

σ=D-1(ε-εT)=D-1(L·FN·u′-αFN′·

T′·a)

(15)

式中：ε為總應變，是彈性體在外載荷作用下產(chǎn)生的彈性應變和初應變εT兩者之和；D為彈性矩陣

(16)

式中：ν為泊松比；E為彈性模量；G為剪切模量。

2.2 單元剛度矩陣和質(zhì)量矩陣

假設(shè)單元發(fā)生虛位移，則虛位移和相應的虛應變?yōu)?/p>

δu=FN·δu′

(17)

δε=Lδu=L·FN·δu′

(18)

故由虛功原理可知，單元中內(nèi)力所做的虛功與外力所做的虛功總和為零

δWI+δWE=0

(19)

單元中內(nèi)力所做的虛功為

(20)

單元中外力所做的虛功為

δWE=δWext+δWine

(21)

式中：δWext為外力(除去慣性力)所做虛功；δWine為慣性力所做虛功

(22)

(23)

式中：PV為作用在單元上的體積力；PS為作用在單元上的面力；PP則為集中力；ρ為單元體的密度。

令

(24)

整理式(17)至式(24)，可得

(25)

化簡可得

(26)

通過上式可計算出熱應力。

2.3 單元組集

對于截面插值梁單元的單元組集過程和傳統(tǒng)有限元計算的過程大致相同，都是將矩陣中的各元素按照結(jié)構(gòu)節(jié)點自由度排列順序“對號入座”疊加而成。但與傳統(tǒng)梁單元不同的是，截面插值梁單元不僅能在梁的軸向方向上組集，還能在截面內(nèi)通過共用特征點的方式組集。這能很準確地描述復雜的截面形狀，同時也能通過局部增加截面插值梁單元的密度來提高有限元計算的精度。

2.4 求解熱模態(tài)

利用計算得到的熱應力推導熱模態(tài)的計算公式，這需要更新剛度矩陣，把熱應力的影響考慮進去?？紤]在實際情況下，梁模型在受熱后主要受力為軸向的熱應力σaxial，故類比于受軸力作用的梁，受熱梁振動時的內(nèi)力虛功包括微小振動帶來的虛功和考慮由軸向熱應力引起的橫向剪切力在對應剪應變上做的虛功。而由軸向熱應力引起的橫向剪切力為：σaxial·usec1,y和σaxial·usec3,y對應的剪應變?yōu)棣膗sec1,y和δusec3,y。其中下標“sec1”，“sec3”表示截面內(nèi)的兩個坐標軸方向，下標“,y”表示對y求偏導。則內(nèi)力虛功變?yōu)椋?/p>

FNsec3,ydV)u

(27)

則更新后的剛度矩陣為

(28)

最后利用更新后的剛度矩陣和質(zhì)量矩陣就可計算得到熱模態(tài)。

3 算例1

考察如圖2所示的矩形截面梁，其尺寸為：a=20 mm，b=100 mm，長度L=2 000 mm。彈性模量E=200 GPa，泊松比ν=0.25，材料密度ρ=2 700 kg/m3，熱膨脹系數(shù)α=11×10-6。該梁的約束情況為一端固支一端自由。

圖2 矩形截面梁Fig.2 Diagram of rectangular section beam

對矩形截面梁進行單元劃分，在截面上劃分兩個單元，每個單元截面采用九點插值(L9),軸向取四個節(jié)點。同時在矩形截面梁上作用一個溫度梯度T=180(20z+1)℃。利用上面的公式進行單元分析，將得到的結(jié)果與Nastran軟件中的Solid單元計算結(jié)果進行比對分析，注意Solid單元在每個特征點處均有節(jié)點，以此保證截面插值梁單元與Solid單元有相同的自由度。

圖3給出了隨著軸向單元數(shù)的增加，截面插值梁單元2L9與Solid單元自由端撓度的變化。從圖中可以看出，截面插值梁單元的撓度計算結(jié)果比Solid單元的撓度計算結(jié)果更快收斂。由此可以看出截面插值梁模型在計算細長梁時比Solid單元更快捷，更有優(yōu)勢。

圖3 自由端撓度隨軸向單元數(shù)變化曲線Fig.3 Deflection varies with axial element number

4 算例2

為驗證截面插值梁模型對復雜截面的處理的準確性，以圖4所示的工字梁為例進行計算。其尺寸為：W=100 mm，H=100 mm，t1=8 mm，t2=5 mm，長度L=1 000 mm。材料的彈性模量E=210 GPa，泊松比ν=0.29，密度ρ=2 700 kg/m3，熱膨脹系數(shù)α=11×10-6。

圖4 工字梁截面Fig.4 I-shaped profile diagram

對工字型截面采取如圖4的單元劃分形式，共劃分14個單元，每個單元的截面都采用九點插值(L9)，同時在軸向上劃分20個單元，每個單元軸向取四個節(jié)點。該梁的約束情況為兩端固支。

圖5(a)是在工字梁上作用均勻溫度場T=10 ℃的截面變形圖(變形放大1 000倍)。圖5(b)是在T=100 ℃時的截面變形圖(變形放大100倍)。圖5(c)是在T=1 200z+60 ℃時的截面變形圖(變形放大100倍)。圖6(a)和(c)分別是均勻溫度場T=100 ℃和非均勻溫度場T=1 200z+60 ℃下點1(0.05,0,0.05)所在軸線的x向撓度，圖6(b)和(d)分別是均勻溫度場和非均勻溫度場下點2(0,0,0)所在軸線的y向位移。均勻溫度場下比較了截面插值梁模型(14L9)、Nastran的Solid單元和Beam單元的計算結(jié)果；而由經(jīng)典梁理論導出計算的Beam單元無法施加非均勻溫度場T=1 200z+60 ℃，所以非均勻溫度場下只比較了截面插值梁模型(14L9)和Nastran的Solid單元的計算結(jié)果。從圖中可以看出無論是作用均勻溫度場還是非均勻溫度場，截面插值梁模型的計算結(jié)果與Solid單元計算結(jié)果一致，精度很高，與Beam單元相比，截面插值梁模型能更好地反映梁截面上的變形。

(a) T=10 ℃

(b) T=100 ℃

(c) T=1 200z+60 ℃圖5 截面變形圖Fig.5 Cross-section deformation

(a) 點1，T=100 ℃

(b) 點2，T=100 ℃

(c) 點1，T=1 200z+60 ℃

(d) 點2，T=1 200z+60 ℃圖6 各向位移曲線圖Fig.6 The displacement curve

5 算例3

同樣考慮算例1的矩形截面梁，梁長度改為L=1 000 mm，在截面上還是劃分兩個單元，在軸向上劃分20個單元。約束情況改為兩端固支。另外用商業(yè)有限元軟件Nastran里的Solid單元和Beam單元進行建模，建模中軸向取相同的節(jié)點數(shù)。

5.1 靜力學分析

靜力學分析中考察的情況是在矩形截面梁上作用溫度載荷，分為兩種工況，第一種是均勻溫度場T=100 ℃，第二種是非均勻溫度場T=1 800z+90 ℃。圖7是點1(0.01,0.05,0.467)(即截面右上角頂點)所在截面的變形圖對比(變形放大100倍)。圖8是點1和點2(0.01,0,-0.05)(即截面右下角頂點)所在軸線的軸向應力對比。表1是兩種溫度場下不同建模方法計算的點1的位移和應力值，從圖中可以看出截面插值梁模型不僅準確地描述了截面變形，也能準確地計算軸向應力，尤其是在兩端約束處的應力激增位置，截面插值梁模型仍可以準確反映這種變化。

(a) T=100 ℃

(b) T=1 800z+90 ℃圖7 截面變形圖Fig.7 Cross-section deformation

(a) T=100 ℃

(b) T=1 800z+90 ℃圖8 不同溫度下正應力分布曲線Fig.8 Positive stress distribution curves under different temperatures

表1 點1位移和應力結(jié)果Tab.1 Displacement and stress result of Point 1

5.2 動力學分析

對矩形截面梁進行模態(tài)分析，表2是不考慮溫度的模態(tài)對比。表3是考慮均勻溫度場的模態(tài)對比。表4是考慮非均勻溫度場的模態(tài)對比。從表中可以看出，Beam單元無法得到扭轉(zhuǎn)模態(tài)。而截面插值梁單元則不會缺失扭轉(zhuǎn)模態(tài)，而且其余模態(tài)的精度相對更高。圖9給出了Solid單元和截面插值梁單元的模態(tài)振型對比(考慮非均勻溫度場)，通過振型比較可以發(fā)現(xiàn)，兩種方法的模態(tài)振型是一致的。圖10是對應的MAC置信度柱狀圖，對角線上前10階的模態(tài)置信度均大于0.99，說明截面插值梁單元的模態(tài)計算結(jié)果與Solid單元的計算結(jié)果吻合程度很好。

(a) 第1階

(b) 第2階

(c) 第3階

(d) 第4階

(e) 第5階

(f) 第6階圖9 非均勻溫度場下模態(tài)振型對比Fig.9 Comparison of mode shapes under non-uniform temperature

圖10 非均勻溫度下MAC柱狀圖Fig.10 MAC histogram under non-uniform temperature

表2 不考慮溫度的模態(tài)結(jié)果Tab.2 Modal results regardless of temperature

表3 均勻溫度場下模態(tài)結(jié)果Tab.3 Modal results under uniform temperature

表4 非均勻溫度場下模態(tài)結(jié)果Tab.4 Modal results under non-uniform temperature

從靜力學分析與動力學分析結(jié)果來看，截面插值梁模型可以很好地應用于復雜截面梁在復雜溫度情況下的響應分析。

6 結(jié) 論

本文提出了一種基于截面插值的空間梁模型，并研究了在考慮溫度效應的情況下模型的運用。與傳統(tǒng)梁模型相比，截面插值梁模型能很好地反映截面內(nèi)的各種情況，比如截面的變形、截面內(nèi)的溫度梯度、截面的復雜形狀，有著更高的計算精度。另外與Solid單元相比，在處理細長結(jié)構(gòu)上，截面插值梁模型需要更少的節(jié)點，更少的自由度，同時能保證計算結(jié)果的精度與可靠性。