基于分數階理論的兩端固定桿熱彈動力響應

郭 穎,杜 翠,馬建軍

(1.河南科技大學 a.土木建筑學院;b.機電工程學院,河南 洛陽 471023;2.成都信息工程大學 軟件工程學院,四川 成都 610000)

0 引言

Biot[1]提出的耦合熱彈性理論存在2個不足,一是該理論忽略了溫度對彈性變化的影響;二是由于熱方程是拋物線形的,其預測的熱波傳播速度是無限大的,這與事實不符。為了彌補這些不足,學者們提出了不同的熱彈性理論,主要有:Lord-Shulman (L-S)[2]廣義熱彈性理論、Green-Lindsay (G-L)[3]廣義熱彈性理論以及Green-Naghdi (G-N)[4-6]廣義熱彈性理論?；谏鲜隼碚?學者們進行了大量研究[7-13],分析了彈性、黏彈性介質以及功能梯度材料中多維熱彈應力應變問題,探討了材料物性參數以及外荷載作用對各物理量的影響。

移動熱源在冶金過程中充當著重要的角色,特別是在焊接過程中。文獻[14-15]基于廣義熱彈性理論,研究了考慮移動熱源影響的兩端固定桿的熱彈/電磁熱彈問題。文獻[16-17]考慮移動熱源的影響,分析了中空圓柱的多物理場耦合動力響應問題。文獻[18]研究了二維軸對稱無限圓柱體在移動熱源影響下的熱彈耦合問題。文獻[19]分析了運動熱源和簡諧熱作用下無限長中空圓柱的熱彈性耦合響應。文獻[20]對柱體類工件在移動熱源作用下溫度場特性進行了分析。上述研究雖然考慮了移動熱源,但沒有考慮材料本身的特性以及溫度對材料特性的影響。

通常情況下材料的彈性模量、熱傳導率等特性參數會與溫度有一定的相關性,從而影響介質的熱彈力學行為。文獻[21]研究了具有溫度相關材料性能的功能梯度復合材料的熱應力問題。文獻[22]基于G-L廣義熱彈性理論,研究了熱沖擊作用下材料特性參數與溫度相關特性之間的關系。文獻[23]基于G-N理論,研究了平面波在各向同性且與溫度相關的半無限大介質中傳播的問題。文獻[24]基于三相滯后熱彈性模型分析了熱沖擊作用下平面波在彈性半空間自由表面與溫度相關的反射問題。文獻[25]建立了可預測納米結構在極端環境下的熱彈性行為的非局部熱彈性模型,分析了瞬態熱對非局部參數和溫度相關導熱系數的影響。

Abel將分數階積分算子應用于求解等時曲線的積分問題中[26],隨后,學者們發現分數階與整數階具有一定的自相似性,分數階廣義熱彈性理論得以快速發展。常用的分數階熱彈性理論有:Sherief[27]型分數階廣義熱傳導理論、Youssef型分數階廣義熱傳導理論[28]、Ezzat型分數階雙相滯后廣義熱傳導理論[29]和Ezzat型分數階三相滯后廣義熱傳導理論[30]。文獻[31]推導出了考慮2個不同熱松弛時間因子的分數階廣義熱彈性理論。文獻[32-33]分析了分數階理論下彈性/黏彈性介質的動力響應問題。文獻[34]采用Sherief型分數階廣義熱傳導理論研究了受移動熱源作用的兩端固定桿熱彈耦合問題。文獻[35-36]基于不同分數階廣義熱彈性理論,研究了多種材料特性參數與溫度相關的熱彈/電磁熱彈耦合動力響應問題。文獻[37-38]考慮了受移動熱源作用的三維熱彈耦合問題,了解了熱源移動速度對耦合問題的影響。文獻[39-41]基于不同分數階理論,結合移動熱源以及多種復雜邊界條件,研究了3種桿件的多物理場耦合動力響應。文獻[42]提出了考慮材料記憶依賴效應和空間非局部效應的非局部廣義熱彈擴散耦合理論。

上述研究中,已有采用不同的分數階積分算子分析介質的熱彈動力響應問題和材料特性對介質影響的研究,但同時考慮分數階積分和2個熱松弛時間因子分析材料特性與溫度相關介質的熱彈耦合問題尚不多見,特別是將2個熱松弛時間因子同時引入分數階廣義熱彈性理論對問題進行求解更為少見,而分數階積分算子的引入可以更為清晰地描述熱波在介質中的傳播情況,消除原本廣義熱彈性理論難以描述的問題。為了能更好地分析不同參數對各無量綱物理量的影響,本文將Ezzat型分數階廣義熱彈性理論和G-L廣義熱彈性理論相結合,建立同時考慮分數階參數和2個熱松弛時間因子的分數階廣義熱彈動力理論,基于該理論建立可描述兩端固定桿溫度場、應力場和變形場等多物理場耦合的數學模型,借助拉普拉斯(Laplace)積分變換及其數值反變換的方法得到無量綱位移、應力和溫度的分布規律,分析分數階參數、2個不同的熱松弛時間因子、溫度相關性參數以及移動熱源傳播速度對無量綱位移、應力和溫度的影響。

