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

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

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

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

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

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

同時,本文分別采用通解法和精細積分法求解狀態(tài)空間方程,以研究何種方法更適用于固支問題。本文數(shù)值算例表明,位移的數(shù)值結(jié)果在2種離散方案及2種求解方案情況下均較為精準,而在計算應力時,B樣條結(jié)果精度更高。雖然精細積分法和通解法均能較快完成二維狀態(tài)空間方程的計算,其中精細積分法由于采用高效迭代的技術(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è)均受固支約束。由于材料的幾何形狀及溫度荷載均對稱,故只截取左半邊進行分析。邊界條件:x=0的邊的位移u=v=0;x=l/2的邊的位移u=0,應力σ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]變分原理,以位移和應力滿足邊界條件的假設對疊層結(jié)構(gòu)中的任意一層建立二維彈性力學中的平衡方程式(1)和幾何方程式(2)。

(1)

(2)

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

Lagrange插值公式為

(3)

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

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

(4)

(5)

(6)

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

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

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

(14)

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

(15)

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

將式(15)計算所得的關(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層下端界面的應力分量。因此,

Ri+1(yi)=Ri(yi)

(22)

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

Aax=b

(23)

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

2.2 精細積分法

若疊層結(jié)構(gòu)的第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層各個子層進行平均近似。同時注意到各子層之間的連續(xù)條件為

(28)

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

(29)

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

(30)

3 數(shù)值算例

本算例中所用模型為各向同性材料組成的疊層結(jié)構(gòu)(如圖1),長l=2m,高h=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)空間有限元法計算,并將結(jié)果與有限元Abaqus的數(shù)值結(jié)果作對比。

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

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

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

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

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

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

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

σxx=D-1(ER-GQ)

(31)

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

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

圖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樣條單元離散法的計算結(jié)果均能很好吻合有限元法結(jié)果。

4 結(jié) 論

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

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

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

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