狀態(tài)空間有限元法分析固支疊層結(jié)構(gòu)熱應(yīng)力

沈宇昂，梁擁成，鄭興偉，韓志林

(東華大學(xué) 理學(xué)院, 上海 201620)

近年來,強(qiáng)度高且性能好的復(fù)合疊層材料被廣泛用于工程界,但疊層板之間為協(xié)調(diào)變形,各層之間必然存在相互制約的作用力,這便是層間應(yīng)力,準(zhǔn)確計(jì)算層間應(yīng)力對于工程界而言意義重大。雖然有限差分法、有限元法、邊界元法等都是可以計(jì)算層間應(yīng)力的數(shù)值方法,但是使用狀態(tài)空間法分析計(jì)算疊層結(jié)構(gòu)的力學(xué)特征,具有計(jì)算精度較高、計(jì)算量較小和計(jì)算所需的時(shí)間較短等優(yōu)點(diǎn)。輸入變量、狀態(tài)轉(zhuǎn)移矩陣和輸出變量是狀態(tài)空間法的3部分。疊層結(jié)構(gòu)交界面上已知的位移和應(yīng)力分量組成了輸入變量,材料的屬性決定了狀態(tài)轉(zhuǎn)移矩陣,而各個(gè)交界面上未知的位移、應(yīng)力分量組成了輸出變量。狀態(tài)空間法已被用于許多力學(xué)問題:范家讓[1]提出將狀態(tài)空間法運(yùn)用于板殼問題;謝霽明[2]和郭增偉等[3]將狀態(tài)空間法用于計(jì)算橋梁顫振分析;歐陽光等[4]將狀態(tài)空間法用于研究吊桿張拉力分析。狀態(tài)空間法不僅可以計(jì)算力學(xué)物理量,還可以用于檢驗(yàn)材料性質(zhì)或材料建模等,如:Wu等[5]使用狀態(tài)空間法來檢驗(yàn)彈性層表面的不穩(wěn)定性;王延靈等[6]將狀態(tài)空間法用于大迎角氣動(dòng)力學(xué)建模。

受有限元思想的啟發(fā),Sheng等[7]提出狀態(tài)空間有限元法,此方法在疊層結(jié)構(gòu)的水平方向采用有限元法離散位移、應(yīng)力分量,而在垂直方向仍然采用狀態(tài)空間法求解;Wu等[8]提出沿疊層結(jié)構(gòu)的水平方向采用樣條插值法,而在垂直方向使用狀態(tài)空間法,這也可以較為準(zhǔn)確地計(jì)算出疊層結(jié)構(gòu)的物理量。

現(xiàn)有研究已將狀態(tài)空間法用于分析熱應(yīng)力問題。例如:鄒貴平等[9]在分析層合圓柱厚殼熱應(yīng)力問題時(shí)將狀態(tài)空間法與打靶法結(jié)合;盛宏玉等[10]引入狀態(tài)空間理論并應(yīng)用差分格式分析了一維桿的縱向熱彈性動(dòng)力響應(yīng)瞬態(tài)問題;韓志林等[11]將狀態(tài)空間法用于分析疊層結(jié)構(gòu)熱應(yīng)力問題;陳新宇等[12]使用狀態(tài)空間法研究受軸對稱溫度載荷的功能梯度圓筒問題,得到正確、有效的數(shù)值解。此外,狀態(tài)空間法可以與有限元法結(jié)合,其他計(jì)算方法也可以與狀態(tài)空間法耦合應(yīng)用,例如:Benedetti等[13]將狀態(tài)空間法與快速邊界元法結(jié)合,求得壓電黏結(jié)材料的應(yīng)力場;Cheng等[14]將狀態(tài)空間法與邊界元法結(jié)合,求解功能梯度材料的彈性力學(xué)問題。

通常研究人員會(huì)選擇在水平面使用Lagrange單元來離散物理場,例如Wu等[15]使用高階Lagrange插值來計(jì)算鹽度,但由于本文研究的疊層結(jié)構(gòu)兩固支端的應(yīng)力變化較為劇烈,使用Lagrange單元離散物理場不一定能準(zhǔn)確描述應(yīng)力變化。因此本文采用更光滑的B樣條單元離散物理場,從而提出等幾何的B樣條狀態(tài)空間法。B樣條插值也廣泛用于力學(xué)分析,如:Simpson等[16]結(jié)合B樣條插值與邊界元法,進(jìn)行彈性力學(xué)分析;范紀(jì)華等[17]使用B樣條插值對旋轉(zhuǎn)懸臂梁進(jìn)行動(dòng)力學(xué)分析。

