分?jǐn)?shù)階熱彈性理論下中空黏彈性圓柱熱彈耦合動(dòng)力響應(yīng)

郭 穎，熊春寶，虞跨海，梁 斌

(1.河南科技大學(xué) 土木工程學(xué)院，河南 洛陽(yáng) 471023；2.天津大學(xué) 建筑工程學(xué)院，天津 300072)

地下隧洞在軌道交通、石油及天然氣開(kāi)采、地?zé)豳Y源開(kāi)發(fā)、熱能儲(chǔ)存、核廢料處置等能源運(yùn)輸以及軍事等工程領(lǐng)域應(yīng)用廣泛，特別是當(dāng)兩個(gè)隧洞很近時(shí)，相互作用對(duì)間距具有敏感性，其耦合動(dòng)力響應(yīng)問(wèn)題直接影響其變形和位移，甚至可能造成失穩(wěn)。Biot[1]首次提出了耦合的熱彈性理論，該理論依舊認(rèn)為熱在介質(zhì)中的傳輸速度是無(wú)限大的，這與事實(shí)不符。為了消除這一現(xiàn)象學(xué)者們紛紛提出了可解釋熱的波動(dòng)效應(yīng)的熱彈性理論以及相應(yīng)的數(shù)學(xué)模型，常用的有：Lord-Shulman (L-S)[2]廣義熱彈性理論、Green-Lindsay (G-L)[3]廣義熱彈性理論以及Green-Naghdi (G-N)[4-6]廣義熱彈性理論。但上述理論不能很好的解釋有些黏彈性和反常擴(kuò)散問(wèn)題。

Abel首次將分?jǐn)?shù)階微積分引入到求解等時(shí)曲線問(wèn)題的積分方程當(dāng)中，并得到了很好的結(jié)果。直到1982年，Mandelbrot[7]發(fā)現(xiàn)并指出了自然界和工程領(lǐng)域中存在著大量的分?jǐn)?shù)維的事實(shí)，并提出該分?jǐn)?shù)維與原有的整數(shù)維之間存在自相似的現(xiàn)象，更好地描述了熱的波動(dòng)效應(yīng)中不好解釋的問(wèn)題。常用的分?jǐn)?shù)階廣義熱彈性理論有：Sherief等[8]型分?jǐn)?shù)階廣義熱傳導(dǎo)理論，Youssef[9]型分?jǐn)?shù)階廣義熱傳導(dǎo)理論，Ezzat等[10]型分?jǐn)?shù)階雙相滯后廣義熱傳導(dǎo)理論，Ezzat等[11]型分?jǐn)?shù)階三相滯后廣義熱傳導(dǎo)理論。聞敏杰等[12-13]基于分?jǐn)?shù)階熱彈性理論，研究了熱/力源作用下無(wú)限大土體中圓形隧洞的熱彈性耦合瞬態(tài)動(dòng)力響應(yīng)以及隧洞內(nèi)襯砌與土體的相互作用。分?jǐn)?shù)階微積分的引入對(duì)于黏彈性介質(zhì)的動(dòng)力特性有了更明確的解釋。

彈性介質(zhì)和黏彈性介質(zhì)往往被認(rèn)為均是可以滿足連續(xù)介質(zhì)的基本假設(shè)的?，F(xiàn)有的廣義和分?jǐn)?shù)階廣義熱彈性理論學(xué)者做了大量研究[14-17]，發(fā)現(xiàn)其研究結(jié)果并不能很好的描述介質(zhì)的松弛和蠕變等黏彈性特性。采用狀態(tài)空間法，Ezzat等[18]推導(dǎo)了二維黏彈性材料的熱黏彈基本方程，并借助Laplace-Fourier雙重反變換的方法對(duì)問(wèn)題進(jìn)行了求解，在已有方程的基礎(chǔ)上Ezzat等[19]又結(jié)合Kirchhoff理論研究了熱導(dǎo)率變換對(duì)無(wú)限長(zhǎng)粘彈性中空?qǐng)A柱的影響。Kar等[20]在廣義熱粘彈理論下，研究了均質(zhì)各向同性黏彈性球殼邊界受到溫度荷載作用的動(dòng)力響應(yīng)問(wèn)題，對(duì)比和分析了兩種不同理論中考慮黏彈性與不考慮黏彈性時(shí)各物理量之間的差異。Othman等[21]采用Laplace正、反變換法描述了雙溫度理論下邊界受到非高斯激光脈沖的作用的均質(zhì)各向同性Kelvin-Voigt黏彈半無(wú)限大固體的動(dòng)力響應(yīng)問(wèn)題。這些研究大多只對(duì)介質(zhì)的外表面或某個(gè)接觸面進(jìn)行了約束，未能考慮兩側(cè)均被約束時(shí)的動(dòng)力響應(yīng)問(wèn)題。