1 基本控制方程

基于Ezzat型分數階廣義熱彈性理論和G-L廣義熱彈性理論的基本控制方程[31,34]為:

(1)

其中:σij為應力分量;λ和G為拉梅(Lamé)常數;β1=(3λ+2G)αt,αt為線性熱膨脹系數;εij為應變分量;τ0為第1個熱松弛時間因子;θ=T-T0,T為絕對溫度,T0為參考溫度;δij為克羅內克(Kronecker)記號;α為分數階參數且0<α≤1。

(2)

其中:逗號后的下標表示對坐標求導;字母上方的點表示對時間求導;ρ為質量密度;ui為位移分量。

(3)

(4)

其中:qi為熱流向量分量;η為熵密度;Q為熱源。

(5)

其中:CE為比熱。

(6)

其中:κ為熱傳導系數;τ1為另1個熱松弛時間因子。

上述方程中考慮了如下的溫度相關性參數λ=λ0f(T),G=G0f(T),κ=κ0f(T),β1=β10f(T),其中,λ0,G0,β10和κ0均為與溫度無關的特性參數,f(T)為與溫度相關的一維函數,當f(T)=1時,則表明與溫度不相關[43]:

f(T)=1-ζT。

(7)

f(T)≈1-ζT0,

(8)

其中:ζ為與溫度相關的經驗常數。

2 問題描述

基于Ezzat型分數階廣義熱彈性理論和G-L廣義熱彈性理論,研究一維兩端固定桿受到移動熱源影響的熱彈性瞬態響應問題。對于一維問題,沿桿軸線方向建立一維坐標系,則有非零位移分量ux=u(x,t),方程(1)～方程(5)可簡化為如下形式:

(9)

結合方程(2)和方程(9),可得:

(10)

考慮分數階參數影響的廣義能量方程:

(11)

為了后續計算方便,引入如下無量綱量:

(12)

利用上述無量綱量,對方程(9)～方程(11)進行無量綱化,為了方便起見,將無量綱化后方程中字母上標的星號去掉,可得:

(13)

(14)

(15)

問題的初始條件:

(16)

問題的邊界條件(假設桿的長度為l且兩端絕熱):

(17)

桿在移動熱源的作用下需要考慮熱源移動速度和熱量,其中,υ為熱源移動速度,是1個沿x軸正方向移動的常數,經過無量綱得:

Q=Q0δ(x-υt),

(18)

其中:Q0為移動熱源大小;δ為delta函數。

3 拉普拉斯域中問題求解

引入拉普拉斯變換公式:

(19)

其中:s為拉普拉斯變換中的變換因子。

利用上述拉普拉斯變換公式,方程(13)～方程(15)可變為:

(20)

(21)

(22)

對邊界條件進行拉普拉斯變換得:

(23)

(24)

方程(24)的通解為:

(25)

其中:Ci(i=1,2,3,4)是與變換因子s相關的1組待定系數。C5的表達式如下:

(26)

其中:±k1和±k2是方程(27)的特征根。

k4-m1k2+m2=0,

(27)

可表達為:

(28)

(29)

方程(29)的通解為:

(30)

其中:Cii(i=1,2,3,4)是與變換因子s相關的另1組待定系數。

將方程(25)和方程(30)代入方程(21),可得到如下關系:

(31)

根據邊界條件方程(23),可得:

C1+C2+C3+C4=-C5;

(32)

C1e-k1l+C2ek1l+C3e-k2l+C4ek2l=-C5e-(s/υ)l;

(33)

-C11k1+C22k1-C33k2+C44k2=(s/υ)C55;

(34)

-C11k1e-k1l+C22k1ek1l-C33k2e-k2l+C44k2ek2l=(s/υ)C55e-(s/υ)l。

(35)

聯立方程(32)～方程(35),可得:

(36)

將方程(36)代入方程(25),可得:

(37)

結合方程(31)和方程(36),可得:

(38)

將方程(38)代入方程(30),可得:

(39)

將方程(37)和方程(39)代入方程(20),可得:

(40)

4 數值反變換

(41)

其中:Re為實部;i為虛數單位。為了能加速收斂,結合大量已有試驗和現有計算結果,β應滿足βt≈4.7[44]。

5 數值結果及分析

為了進行數值計算,除了引入方程(40)進行數值拉普拉斯反變換,得到時域內桿的無量綱位移、應力和溫度,還需要引入材料特性相關的參數:λ=7.76×1010Nm-2,G=3.86×1010Nm-2,ρ=8 954 kgm-3,αt=1.78×10-5K-1,CE=383.1 J kg-1K-1,κ=386 Wm-1K-1,T0=293 K,ζ=0.000 5 K-1,Q0=10,l=10。

