狀態空間有限元法分析固支疊層結構熱應力

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

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

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

受有限元思想的啟發,Sheng等[7]提出狀態空間有限元法,此方法在疊層結構的水平方向采用有限元法離散位移、應力分量,而在垂直方向仍然采用狀態空間法求解;Wu等[8]提出沿疊層結構的水平方向采用樣條插值法,而在垂直方向使用狀態空間法,這也可以較為準確地計算出疊層結構的物理量。

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

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

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

1 二維模型變分原理

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

圖1 Lagrange單元和B樣條單元離散疊層結構Fig.1 Discretization of the laminated structure by Lagrange and B-spline elements

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

(1)

(2)

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

Lagrange插值公式為

(3)

式中:w為lagrange公式局部坐標系的坐標,wi為第i個lagrange節點的坐標zi為對應的函數值;pk(w)為第k個基函數,在wk處值為1,其余插值點處值為0;Bk為Lagrange插值法產生的所有節點的集合;Ln(w)為總體插值函數,其經過每一個插值點。

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

(4)

(5)

(6)

式中:wi+1和wi+p分別表示第i+1和第i+p個控制節點的坐標;Ni,0為0階B樣條的第i個基函數;Ni,p表示p階B樣條的第i個基函數。由式(4)和(5)可知,p階B樣條基函數可由0階基函數迭代得到。

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

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

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

ε=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)進行合并得到式(11)。

(10)

(11)

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

Dfsf=EfRf-GfK

(12)

(13)

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

(14)

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

(15)

式中:K1是節點1處的溫度。

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

(16)

2 狀態方程的求解

2.1 通解法

對于疊層結構中的第i層,得到如式(16)的非齊次狀態方程為

(17)

式中:Ri為狀態方程中的狀態向量,是僅關于y的函數;Ni為狀態方程中的狀態矩陣;Mi為狀態方程中的輸入矩陣。

根據狀態方程的通解,并令y=yi,可得:

(18)

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

Aixi=bi

(19)

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

同樣的,第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層下端界面的應力分量。因此,

Ri+1(yi)=Ri(yi)

(22)

聯立各材料的控制方程,最終得到系統方程式(23)。

Aax=b

(23)

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

2.2 精細積分法

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

(24)

(25)

整理、移項后得

(26)

(27)

對第i層各個子層進行平均近似。同時注意到各子層之間的連續條件為

(28)

從第一個子層迭代到最后一個子層,可得

(29)

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

(30)

3 數值算例

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

由于結構的形狀、性質和溫度載荷都是對稱的,所以本文中僅分析左邊一半結構。使用狀態空間有限元法計算,并將結果與有限元Abaqus的數值結果作對比。

算例中分別使用二次Lagrange單元離散和二次B樣條單元離散該疊層結構。B樣條離散法取節點向量為{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},控制點的坐標為{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單元離散法用以比較的節點坐標與B樣條單元離散法的控制點坐標相同。而有限元方法采用Quad(四邊形)單元劃分網格,節點的間距為0.02 m。

進行不同求解算法的收斂性和運算速度的分析,分別使用二次Lagrange單元離散法和二次B樣條單元離散法進行對比。表1為不同求解算法下Lagrange單元離散法計算第1、2層交界面點P1(0.3,0.2)的位移u的結果。由表1可以看出,不同算法條件下,Lagrange單元離散法均不需要劃分很多子層就能收斂。隨著子層數增多,通解法計算速度明顯變慢,而精細積分法計算所需時間穩定在0.01 s左右,速度顯著快于通解法。通解法和精細積分法結果均與有限元參考解吻合良好。

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

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

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

算例均采用精細積分法求解狀態方程,k取值為4,圖2(a)為使用Lagrange單元離散法和B樣條單元離散法計算得到的第1、2層間界面位移u與有限元計算結果的對比。圖2(b)同樣為兩種離散法計算的相同點的位移v的結果與有限元法的對比。圖3為使用Lagrange單元離散法和B樣條單元離散法計算得到的第1層上界面的位移u和v與有限元法的對比。圖2、3中:d均表示節點和疊層板最左邊的距離,單位為m;FEM表示有限元計算結果,SSM-L表示Lagrange離散單元狀態空間法計算的結果;SSM-B表示B樣條離散單元狀態空間法計算的結果。由圖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樣條單元離散法計算得到的第1、2層界面應力σxy的結果與有限元結果的比較。圖4(b)同樣為兩種離散法計算的應力σyy與有限元結果的對比。由圖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變化較為平穩,但在固支端應力變化較為劇烈。兩種離散方案的數值結果均能夠較好與有限元結果吻合,Lagrange單元離散法計算σyy(圖4(b))的結果比B樣條單元離散法的計算結果精度差。但B樣條單元離散法在節點處的連續性更優。

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

狀態空間法求正應力σxx需要利用已求得的分量進一步計算。

σxx=D-1(ER-GQ)

(31)

式中:常數矩陣D、E、G依據材料參數確定,即第1層材料的下界面和第2層材料的上界面各節點受到的正應力σxx有差異。

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

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

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

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

4 結 論

本文使用狀態空間法計算溫度應力問題,研究表明:

——相較于傳統的有限元方法,狀態空間法在劃分單元較少的情況下仍然得到較為精確的結果;計算時即使劃分很多子層,仍然保持較快的計算速度。

——使用不同離散法,得到的結果略有差異,在B樣條單元離散法的控制節點數量和Lagrange單元離散法的節點數量一致的條件下,B樣條單元離散法得到的應力結果更為精確。

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