同時(shí),本文分別采用通解法和精細(xì)積分法求解狀態(tài)空間方程,以研究何種方法更適用于固支問題。本文數(shù)值算例表明,位移的數(shù)值結(jié)果在2種離散方案及2種求解方案情況下均較為精準(zhǔn),而在計(jì)算應(yīng)力時(shí),B樣條結(jié)果精度更高。雖然精細(xì)積分法和通解法均能較快完成二維狀態(tài)空間方程的計(jì)算,其中精細(xì)積分法由于采用高效迭代的技術(shù),其計(jì)算時(shí)長少于通解法。

1 二維模型變分原理

Lagrange單元離散和B樣條單元離散疊層結(jié)構(gòu)如圖1所示,由M(M=3)層各向同性材料所組成的疊層結(jié)構(gòu),不同層的彈性模量E、泊松比μ以及熱膨脹系數(shù)α均不同。橫向位移u和縱向位移v分別為x和y方向的位移;l和h分別為整體材料的長度和高度;ΔQ為溫度荷載,疊層結(jié)構(gòu)兩側(cè)均受固支約束。由于材料的幾何形狀及溫度荷載均對稱,故只截取左半邊進(jìn)行分析。邊界條件:x=0的邊的位移u=v=0;x=l/2的邊的位移u=0,應(yīng)力σxy=0。

圖1 Lagrange單元和B樣條單元離散疊層結(jié)構(gòu)Fig.1 Discretization of the laminated structure by Lagrange and B-spline elements

根據(jù)Hellinger-Reissner[18]變分原理,以位移和應(yīng)力滿足邊界條件的假設(shè)對疊層結(jié)構(gòu)中的任意一層建立二維彈性力學(xué)中的平衡方程式(1)和幾何方程式(2)。

(1)

(2)

本文僅在疊層結(jié)構(gòu)第一層最上方劃分單元,分別通過Lagrange插值和B樣條插值來劃分單元。以Lagrange二次元為例,在第一層材料的上方劃分n個(gè)單元,產(chǎn)生了2n+1個(gè)節(jié)點(diǎn)(見圖1(a))。

Lagrange插值公式為

(3)

式中:w為lagrange公式局部坐標(biāo)系的坐標(biāo),wi為第i個(gè)lagrange節(jié)點(diǎn)的坐標(biāo)zi為對應(yīng)的函數(shù)值;pk(w)為第k個(gè)基函數(shù),在wk處值為1,其余插值點(diǎn)處值為0;Bk為Lagrange插值法產(chǎn)生的所有節(jié)點(diǎn)的集合;Ln(w)為總體插值函數(shù),其經(jīng)過每一個(gè)插值點(diǎn)。

B樣條曲線中的遞推形B樣條基函數(shù)[19]多被用于CAD(computer aided design)。B樣條與數(shù)值算法的結(jié)合可以直接套用CAD模型,從而大幅減少劃分單元的時(shí)間。B樣條插值法與Lagrange插值法有所不同,最顯著的特點(diǎn)在于B樣條不滿足Kronecker-Delta特性。B樣條的基函數(shù)遞推式如式(4)和(5)所示,其一階導(dǎo)函數(shù)如式(6)所示。

(4)

(5)

(6)

式中:wi+1和wi+p分別表示第i+1和第i+p個(gè)控制節(jié)點(diǎn)的坐標(biāo);Ni,0為0階B樣條的第i個(gè)基函數(shù);Ni,p表示p階B樣條的第i個(gè)基函數(shù)。由式(4)和(5)可知,p階B樣條基函數(shù)可由0階基函數(shù)迭代得到。

B樣條插值法的一大特點(diǎn)在于B樣條插值需要節(jié)點(diǎn)向量和控制點(diǎn)共同作用,節(jié)點(diǎn)向量確定B樣條的基函數(shù),而控制點(diǎn)確定B樣條曲線的形狀。本文所求物理量反映控制點(diǎn)所具有的各分量,而真實(shí)物理量需根據(jù)B樣條形函數(shù)進(jìn)一步插值得到。k1,k2,…,kn+1為B樣條的節(jié)點(diǎn),每兩個(gè)相鄰節(jié)點(diǎn)之間為一個(gè)單元;c1,c2,…,cn-2為控制點(diǎn),對二次插值來說,相鄰的3個(gè)控制點(diǎn)作用在一個(gè)單元內(nèi)(見圖1(b))。