雖然分?jǐn)?shù)階廣義熱彈性理論已取得了較大發(fā)展，但將該理論應(yīng)用于研究外表面作用有熱沖擊的中空黏彈性柱體熱彈性分析的內(nèi)容鮮有報(bào)道。本文基于Ezzat型分?jǐn)?shù)階雙相滯后廣義熱彈性理論，引入Kelvin-Voigt黏彈性模型從而建立了黏彈性中空?qǐng)A柱熱彈耦合動(dòng)力模型，考慮軸對(duì)稱(chēng)各向同性均質(zhì)黏彈性中空柱外表面作用有熱沖擊的情況，利用Laplace變換以及特征值法，給出各物理量解的形式，并通過(guò)Laplace數(shù)值方便換得到各物理量的數(shù)值解。分析了黏彈性中空?qǐng)A柱熱彈耦合問(wèn)題。分析了荷載作用時(shí)間、分?jǐn)?shù)階系數(shù)和黏彈性松弛時(shí)間因子對(duì)中空?qǐng)A柱中無(wú)量綱位移、溫度、徑向應(yīng)力和環(huán)向應(yīng)力的影響。

1 熱黏彈模型建立及問(wèn)題描述

本文基于Ezzat型分?jǐn)?shù)階雙相滯后廣義熱彈性理論，引入Kelvin-Voigt黏彈性模型建立了均質(zhì)、各向同性黏彈性中空?qǐng)A柱熱彈耦合動(dòng)力模型。t=0時(shí)刻在圓柱外表面施加熱沖擊，其外邊界自由，內(nèi)邊界絕熱且固定(如圖1所示)。其不計(jì)體力和內(nèi)熱源的基本控制方程

(1)

(2)

(3)

(4)

圖1 黏彈性中空?qǐng)A柱熱彈耦合問(wèn)題模型示意圖Fig.1 Schematic diagram of the coupled thermoelastic problem of viscoelastic hollow cylindrical

Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)定義如下

(5)

根據(jù)研究對(duì)象的幾何結(jié)構(gòu)對(duì)稱(chēng)性，位移分量具有以下形式

ur=u(r,t),uφ=0,uz=0

(6)

基于式(6)，應(yīng)變分量為

(7)

結(jié)合柱坐標(biāo)，體積膨脹率e為

(8)

由式(1)～式(5)可得

(9)

(10)

(11)

(12)

(13)

引入如下無(wú)量綱量[22]

(14)

式(9)～式(13)采用式(14)無(wú)量綱化處理后仍使用略去各物理量右上方撇號(hào)的形式表示，則可進(jìn)一步寫(xiě)為

(15)

(16)

(17)

(18)

(19)

2 問(wèn)題求解

2.1 拉普拉斯域方程求解

將Riemann-Liouville拉普拉斯變換式[23]

(20)

將式(20)應(yīng)用于式(15)～式(19)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

消去式(26)和式(27)中的σ，得到滿足e的微分方程

(31)

同理，滿足θ的微分方程

(32)

式(31)和式(32)可被分解成如下形式

(33)

式(33)的根可展開(kāi)成