數值計算得到了桿中無量綱位移、應力和溫度的分布規律。計算中,研究了分數階參數α、2個不同的熱松弛時間因子τ0和τ1、溫度相關性參數?以及移動熱源傳播速度υ對無量綱位移、應力和溫度分布規律的影響效應,數值計算分析分別考慮4種不同的情形,每種情況分數階參數均為α=0.25和α=1.0。4種不同的情形如下:(1)τ1=0.05時,τ0分別取值τ0=0.03,τ0=0.06,τ0=0.09;(2)τ0=0.03時,τ1分別取值τ1=0.05,τ1=0.1,τ1=0.15;(3)τ0=0.03和τ1=0.05時,?分別取值?=0.5,?=1.0,?=1.5;(4)τ0=0.03和τ1=0.05時,υ分別取值υ=1.0,υ=2.0,υ=3.0。所得結果如圖1～圖4所示。

(a) 分數階參數下τ0變化對無量綱位移的影響

圖1中給出了分數階參數和熱松弛時間因子τ0變化對無量綱位移、應力和溫度的影響。從圖1a和圖1b可以看出,無量綱位移和應力隨分數階參數增大逐漸減小,而隨熱松弛時間因子的增加而增大。隨著熱源的移動,圖1a中桿產生了沿桿軸向的熱膨脹變形,熱擾動區域基本保持不變。隨著分數階增大,熱松弛時間因子對位移和應力產生的影響逐漸不明顯。圖1c中僅可以在曲線峰值出看出分數階參數變化對無量綱溫度的影響,但熱松弛時間因子τ0變化對無量綱溫度的影響不明顯。

圖2給出了分數階參數和熱松弛時間因子τ1變化對無量綱位移、應力和溫度的影響。從圖2中可以看出:當分數階參數α=1.0時,熱松弛時間因子τ1變化對各無量綱物理量均有一定影響,但影響不明顯。圖2a形式分數階參數α=0.25時,熱松弛時間因子τ1變化對無量綱位移影響也比較小。從圖2b和圖2c中可以看出:當分數階參數α=0.25時,無量綱應力和溫度曲線峰值和谷值處,熱松弛時間因子τ1變化影響明顯,隨著熱松弛時間因子τ1增大,無量綱應力和溫度逐漸增大。

(a) 不同τ1和α時無量綱位移分布情況

圖3中給出了分數階參數和溫度相關性參數?變化對無量綱位移、應力和溫度的影響。除無量綱溫度外的所有物理量受分數階參數和溫度相關性參數均非常明顯。圖3a中,隨著溫度相關性參數?增大,位移逐漸減小,其峰值逐漸向著z=0處移動,且產生的熱擾動區域也在逐漸減小。圖3b顯示除了熱擾動區域逐漸減小外,隨著溫度相關性參數增大曲線在z=0處明顯減小。從圖3c中可以看出:分數階參數和溫度相關性參數對溫度的影響不大,隨著分數階參數增大,溫度相關性參數對無量綱溫度的影響更不明顯了。

(a) 考慮α影響時?變化對無量綱位移的影響

圖4給出了分數階參數和移動熱源傳播速度υ變化對無量綱位移、應力和溫度的影響。與圖1～圖3類似的是,分數階參數對無量綱溫度影響不明顯。圖4a和圖4b顯示隨著移動熱源傳播速度υ增大,無量綱位移和應力曲線逐漸減小,但擾動區域逐漸增大,且峰值向后移動。從圖4c可以看出:隨熱源移動速度增大,無量綱溫度逐漸減小,主要是因為在相同的時間段內,熱源釋放出相同的熱量,然而隨著熱源速度的不斷增加,單位長度所釋放出的熱量密度降低造成的。

(a) 分數階參數下不同速度對無量綱位移的影響

此外,本文公式和程序中取τ0=0和ζ=0時,可將模型退化成與參考文獻[34]相同的模型,并將本文退化結果與文獻[34]的相應結果進行對比。結果見圖5。圖5a～圖5c表明:當τ0=0、ζ=0和α=0.25時,可得到與文獻[34]吻合較好的計算結果,有微小差異是由于數值計算方法的不同引起,其中最大偏差為9.98%,平均偏差為1.8%。從而進一步驗證了本文計算結果的正確性和合理性。

(a) 無量綱位移與文獻[34]結果的比較

6 結論

本文基于Ezzat型分數階廣義熱彈性理論和G-L廣義熱彈性理論,借助拉普拉斯積分變換及其反變換法研究了移動熱源作用下兩端固定的溫度相關性桿一維熱彈動力響應問題,得到了無量綱位移、應力和溫度沿x方向的分布情況,研究了分數階參數、2個不同的熱松弛時間因子τ0和τ1、溫度相關性參數以及移動熱源傳播速度對各無量綱量的影響。結果表明:

(1)分數階參數對各物理量均有一定的影響,但對無量綱位移和應力的影響更明顯,分數階參數會對無量綱量產生一定的滯后影響。

(2)熱松弛時間因子τ0對無量綱溫度外的其他2個物理量影響顯著,而熱松弛時間因子τ1僅對無量綱溫度的影響最為顯著,分數階參數增大使得2個熱松弛時間因子對3個無量綱物理量的影響逐漸減小。

(3)溫度相關性參數和移動熱源傳播速度對3個無量綱量影響非常明顯,不僅影響曲線峰值的大小,對桿的擾動區域影響也十分明顯。