兩種單元對應(yīng)的離散方法類似,以位移u及二次元為例,u=N1(x1)u1+N2(x2)u2+N3(x3)u3。

對于本文所探討的熱應(yīng)力問題,其物理場方程為

ε=Sσ+J

(7)

式中:S為各向同性材料的柔度矩陣;J=[αQαQ0]T,Q為每層材料的溫度,℃。

將式(7)代入式(2),并將σxx和σxy、σyy拆分,可得到

(8)

αK)dΩ=0

(9)

式中:S1、S2、S3分別為柔度矩陣S中的元素。將離散每一層材料的界面物理量分別代入式(9)得到式(10),代入式(8)并與式(1)進(jìn)行合并得到式(11)。

(10)

(11)

引入邊界條件至式(10)和(11),并由變分理論等價(jià)得式(12)。

Dfsf=EfRf-GfK

(12)

(13)

聯(lián)立方程式(12)和(13),消去sf得到式(14)。

(14)

計(jì)算由式(13)得到微分方程時(shí),需要保障最后的系數(shù)矩陣是正定的,故要求未知狀態(tài)量中位移和應(yīng)力分量的個(gè)數(shù)相同。然而通過固支邊界條件可知,每層固支邊都存在多余未知量σxy和σyy。為保證正定性,利用式(15)消除這兩個(gè)分量。

(15)

式中:K1是節(jié)點(diǎn)1處的溫度。

將式(15)計(jì)算所得的關(guān)系式代入式(14)中,形成任意一層的非齊次狀態(tài)方程如式(16)所示。

(16)

2 狀態(tài)方程的求解

2.1 通解法

對于疊層結(jié)構(gòu)中的第i層,得到如式(16)的非齊次狀態(tài)方程為

(17)

式中:Ri為狀態(tài)方程中的狀態(tài)向量,是僅關(guān)于y的函數(shù);Ni為狀態(tài)方程中的狀態(tài)矩陣;Mi為狀態(tài)方程中的輸入矩陣。

根據(jù)狀態(tài)方程的通解,并令y=yi,可得:

(18)

式中:Di(y)為狀態(tài)方程中的狀態(tài)轉(zhuǎn)移矩陣;hi為第i層材料的厚度。

Aixi=bi

(19)

式中:Ai為第i層材料的系數(shù)矩陣;bi為已知狀態(tài)分量組成的向量;xi為未知狀態(tài)分量組成的向量。

同樣的,第i+1層材料也可以形成線性方程式(20)。

Ai+1xi+1=bi+1

(20)

由第i層和第i+1層間界面的連接條件表示為

pi+1(yi)=pi(yi),qi+1(yi)=qi(yi)

(21)

式中:pi+1(yi)和pi(yi)分別為第i+1層上端界面和第i層下端界面的位移分量;qi+1(yi)和qi(yi)分別表示第i+1層上端界面和第i層下端界面的應(yīng)力分量。因此,

Ri+1(yi)=Ri(yi)

(22)

聯(lián)立各材料的控制方程,最終得到系統(tǒng)方程式(23)。

Aax=b

(23)

式中:Aa為各層材料的系數(shù)矩陣組成的總系數(shù)矩陣;x為全部的未知分量組成的向量;b為全部已知分量組成的向量。

2.2 精細(xì)積分法

若疊層結(jié)構(gòu)的第i層足夠薄,式(17)中矩陣Ri(y)可表示為(Ri(v1)+Ri(v2))/2。其中,v1和v2分別為第i層材料上、下表面沿y軸的坐標(biāo)。若第i層并非足夠薄,可以將該層材料分成Ki=2k層材料,各個(gè)子層的厚度為Δi=hi/2k,則Ri(y)可以在足夠薄的子層中取平均值來處理。例如,對i層進(jìn)行過2k次劃分之后的第一個(gè)子層在式(14)的左右兩端同時(shí)做積分,可得

(24)

(25)

整理、移項(xiàng)后得

(26)

(27)

對第i層各個(gè)子層進(jìn)行平均近似。同時(shí)注意到各子層之間的連續(xù)條件為