(34)

(35)

(36)

(37)

(38)

根據(jù)式(8)和式(24)，有

(39)

根據(jù)應(yīng)變分量式(7)和式(39)，可得

(40)

將式(34)、式(38)～式(40)代入式(28)～式(30)中，則有徑向應(yīng)力、環(huán)向應(yīng)力和軸向應(yīng)力的表達(dá)式

(41)

(42)

(43)

為了確定待定系數(shù)A1和Di，須引入問(wèn)題的初始條件和邊界條件，初始條件和邊界條件如下：

初始條件

(44)

經(jīng)過(guò)拉普拉斯變換后的邊界條件

(45)

(46)

根據(jù)初邊值條件，結(jié)合式(38)、式(39)和式(41)，有如下的方程組

(47)

(48)

(49)

(50)

求解Ai和Di的過(guò)程較為繁雜，此處不再贅述。

2.2 拉普拉斯數(shù)值反變換

為得到時(shí)間域中溫度、位移、徑向應(yīng)力和環(huán)向應(yīng)力的解答及分布規(guī)律，需對(duì)拉氏域中得到的各物理量的解進(jìn)行拉普拉斯反變換。由于拉氏域中各物理量解的復(fù)雜性，要對(duì)其進(jìn)行解析反變換幾乎是不可能的，只能采取數(shù)值反變換的方法。為此，引入如下拉普拉斯數(shù)值反變換方法[24]

(51)

在進(jìn)行拉氏數(shù)值反變換時(shí)，應(yīng)用了下面關(guān)于貝塞爾函數(shù)的漸近展開(kāi)式

(52)

式中，μ1=4v2。

3 數(shù)值結(jié)果及討論

本文以黏彈性中空?qǐng)A柱為例，在圓柱外表面施加熱沖擊。結(jié)合邊界條件最終得到了作用時(shí)間、分?jǐn)?shù)階系數(shù)和黏彈性松弛時(shí)間因子變化等對(duì)該柱體中無(wú)量綱位移、溫度、將向應(yīng)力和環(huán)向應(yīng)力的影響。本算例中所需參數(shù)同參考文獻(xiàn)[25]

其他參數(shù)如下：

θ0=1,r1=3,r2=1

當(dāng)α=1，α0=0，α1=0時(shí)，分?jǐn)?shù)階黏彈性問(wèn)題將退化為對(duì)應(yīng)的廣義熱彈性理論。此時(shí)，式(37)和式(38)與文獻(xiàn)[26]研究圓柱體熱彈性問(wèn)題所求得的解形式相同，這驗(yàn)證了本文根據(jù)特征值法在Laplace變換域內(nèi)所求得結(jié)果的正確性。

表1 考慮熱彈耦合效應(yīng)的黏彈性中空計(jì)算參數(shù)Tab 1 Calculation parameters of saturated porous viscoelastic foundation considering the coupled thermoelastic effect

此外，本文還將MATLAB程序中取α=1，α0=0，α1=0即可退化為整數(shù)階廣義熱彈性問(wèn)題，并結(jié)合朱海陶等研究中的邊界條件對(duì)問(wèn)題進(jìn)行求解，并將無(wú)量綱位移和溫度與朱海陶等研究的相應(yīng)結(jié)果進(jìn)行對(duì)比。結(jié)果表明(如圖2所示)，當(dāng)α=1，α0=0，α1=0時(shí)，結(jié)果吻合較好，有微小差異是由于數(shù)值計(jì)算方法的不同引起，其中最大偏差為5.88%，平均偏差為3.26%。從而進(jìn)一步驗(yàn)證了本文計(jì)算結(jié)果的正確性和合理性。

圖2 無(wú)量綱量與朱海陶等研究結(jié)果的比較Fig.2 Comparison between the non-dimension physical variables results in this paper and those in zhuhaitao,et al

