考慮溫度梯度影響的偶應(yīng)力理論的無網(wǎng)格方法

谷小強， 常利武， 孫玉周， 張旺璽

(1.中原工學(xué)院 建筑工程學(xué)院， 河南 鄭州 450007； 2.中原工學(xué)院 材料與化工學(xué)院， 河南 鄭州 450007)

隨著科技的進步和問題研究的深入，許多原來能滿足精度要求的理論由于限制條件的增加而不再能滿足，進而發(fā)展為新的理論[1]。如伴隨著熱以有限速度傳播的熱機耦合理論的提出，經(jīng)典熱彈性理論發(fā)展為廣義熱彈性理論[2-3]；考慮材料的微結(jié)構(gòu)特征對熱彈性介質(zhì)的影響，在控制方程中引入偶應(yīng)力張量，廣義熱彈性理論發(fā)展為熱彈性偶應(yīng)力理論[4]。新近大量的研究表明，熱彈性介質(zhì)材料在超高溫度梯度、超高溫升速率等特殊環(huán)境中表現(xiàn)出復(fù)雜的熱力耦合特征，溫度梯度的存在必然對溫度梯度固體的應(yīng)力場產(chǎn)生影響[5]。由于溫度梯度會使物質(zhì)結(jié)構(gòu)產(chǎn)生偶應(yīng)力，故可將溫度梯度產(chǎn)生的偶應(yīng)力引入到熱彈性本構(gòu)關(guān)系中,以研究其對結(jié)構(gòu)力學(xué)性能的影響。

熱彈性偶應(yīng)力理論雖然考慮了應(yīng)變梯度對結(jié)構(gòu)的影響，但并沒有考慮溫度梯度對應(yīng)力場的影響。本文將溫度梯度引入到偶應(yīng)力熱彈性理論的本構(gòu)關(guān)系中，建立了考慮溫度梯度影響的無網(wǎng)格數(shù)值計算方法，并編寫了相應(yīng)的Fortran計算程序，對溫度梯度影響下的結(jié)構(gòu)性能進行了研究。

1 無網(wǎng)格法計算框架

1.1 本構(gòu)方程

若只考慮穩(wěn)態(tài)熱傳導(dǎo)問題，可以認為溫度梯度只與偶應(yīng)力有關(guān)。將溫度梯度引入到偶應(yīng)力理論的本構(gòu)關(guān)系中，即可得到考慮溫度梯度影響的偶應(yīng)力理論的本構(gòu)方程。對于各向同性的線彈性熱傳導(dǎo)固體，Helmholtz自由能函數(shù)Φ可表示為應(yīng)變張量、應(yīng)變梯度張量、溫度以及溫度梯度的函數(shù)[4-5]

(1)

式中：λ、μ是拉梅常數(shù)；η、η,i為偶應(yīng)力介質(zhì)的材料常數(shù)；β是熱力耦合系數(shù)，β=(3λ+2μ)α；α是材料的線性膨脹系數(shù)；εji為應(yīng)變張量；χji為應(yīng)變梯度張量；ΔT表示溫度的變化量，ΔT=T-T0；c是材料的比熱容。

參考張曉敏等提出的方法[5]，考慮溫度梯度影響的偶應(yīng)力理論的本構(gòu)方程可表示為：

(2)

(3)

式中,γ為介質(zhì)常量。本文取γ=1.68×10-4N/K。

1.2 無網(wǎng)格法計算框架

與傳統(tǒng)數(shù)值計算方法如有限元方法、邊界元法相比，無網(wǎng)格方法具有許多突出的優(yōu)點，如移動最小二乘近似的形函數(shù)具有高階連續(xù)的特征[6-8]。利用形函數(shù)的高階連續(xù)特征，溫度梯度、結(jié)構(gòu)的應(yīng)變和旋轉(zhuǎn)曲率可直接由節(jié)點的溫度或位移的插值得到，這樣，問題大大簡化。

在無熱源的情況下，由于溫度發(fā)生改變，彈性結(jié)構(gòu)體內(nèi)部產(chǎn)生的勢能由結(jié)構(gòu)的彈性應(yīng)變能提供。系統(tǒng)的彈性應(yīng)變能為：

(4)

其中:

[εT]T=[αΔTαΔT0αT,iαT,i]

(5)

[εE]T=[ε11ε22γ12χ1χ2]

(6)

[ε]T=[εE]+[εT]