(28)

從第一個(gè)子層迭代到最后一個(gè)子層,可得

(29)

最后和通解法一樣聯(lián)立控制方程,可得到系統(tǒng)方程式(30)。

(30)

3 數(shù)值算例

本算例中所用模型為各向同性材料組成的疊層結(jié)構(gòu)(如圖1),長l=2m,高h(yuǎn)=0.3m,每一層材料的高度均為0.1 m。x=0 m和x=1 m的兩邊均為固支約束。沿y軸正方向的材料參數(shù)分別為彈性模量E1=2.1×1111Pa,E2=1.05×1011Pa,E3=2.1×1111Pa;泊松比均為μ=0.3;熱膨脹系數(shù)α1=1×10-5℃,α2=2×10-5℃,α3=1×10-5℃。疊層結(jié)構(gòu)受到的溫度荷載為ΔQ=-10x2+20x。

由于結(jié)構(gòu)的形狀、性質(zhì)和溫度載荷都是對稱的,所以本文中僅分析左邊一半結(jié)構(gòu)。使用狀態(tài)空間有限元法計(jì)算,并將結(jié)果與有限元Abaqus的數(shù)值結(jié)果作對比。

算例中分別使用二次Lagrange單元離散和二次B樣條單元離散該疊層結(jié)構(gòu)。B樣條離散法取節(jié)點(diǎn)向量為{0,0,0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1,1},控制點(diǎn)的坐標(biāo)為{0,0.025,0.075,0.125,0.175,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1};Lagrange單元離散法用以比較的節(jié)點(diǎn)坐標(biāo)與B樣條單元離散法的控制點(diǎn)坐標(biāo)相同。而有限元方法采用Quad(四邊形)單元?jiǎng)澐志W(wǎng)格,節(jié)點(diǎn)的間距為0.02 m。

進(jìn)行不同求解算法的收斂性和運(yùn)算速度的分析,分別使用二次Lagrange單元離散法和二次B樣條單元離散法進(jìn)行對比。表1為不同求解算法下Lagrange單元離散法計(jì)算第1、2層交界面點(diǎn)P1(0.3,0.2)的位移u的結(jié)果。由表1可以看出,不同算法條件下,Lagrange單元離散法均不需要?jiǎng)澐趾芏嘧訉泳湍苁諗?。隨著子層數(shù)增多,通解法計(jì)算速度明顯變慢,而精細(xì)積分法計(jì)算所需時(shí)間穩(wěn)定在0.01 s左右,速度顯著快于通解法。通解法和精細(xì)積分法結(jié)果均與有限元參考解吻合良好。

表1 Lagrange單元離散法使用不同解法計(jì)算點(diǎn)P1(0.3,0.2)的位移u的結(jié)果Table 1 Results of displacements u at point P1 (0.3,0.2) by Lagrange discretization with different solving schemes

表2為不同求解算法下B樣條單元離散法計(jì)算第1、2層交界面點(diǎn)P2(0.25,0.20)的位移u結(jié)果。使用通解法計(jì)算時(shí),B樣條單元離散法并不能隨著子層數(shù)的增多收斂到一定值,并且計(jì)算結(jié)果相較有限元法有較大差距。而使用精細(xì)積分法計(jì)算時(shí),B樣條單元離散法的計(jì)算結(jié)果能收斂到一個(gè)較為準(zhǔn)確的結(jié)果。同樣精細(xì)積分法計(jì)算的速度顯著快于通解法。

對兩種插值法進(jìn)行比較分析,由于Lagrange單元離散法在單元邊界處僅具有C0連續(xù),而本文所采用的二次B樣條單元離散法的連續(xù)性為C2,所以B樣條單元離散法的結(jié)果較Lagrange單元離散法更為精確。

算例均采用精細(xì)積分法求解狀態(tài)方程,k取值為4,圖2(a)為使用Lagrange單元離散法和B樣條單元離散法計(jì)算得到的第1、2層間界面位移u與有限元計(jì)算結(jié)果的對比。圖2(b)同樣為兩種離散法計(jì)算的相同點(diǎn)的位移v的結(jié)果與有限元法的對比。圖3為使用Lagrange單元離散法和B樣條單元離散法計(jì)算得到的第1層上界面的位移u和v與有限元法的對比。圖2、3中:d均表示節(jié)點(diǎn)和疊層板最左邊的距離,單位為m;FEM表示有限元計(jì)算結(jié)果,SSM-L表示Lagrange離散單元狀態(tài)空間法計(jì)算的結(jié)果;SSM-B表示B樣條離散單元狀態(tài)空間法計(jì)算的結(jié)果。由圖2、3可知:頂層和交界面的橫向位移u均小于0,即疊層板在橫向上受到擠壓后收縮,形變程度先大后小;頂層和交界面的縱向位移v均大于0,即疊層板在縱向上受到擠壓后凸起,形變程度越來越大;兩種插值方法得到的位移都能與有限元法吻合。