為分析溫度、位移、徑向應(yīng)力和環(huán)向應(yīng)力隨時(shí)間t、分?jǐn)?shù)階系數(shù)α和黏彈性松弛時(shí)間因子α0和α1的變化規(guī)律，計(jì)算時(shí)考慮了兩個(gè)時(shí)間，即t=0.2，t=0.25，三個(gè)不同分?jǐn)?shù)階系數(shù)和黏彈性松弛時(shí)間因子α，α0和α1，分別是α=0.1，α=0.5，α=1；α0=0，α0=0.06，α0=0.18；和α1=0，α1=0.09，α1=0.27時(shí)幾種情況的不同組合。得到了無(wú)量綱溫度θ、位移u、徑向應(yīng)力σrr和環(huán)向應(yīng)力σφφ的分布規(guī)律，如圖3～圖5所示。當(dāng)α0=0，α1=0時(shí)，黏彈性模型就可以退化成為彈性模型；當(dāng)α=1時(shí)，分?jǐn)?shù)階模型可以退化成為整數(shù)階模型。

圖3 α0變化對(duì)各無(wú)量綱量的影響Fig.3 Influence of different α0 on non-dimension physical variables

圖3描述了熱沖擊時(shí)，黏彈性松弛時(shí)間因子α0變化對(duì)無(wú)量綱溫度θ、位移u、徑向應(yīng)力σrr和環(huán)向應(yīng)力σφφ的影響。文中正值表示介質(zhì)處于膨脹狀態(tài)，負(fù)值表示介質(zhì)處于壓縮狀態(tài)。圖3(a)是中空?qǐng)A柱中的位移分布圖，由于柱體受到熱沖擊，外表面處發(fā)生熱膨脹變形，柱體內(nèi)的位移從膨脹狀態(tài)到壓縮狀態(tài)動(dòng)態(tài)變化，到內(nèi)表面時(shí)恢復(fù)到0。當(dāng)α0和α1都為0時(shí)，曲線在壓縮區(qū)域的峰值分別在r=2.6和r=2.7處出現(xiàn)。同一作用時(shí)間下，隨著α0增大，曲線峰值向內(nèi)移動(dòng)，這主要是由于黏彈性介質(zhì)的遲滯作用引起的，黏彈性松弛時(shí)間因子α0對(duì)位移的影響主要在柱體外表面到柱體內(nèi)外表面之間處最為明顯，當(dāng)r=1.75時(shí)，幾條曲線基本重合并衰減為零，擾動(dòng)基本消失。此外，對(duì)于彈性介質(zhì)而言，作用時(shí)間增大使得曲線壓縮區(qū)域的峰值向內(nèi)移動(dòng)，但數(shù)值差異不大，而對(duì)于黏彈性介質(zhì)而言，作用時(shí)間僅對(duì)曲線峰值的幅值大小有微弱的影響，而對(duì)峰值的位置影響不大。

圖3(b)是中空?qǐng)A柱的溫度變化圖，由于柱體外側(cè)受到熱荷載作用，熱擾動(dòng)區(qū)域從柱體外表面向內(nèi)部遷移，且溫度從1～0依次減小。黏彈性松弛時(shí)間因子α0對(duì)溫度的變化影響不大，隨著作用時(shí)間的增大，溫度衰減速度變快，擾動(dòng)深度向柱體內(nèi)部延伸。

圖3(c)和圖3(d)分別是徑向應(yīng)力和環(huán)向應(yīng)力分布圖。圖3(c)中徑向應(yīng)力在柱體內(nèi)外表面處均為零，并隨著徑向深度增大，應(yīng)力先表現(xiàn)處壓縮狀態(tài)，直到r>2.6后慢慢進(jìn)入膨脹狀態(tài)。隨著黏彈性松弛時(shí)間因子α0增大，曲線在壓縮區(qū)域的峰值向柱體內(nèi)徑方向移動(dòng)，產(chǎn)生明顯的滯后現(xiàn)象。不考慮黏彈性時(shí)，曲線在壓縮區(qū)域峰值基本在r=2.7時(shí)出現(xiàn)，而在膨脹區(qū)域的峰值出現(xiàn)的位置則有明顯的不同。