(7)

式中，[D]為剛度矩陣。

問題域上的位移函數(shù)u(x)可表示為

u(x)=uh(x)=φ(x)u

(8)

其中:φ(x)為形函數(shù)，u為節(jié)點的位移向量。

應(yīng)變分量ε可寫成：

(9)

uT=[u1v1u2v2…unvn]

(10)

利用移動最小二乘近似，溫度場可表示為：

T(x)=Th(x)=φ(x)T

(11)

溫度梯度可表示為：

(12)

TT=[T1xT1yT2xT2y…TnxTny]

(13)

將式(13)、式(9)和式(5)代入式(4),可得

(15)

其中：

(16)

(17)

(18)

式中，C只與溫度有關(guān)，不產(chǎn)生勢能，可略去。式(15)可簡化為

(19)

令δΠ=0，得無網(wǎng)格法的離散控制方程為

[K][u]=[F]

(20)

式中：

(21)

(22)

2 算例及結(jié)果分析

圖1所示為矩形薄板，長度L=20 cm，高度H=5 cm，厚度t=0.01 cm，材料的彈性模量E=2.6 GPa，泊松比ν=0.32，線性膨脹系數(shù)α=1.0×10-5/K。對薄板施加的連續(xù)溫度場的邊界條件分別取TL=0，TU=10°和TL=0，TU=40°兩種工況。數(shù)值計算中，薄板模型的總節(jié)點分布為6×21，背景網(wǎng)格采用3×3的高斯積分方案，且溫度場和位移場采用相同的節(jié)點分布和背景網(wǎng)格積分方案。

圖1 薄板模型

考慮溫度梯度影響的偶應(yīng)力熱彈性力學(xué)問題的研究是在偶應(yīng)力的基礎(chǔ)上研究溫度梯度對結(jié)構(gòu)力學(xué)性能的影響，故還需要考慮微結(jié)構(gòu)的尺度效應(yīng)問題。本文通過含節(jié)理特征的層狀巖體的尺度因子(l=0.4t)來表述微結(jié)構(gòu)的影響[9-10]。

圖2和圖3分別為偶應(yīng)力熱彈性力學(xué)下的薄板彎曲變形圖和考慮溫度梯度影響的偶應(yīng)力理論下的薄板彎曲變形圖。從圖中可以看出，兩種情況下結(jié)構(gòu)內(nèi)部的彎曲程度從邊緣到中間均不斷加大，且考慮溫度梯度影響的薄板的彎曲變形得到了加強。

圖2 偶應(yīng)力理論的熱彈性位移云圖(TL=0，TU=10 ℃)

表1所示為邊界條件TL=0，TU=10°和TL=0，TU=40°兩種工況下薄板上邊緣5個節(jié)點的Y方向位移值(括號里的數(shù)值表示邊界條件為TL=0，TU=40°時的位移值)。由表1可知，薄板彎曲變形從邊緣部位向薄板中間隆起并在薄板中間達到峰值，考慮溫度梯度影響的偶應(yīng)力理論下得到的位移峰值較偶應(yīng)力熱彈性理論下的結(jié)果，在兩種溫度梯度下分別提高了2.5%和5.8%。這說明溫度梯度的存在增加了薄板的彎曲變形值，其提高幅度隨著溫度梯度的增大而增大。

圖3 溫度梯度偶應(yīng)力位移云圖(TL=0，TU=10 ℃)

表1 Y方向位移值

3 結(jié)論

在假設(shè)溫度梯度只產(chǎn)生偶應(yīng)力的基礎(chǔ)上，將其引入到偶應(yīng)力熱彈性理論的本構(gòu)方程中，得到了一個簡化的考慮溫度梯度影響的偶應(yīng)力本構(gòu)方程，建立了無網(wǎng)格數(shù)值計算框架，并編寫了相應(yīng)的Fortran計算程序，對溫度梯度影響下薄板的彎曲變形規(guī)律進行了研究。研究結(jié)果表明，利用移動最小二乘近似建立起來的無網(wǎng)格方法可以較好地模擬考慮溫度梯度和應(yīng)變梯度的高階連續(xù)問題，薄板在考慮溫度梯度影響時有顯著的變形，且變形值隨著溫度梯度的增大而增大。本文結(jié)論為考慮應(yīng)變梯度場和溫度梯度場影響的結(jié)構(gòu)的數(shù)值模擬提供了一種新的思路。