圖3 不同離散法頂層表面位移u和vFig.3 Displacements u and v on the top layer by different discretizations

圖4(a)為使用Lagrange單元離散法和B樣條單元離散法計(jì)算得到的第1、2層界面應(yīng)力σxy的結(jié)果與有限元結(jié)果的比較。圖4(b)同樣為兩種離散法計(jì)算的應(yīng)力σyy與有限元結(jié)果的對比。由圖4(a)可知:0≤d≤0.1m則σxy變大;d≥0.1 m后σxy變小。由圖4(b)可知:0≤d≤0.3 m則σyy變化劇烈;而d≥0.3 m后σyy變化較為平穩(wěn),但在固支端應(yīng)力變化較為劇烈。兩種離散方案的數(shù)值結(jié)果均能夠較好與有限元結(jié)果吻合,Lagrange單元離散法計(jì)算σyy(圖4(b))的結(jié)果比B樣條單元離散法的計(jì)算結(jié)果精度差。但B樣條單元離散法在節(jié)點(diǎn)處的連續(xù)性更優(yōu)。

圖4 不同離散法下第1、2層交界面應(yīng)力σxy和σyyFig.4 Stresses σxy and σyy on the interface of 1st and 2nd layer by different discretizations

狀態(tài)空間法求正應(yīng)力σxx需要利用已求得的分量進(jìn)一步計(jì)算。

σxx=D-1(ER-GQ)

(31)

式中:常數(shù)矩陣D、E、G依據(jù)材料參數(shù)確定,即第1層材料的下界面和第2層材料的上界面各節(jié)點(diǎn)受到的正應(yīng)力σxx有差異。

圖5為使用Lagrange單元離散法和B樣條單元離散法計(jì)算的頂層上表面應(yīng)力σxx的結(jié)果。圖6為使用兩種離散方法計(jì)算第1、2層界面處的σxx,并與有限元結(jié)果對比。由圖5可知,頂層面的σxx也是在邊界處的變化即為劇烈,而其他區(qū)域較為平穩(wěn)。

圖5 不同離散法下最頂層面應(yīng)力σxxFig.5 Stresses σxx on the top layer by different discretizations

圖6 不同離散法下第1、2層交界面應(yīng)力σxxFig.6 Stresses σxx on the interface of 1st and 2nd layer by different discretizations

由圖6可知,第1、2層交界面處,不同層所受到的應(yīng)力σxx確實(shí)是不同的,但總體仍是邊界處的變化更為劇烈。由圖5和6可知,Lagrange單元離散法與B樣條單元離散法的計(jì)算結(jié)果均能很好吻合有限元法結(jié)果。

4 結(jié) 論

本文使用狀態(tài)空間法計(jì)算溫度應(yīng)力問題,研究表明:

——相較于傳統(tǒng)的有限元方法,狀態(tài)空間法在劃分單元較少的情況下仍然得到較為精確的結(jié)果;計(jì)算時(shí)即使劃分很多子層,仍然保持較快的計(jì)算速度。

——使用不同離散法,得到的結(jié)果略有差異,在B樣條單元離散法的控制節(jié)點(diǎn)數(shù)量和Lagrange單元離散法的節(jié)點(diǎn)數(shù)量一致的條件下,B樣條單元離散法得到的應(yīng)力結(jié)果更為精確。

——比較精細(xì)積分法和通解法的計(jì)算時(shí)長及精度發(fā)現(xiàn),無論何種單元離散方案,精細(xì)積分解法均優(yōu)于通解法,主要體現(xiàn)在計(jì)算時(shí)長更短且精度更高。雖然兩種方法需要計(jì)算的矩陣維度、大小相同,但得到矩陣的過程中精細(xì)積分法不需要計(jì)算eNi(hi),節(jié)省了大量時(shí)間。