圖3(d)中在熱沖擊作用下，環(huán)向產(chǎn)生壓應(yīng)力，同一作用時(shí)間下，黏彈性松弛時(shí)間因子α0增大使得曲線峰谷值均向柱體內(nèi)徑方向移動(dòng)，谷值處尤為明顯?？紤]黏彈性松弛時(shí)間因子作用時(shí)，不同時(shí)間的曲線谷值位置差異很大，但幅值比較接近，不考慮黏彈性松弛時(shí)間因子的情況則反之。而對(duì)于曲線峰值而言，考慮黏彈性時(shí)間因子的曲線峰值大小和位置都相對(duì)接近，不考慮黏彈性松弛時(shí)間因子的則反之。

在圖3(c)和圖3(d)中，當(dāng)作用時(shí)間t=0.2時(shí)，徑向應(yīng)力和環(huán)向應(yīng)力的三條曲線在1≤r<1.5的區(qū)域內(nèi)基本都是零，隨著作用時(shí)間增大，這個(gè)現(xiàn)象逐漸消失，說(shuō)明其擾動(dòng)向柱體內(nèi)部延伸。

綜上所述，作用時(shí)間對(duì)所有物理量均有明顯影響，而黏彈性松弛時(shí)間因子則對(duì)無(wú)量綱溫度外的所有物理量均有明顯影響。

圖4 α1變化對(duì)各無(wú)量綱量的影響Fig.4 Influence of different α1 on non-dimension physical variables

圖4描述了熱沖擊時(shí)，黏彈性松弛時(shí)間因子α1變化對(duì)無(wú)量綱溫度θ、位移u、徑向應(yīng)力σrr和環(huán)向應(yīng)力σφφ的影響。圖4(a)是中空?qǐng)A柱中的位移分布圖，曲線分布趨勢(shì)與圖3(a)黏彈性松弛時(shí)間因子α0變化時(shí)類(lèi)似，但沒(méi)有α0變化時(shí)影響明顯。不同的是，α1變化不僅對(duì)曲線峰值出現(xiàn)的位置有一定的影響，而且還影響了其幅值的大小。圖4(b)與圖3(b)類(lèi)似，黏彈性松弛時(shí)間因子α1對(duì)溫度的變化影響不大。圖4(c)和圖4(d)分別是徑向應(yīng)力和環(huán)向應(yīng)力分布圖。黏彈性松弛時(shí)間因子α1對(duì)徑向應(yīng)力的影響明顯，與圖3(c)不同的是，隨著α1增大，考慮α1的曲線谷值不僅向著圓柱內(nèi)表面方向移動(dòng)，數(shù)值也在逐漸減小。此外，考慮了黏彈性松弛時(shí)間因子作用的曲線在膨脹區(qū)域的數(shù)值非常小，只有微弱的擾動(dòng)。圖4(d)中環(huán)向同樣僅產(chǎn)生壓應(yīng)力，在不考慮黏彈性松弛時(shí)間因子時(shí)，曲線谷值明顯靠近柱體外表面，并且更大一些?？紤]黏彈性松弛時(shí)間因子作用時(shí)，隨著α1增大，曲線谷值大小和位置都非常接近。當(dāng)α1=0.27時(shí)，曲線的峰谷值基本不太明顯，這主要是黏滯性使得曲線的增大和衰減速度明顯減慢造成的。

綜上，由考慮兩種不同黏彈性松弛時(shí)間因子的兩組圖可知，黏彈性松弛時(shí)間因子α0和α1對(duì)無(wú)量綱溫度外的各物理量均有明顯的滯后作用，這個(gè)現(xiàn)象在曲線的峰值和谷值處最為明顯。若僅需計(jì)算無(wú)量綱溫度的分布情況可以不考慮黏彈性松弛時(shí)間因子α0和α1，從而簡(jiǎn)化計(jì)算。但是，往往熱沖擊作用會(huì)產(chǎn)生變形、位移和應(yīng)力，僅考慮溫度變化是不夠的，因此，研究黏彈性松弛時(shí)間因子α0和α1對(duì)各物理量的研究是不可忽略的。

