李慶華，陳莘莘

(華東交通大學(xué) 土木建筑學(xué)院，南昌 330013)

當(dāng)結(jié)構(gòu)受到溫變，一般會(huì)產(chǎn)生熱應(yīng)力，并且熱應(yīng)力是物體破壞的一個(gè)重要因素[1-2]。對(duì)受熱結(jié)構(gòu)進(jìn)行分析時(shí)，解耦方法可先由熱傳導(dǎo)方程求出溫度分布，再由熱彈性方程求解位移和應(yīng)力。但是，解耦方法沒有考慮結(jié)構(gòu)變形對(duì)溫度場(chǎng)的影響[3]。事實(shí)上，熱彈性力學(xué)中最基本的問題就是耦合熱彈性問題。在耦合熱彈性問題中，溫度和變形會(huì)相互影響，溫度場(chǎng)和應(yīng)變場(chǎng)的耦合項(xiàng)必須體現(xiàn)在熱傳導(dǎo)方程中。為了求解溫度、位移和應(yīng)力，必須聯(lián)立求解熱傳導(dǎo)方程和熱彈性運(yùn)動(dòng)方程。相對(duì)于非耦合熱彈性問題，耦合熱彈性問題求解更困難。

熱應(yīng)力問題的數(shù)值方法主要基于發(fā)展較為成熟的有限元法[4-5]和邊界元法[6-8]。近年來，部分學(xué)者嘗試采用無網(wǎng)格法[9-12]求解熱應(yīng)力問題。無網(wǎng)格法不僅能夠避免網(wǎng)格生成的復(fù)雜過程，而且在節(jié)點(diǎn)分布不規(guī)則時(shí)，損失的計(jì)算精度較小，從而日益得到重視[13-14]。近十多年來發(fā)展起來的無網(wǎng)格法―無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法[15-16]不僅允許加權(quán)函數(shù)和試函數(shù)取自不同的函數(shù)空間[17]，而且克服了本質(zhì)邊界條件不易施加的困難。此方法中，所有的積分都在中心為所考慮點(diǎn)的多邊形子域上進(jìn)行，而且多邊形子域的構(gòu)造十分方便。目前，無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法在很多領(lǐng)域都得到廣泛應(yīng)用[18-20]。本文采用自然鄰接點(diǎn)插值對(duì)溫度和位移分別插值，與局部加權(quán)余量法結(jié)合，提出了適用于耦合熱彈性動(dòng)力學(xué)問題的無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法。最后，通過數(shù)值算例驗(yàn)證了本文方法應(yīng)用于耦合熱彈性動(dòng)力學(xué)問題分析的有效性和合理性。

1 自然鄰接點(diǎn)插值

在求解域內(nèi)給定M個(gè)離散節(jié)點(diǎn)，其集合為N={x1,x2,…,xM}。對(duì)任一節(jié)點(diǎn)xI，其一階Voronoi結(jié)構(gòu)可定義為

TI={x∈R2:d(x,xI)

(1)

式中：d(x,xI)表示點(diǎn)x與節(jié)點(diǎn)xI之間的距離。

為計(jì)算Sibson插值形函數(shù)，需進(jìn)一歩定義二階Voronoi結(jié)構(gòu)。

TIJ={x∈R2:d(x,xI)

d(x,xK),?J≠I≠K}

(2)

圖1所示為點(diǎn)x的一階和二階Voronoi結(jié)構(gòu)。根據(jù)Sibson插值的定義[15]，點(diǎn)x的插值函數(shù)為

φI(x)=AI(x)/A(x)

(3)

式中：AI(x)為點(diǎn)x的二階Voronoi結(jié)構(gòu)TxI的面積；A(x)為點(diǎn)x的一階Voronoi結(jié)構(gòu)Tx的面積。

圖1 點(diǎn)x的一次和二次Voronoi結(jié)構(gòu)Fig.1 First-order and second-order Voronoi cells about

定義了各節(jié)點(diǎn)的插值函數(shù)后，點(diǎn)x的溫度函數(shù)類似于有限元法可寫為

(4)

式中：θI(I=1,2,L,n)是點(diǎn)x周圍自然鄰接點(diǎn)的節(jié)點(diǎn)溫度。

2 控制方程的弱形式及其離散化

設(shè)線性耦合熱彈性動(dòng)力學(xué)問題的區(qū)域?yàn)棣福棣傅倪吔?。熱彈性力學(xué)的應(yīng)力、位移與溫度之間的關(guān)系為[1]

σij=2μεij+λεkkδij-βθδij

(5)

