二維耦合熱彈性動力學(xué)問題的無網(wǎng)格自然鄰接點Petrov-Galerkin法

李慶華，陳莘莘

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

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

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

1 自然鄰接點插值

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

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

(1)

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

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

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

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

(2)

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

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

(3)

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

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

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

(4)

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

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

設(shè)線性耦合熱彈性動力學(xué)問題的區(qū)域為Ω，Γ為Ω的邊界。熱彈性力學(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)變εij與位移ui的幾何關(guān)系為

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

(6)

在經(jī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ù)。對式(16)左邊積分的第1項進(jìn)行分部積分并利用散度定理后，可得

(18)

(19)

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

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

(20)

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

(21)

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

(22)

(23)

式中:

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

式中:

(32)

(33)

(34)

(35)

(36)

(37)

(38)

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

(39)

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

3 數(shù)值算例

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

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

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

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

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

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

(40)

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

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

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

4 結(jié)論

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