二維耦合熱彈性動力學問題的無網格自然鄰接點Petrov-Galerkin法

李慶華，陳莘莘

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

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

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

1 自然鄰接點插值

在求解域內給定M個離散節點，其集合為N={x1,x2,…,xM}。對任一節點xI，其一階Voronoi結構可定義為

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

(1)

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

為計算Sibson插值形函數，需進一歩定義二階Voronoi結構。

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

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

(2)

圖1所示為點x的一階和二階Voronoi結構。根據Sibson插值的定義[15]，點x的插值函數為

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

(3)

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

圖1 點x的一次和二次Voronoi結構Fig.1 First-order and second-order Voronoi cells about

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

(4)

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

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

設線性耦合熱彈性動力學問題的區域為Ω，Γ為Ω的邊界。熱彈性力學的應力、位移與溫度之間的關系為[1]

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

(5)

式中：σij為應力;εij為應變，λ和μ為拉梅常數;β為熱應力系數;θ為溫度變化值。小變形情況下，應變εij與位移ui的幾何關系為

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

(6)

在經典的熱彈性理論中，運動方程和熱傳導方程可表示為[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為加權函數。對式(16)左邊積分的第1項進行分部積分并利用散度定理后，可得

(18)

(19)

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

以此類推，式(17)可變為

(20)

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

(21)

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

(22)

(23)

式中:

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

式中:

(32)

(33)

(34)

(35)

(36)

(37)

(38)

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

(39)

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

3 數值算例

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

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

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

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

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

當考慮慣性項時，為了便于對耦合和非耦合情況下的計算結果進行比較，引入如下修正的耦合系數[12]

(40)

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

圖6 耦合效應對溫度的影響Fig.6 Coupling effects on the

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

4 結論

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