圖5描述了熱沖擊時(shí)，分?jǐn)?shù)階系數(shù)變化對(duì)無(wú)量綱溫度θ、位移u、徑向應(yīng)力σrr和環(huán)向應(yīng)力σφφ的影響。分?jǐn)?shù)階系數(shù)對(duì)無(wú)量綱位移的影響最小，圖5(a)中，分?jǐn)?shù)階系數(shù)對(duì)無(wú)量綱位移的影響主要體現(xiàn)在曲線峰值及其附近區(qū)域。當(dāng)r<1.5后，熱沖擊所產(chǎn)生的變形和位移都很小了，這是分?jǐn)?shù)階系數(shù)變化對(duì)其影響也是一個(gè)很小的值，即可忽略其影響。隨著分?jǐn)?shù)階系數(shù)增大，曲線峰值逐漸增大。圖5(b)是圓柱中的溫度分布圖，由于熱沖擊作用在圓柱外表面，熱量向內(nèi)傳遞并逐漸減小，分?jǐn)?shù)階系數(shù)使得溫度曲線衰減明顯變慢，這與分?jǐn)?shù)階微分算子的記憶性有一定的關(guān)系。隨著作用時(shí)間的增大，不同分?jǐn)?shù)階系數(shù)曲線之間的差異明顯增大，這個(gè)現(xiàn)象在1.75

圖5 分?jǐn)?shù)階系數(shù)變化對(duì)各無(wú)量綱變量的影響Fig.5 Influence of fractional order coefficient variation on non-dimension physical variables

4 結(jié) 論

本文基于Ezzat型分?jǐn)?shù)階雙相滯后廣義熱彈性理論，根據(jù)特征值法以及拉普拉斯數(shù)值反變換解法對(duì)中空均質(zhì)各向同性黏彈性柱外表面作用有熱沖擊的問(wèn)題進(jìn)行了研究，得到了在相應(yīng)熱沖擊下無(wú)限長(zhǎng)中空柱體無(wú)量綱位移、溫度、徑向應(yīng)力和環(huán)向應(yīng)力的分布規(guī)律。通過(guò)對(duì)比分析揭示了分?jǐn)?shù)階系數(shù)、作用時(shí)間和黏彈性松弛時(shí)間因子對(duì)所考慮各物理量的影響規(guī)律，結(jié)論如下：

(1) 當(dāng)r<1.75后熱沖擊對(duì)無(wú)量綱位移的擾動(dòng)基本消失，不同作用時(shí)間對(duì)擾動(dòng)深度的影響不大。而對(duì)于無(wú)量綱溫度、徑向應(yīng)力和環(huán)向應(yīng)力，隨著作用時(shí)間增大，擾動(dòng)區(qū)域向柱體內(nèi)部延伸。

(2) 黏彈性松弛時(shí)間因子α0和α1對(duì)無(wú)量綱溫度外的所有物理量具有明顯的影響，主要在曲線的峰值和谷值處有明顯的滯后現(xiàn)象，對(duì)徑向應(yīng)力和環(huán)向應(yīng)力的影響更為明顯，主要體現(xiàn)在1.5

(3) 分?jǐn)?shù)階系數(shù)對(duì)所有物理量均有明顯影響，主要體現(xiàn)在幅值最大處及周?chē)鷧^(qū)域影響最為明顯，當(dāng)r<1.5后，分?jǐn)?shù)階系數(shù)對(duì)各物理量的影響基本消失。考慮分?jǐn)?shù)階系數(shù)后，曲線更平滑，除環(huán)向應(yīng)力外的各曲線峰值同樣有明顯的滯后。此外，本文的研究結(jié)果，可對(duì)分析具體熱彈性行為α?xí)r的取值提供一定的理論參考。