式中：σij為應(yīng)力;εij為應(yīng)變，λ和μ為拉梅常數(shù);β為熱應(yīng)力系數(shù);θ為溫度變化值。小變形情況下，應(yīng)變?chǔ)舏j與位移ui的幾何關(guān)系為

εij=(ui,j+uj,i)/2

(6)

在經(jīng)典的熱彈性理論中，運(yùn)動(dòng)方程和熱傳導(dǎo)方程可表示為[1]

(7)

(8)

(9)

(10)

(11)

(12)

θ(x,0)=θ0,x∈Ω

(13)

u(x,0)=u0,x∈Ω

(14)

v(x,0)=v0,x∈Ω

(15)

式中:θ0為初始溫度;u0和v0分別為初始位移和初始速度。

(16)

(17)

式中:vI為加權(quán)函數(shù)。對(duì)式(16)左邊積分的第1項(xiàng)進(jìn)行分部積分并利用散度定理后，可得

(18)

(19)

圖2 局部多邊形子域Fig.2 The local polygonal

以此類推，式(17)可變?yōu)?/p>

(20)

由于只對(duì)空間域進(jìn)行離散，求解域Ω內(nèi)的位移u(x,t)和溫度θ(x,t)可由式(4)表示為

(21)

將式(21)代入式(19)和式(20)，可得耦合熱彈性動(dòng)力學(xué)問題的離散控制方程為

(22)

(23)

式中:

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

式中:

(32)

(33)

(34)

(35)

(36)

(37)

(38)

對(duì)于平面應(yīng)變問題，需把E換成E/(1-v2)，v換成v/(1-v)，α換成(1+v)α。式(22)和式(23)可合并為

(39)

式(39)即為施加邊界條件后的系統(tǒng)耦合微分方程組，可采用Newmark方法[21]進(jìn)行求解。

3 數(shù)值算例

為了驗(yàn)證所提方法的有效性，考慮如圖3所示的單位面積方板，該問題為平面應(yīng)變狀態(tài)下的一個(gè)經(jīng)典算例。初始時(shí)刻板的溫度和位移均為0，板上邊受到突加的熱載荷，另外3邊均絕熱，且無法向位移。彈性模量E=1，泊松比v=0.3，導(dǎo)熱系數(shù)k=1，密度ρ=1，比熱容c=1，熱膨脹系數(shù)α=0.02。計(jì)算中，采用15×15個(gè)節(jié)點(diǎn)將方板離散，時(shí)間步長(zhǎng)取為2.0×10-3。

圖3 突加熱載荷的方板Fig.3 A suddenly heated unit square

當(dāng)不考慮慣性項(xiàng)和耦合項(xiàng)時(shí)，此問題屬于準(zhǔn)靜態(tài)熱彈性力學(xué)問題。圖4和圖5分別給出了方板y軸上不同坐標(biāo)處的溫度和豎向位移變化情況。從圖4和圖5可以看出，本文數(shù)值解與解析解[22]吻合得很好。

圖4 溫度隨時(shí)間的變化Fig.4 Temporal variation of the

圖5 豎向位移隨時(shí)間的變化Fig.5 Temporal variation of the vertical

當(dāng)考慮慣性項(xiàng)時(shí)，為了便于對(duì)耦合和非耦合情況下的計(jì)算結(jié)果進(jìn)行比較，引入如下修正的耦合系數(shù)[12]

(40)

式中:耦合系數(shù)η的取值范圍一般為0.01～0.2。相關(guān)溫度取為θ0=100，則對(duì)應(yīng)的耦合系數(shù)為η=0.186。圖6和圖7分別為y軸上不同坐標(biāo)處耦合項(xiàng)對(duì)溫度和位移的影響。顯然，耦合項(xiàng)對(duì)溫度的影響很大，但對(duì)位移的影響可忽略不計(jì)。

圖6 耦合效應(yīng)對(duì)溫度的影響Fig.6 Coupling effects on the

圖7 耦合效應(yīng)對(duì)豎向位移的影響Fig.7 Coupling effects on the vertical displacement

4 結(jié)論

無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法是一種簡(jiǎn)單適用，且效率和精度均十分優(yōu)良的無網(wǎng)格分析方法。在采用自然鄰接點(diǎn)插值對(duì)位移和溫度分別插值的基礎(chǔ)上，將FEM三角形線性單元的形函數(shù)作為加權(quán)函數(shù)，采用加權(quán)余量法詳細(xì)推導(dǎo)了二維耦合熱彈性動(dòng)力學(xué)問題的無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法計(jì)算公式。數(shù)值算例表明，所提的二維耦合熱彈性動(dòng)力學(xué)問題求解方法可行。