基于分?jǐn)?shù)階理論的兩端固定桿熱彈動力響應(yīng)

郭 穎,杜 翠,馬建軍

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

0 引言

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

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

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

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

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

1 基本控制方程

基于Ezzat型分?jǐn)?shù)階廣義熱彈性理論和G-L廣義熱彈性理論的基本控制方程[31,34]為:

(1)

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

(2)

其中:逗號后的下標(biāo)表示對坐標(biāo)求導(dǎo);字母上方的點表示對時間求導(dǎo);ρ為質(zhì)量密度;ui為位移分量。

(3)

(4)

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

(5)

其中:CE為比熱。

(6)

其中:κ為熱傳導(dǎo)系數(shù);τ1為另1個熱松弛時間因子。

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

f(T)=1-ζT。

(7)

f(T)≈1-ζT0,

(8)

其中:ζ為與溫度相關(guān)的經(jīng)驗常數(shù)。

2 問題描述

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

(9)

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

(10)

考慮分?jǐn)?shù)階參數(shù)影響的廣義能量方程:

(11)

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

(12)

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

(13)

(14)

(15)

問題的初始條件:

(16)

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

(17)

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

Q=Q0δ(x-υt),

(18)

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

3 拉普拉斯域中問題求解

引入拉普拉斯變換公式:

(19)

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

利用上述拉普拉斯變換公式,方程(13)～方程(15)可變?yōu)?

(20)

(21)

(22)

對邊界條件進(jìn)行拉普拉斯變換得:

(23)

(24)

方程(24)的通解為:

(25)

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

(26)

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

k4-m1k2+m2=0,

(27)

可表達(dá)為:

(28)

(29)

方程(29)的通解為:

(30)

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

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

(31)

根據(jù)邊界條件方程(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)

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

(36)

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

(37)

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

(38)

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

(39)

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

(40)

4 數(shù)值反變換

(41)

其中:Re為實部;i為虛數(shù)單位。為了能加速收斂,結(jié)合大量已有試驗和現(xiàn)有計算結(jié)果,β應(yīng)滿足βt≈4.7[44]。

5 數(shù)值結(jié)果及分析

為了進(jìn)行數(shù)值計算,除了引入方程(40)進(jìn)行數(shù)值拉普拉斯反變換,得到時域內(nèi)桿的無量綱位移、應(yīng)力和溫度,還需要引入材料特性相關(guān)的參數(shù):λ=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。

數(shù)值計算得到了桿中無量綱位移、應(yīng)力和溫度的分布規(guī)律。計算中,研究了分?jǐn)?shù)階參數(shù)α、2個不同的熱松弛時間因子τ0和τ1、溫度相關(guān)性參數(shù)?以及移動熱源傳播速度υ對無量綱位移、應(yīng)力和溫度分布規(guī)律的影響效應(yīng),數(shù)值計算分析分別考慮4種不同的情形,每種情況分?jǐn)?shù)階參數(shù)均為α=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。所得結(jié)果如圖1～圖4所示。

(a) 分?jǐn)?shù)階參數(shù)下τ0變化對無量綱位移的影響

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

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

(a) 不同τ1和α?xí)r無量綱位移分布情況

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

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

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

(a) 分?jǐn)?shù)階參數(shù)下不同速度對無量綱位移的影響

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

(a) 無量綱位移與文獻(xiàn)[34]結(jié)果的比較

6 結(jié)論

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

(1)分?jǐn)?shù)階參數(shù)對各物理量均有一定的影響,但對無量綱位移和應(yīng)力的影響更明顯,分?jǐn)?shù)階參數(shù)會對無量綱量產(chǎn)生一定的滯后影響。

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

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