地下隧洞在城市轨道交通、石油和天然气运输、电力管线、军事等工程领域得到广泛应用。近年来,由于恐怖分子平凡活动,这些结构内部常受到爆炸、冲击等动力荷载破坏作用,严重影响了地下结构的稳定性[1,2]。爆炸冲击时瞬间既引起较大的冲击波,还会产生很高的温度导致隧洞的变形破坏。因此,地下隧洞在热、力冲击作用下的动力响应问题备受许多学者的关注。基于热弹性理论,Abd-alla等[3]研究了磁场和温度场对无限圆柱形孔洞动力响应的影响。基于广义热弹性理论,He等[4]考虑材料的热物性参数随温度的变化,建立了广义非线性热弹性模型,采用有限法研究了半无限体在热冲击作用下的瞬态响应。孙飞等[5]研究了半空间非均匀各向同性功能梯度体的热弹性问题,利用状态空间法建立了热力耦合控制方程,通过特征值分析法和Laplace逆变换得到了各物理场的解析解。Abbas[6]引入两个松弛时间,采用有限元法研究了非均质圆柱形孔洞的热冲击问题,分析了三种理论模型下温度梯度、位移和应力随时间的变化规律。基于广义热扩散理论,Aouadi[7]研究了具有球形孔洞的无限弹性体热-力-化学耦合动力响应,证明了与热弹性模型(C-T模型)的位移、应力幅值存在较大差异。Sherief等[8]研究了热冲击荷载作用下无限圆柱形孔洞的动力响应,分析了不同时刻位移和应力随半径的变化规律。郑荣跃等[9]研究了无限弹性介质中圆形隧洞内壁受时间变化热源作用时的动力响应问题,得到了热、力和化学耦合条件下温度、应力、位移和化学势响应的分布规律。

但是,在急速冲击作用下,经典热弹性理论和广义热弹性理论难以有效描述其热弹性行为。Youssef[10,11]将Riemann-Liouville分数阶积分算子引入到热传导方程中,建立了分数阶广义热弹性理论,并成功应用于分析半无限大体热冲击问题[12,13,14],该理论通过改变延迟效应来影响热冲击的作用效果,从而更好地描述弹性波及热波传播的特性。本文基于分数阶广义热弹性理论,研究了热、力冲击作用下无限大土体中圆形隧洞的瞬态动力响应。利用Laplace变换及其逆变换得到了圆形隧道周围土体的无量纲温度、径向位移、径向应力和环向应力等的解析表达式。在此基础上,分析了分数阶参数对热/力源条件下土体的温度、径向位移、径向应力和环向应力的影响。

1、土体的控制方程及求解过程

1.1 热源作用下的土体控制方程

考虑一含有半径为R的圆形隧洞的无限土体,内壁表面处作用一随时间t变化的径向均布爆炸荷载为f(t)=4q0(e-t/2-e-t)sin(t),受到热冲击作用为T(t)=θ0e-0.5t,q0为冲击力荷载幅值,θ0为热荷载幅值。

采用极坐标系,坐标原点在圆形隧洞的中心处。轴对称作用下,土体的运动方程为[15]

(λ+2μ)∂e∂r−β∂θ∂r=ρ∂2ur∂t2         (1)(λ+2μ)∂e∂r-β∂θ∂r=ρ∂2ur∂t2         (1)

式中:λ、μ为Lame常数;e=1r∂∂r(rur)e=1r∂∂r(rur)为体应变;β=(3λ+2μ)αt,αtβ=(3λ+2μ)αt,αt为线性热膨胀系数;ρ为土体的密度;ur为土体的径向位移;θ=T-T0为温度增量,T,T0分别为绝对温度和初始参考温度,满足|(T−T0)/T0|<<1。|(Τ-Τ0)/Τ0|<<1。

Sherief通过引入Caputo型的分数阶导数,给出了分数阶热传导方程表达式[16]

κ∇2θ=(1+τ0∂α∂tα)(ρcE∂θ∂t+T0β∂e∂t)         (2)κ∇2θ=(1+τ0∂α∂tα)(ρcE∂θ∂t+Τ0β∂e∂t)         (2)

式中,κ为热传导系数;∇2=∂2∂r2+1r∂∂r∇2=∂2∂r2+1r∂∂r为拉普拉斯算子;τ0为热松弛时间;cE为常应变下比热;T0为初始温度。α为分数阶参数,表示为延迟时间影响因子。因此,通过改变延迟效应来影响热冲击的作用效果,从而更好地描述弹性波及热波传播的特性。当α=1时,上述热传导方程式(2)退化为Lord-Shulman(LS)广义热弹性理论。

本构方程为

{σr=2μ∂ur∂r+λe−βθσθ=2μurr+λe−βθ         (3){σr=2μ∂ur∂r+λe-βθσθ=2μurr+λe-βθ         (3)

式中,σr和σθ为径向应力和环向应力。

1.2 方程的求解

引入如下无量纲量

⎧⎩⎨⎪⎪r∗=c0η0r,u∗r=c0η0ur,σ∗r=σr/(λ+2μ)t∗=c20η0t,τ∗0=c20τ0η0,c20=(λ+2μ)/ρη0=ρcE/κ,θ∗=βθ/(λ+2μ)         (4){r*=c0η0r,ur*=c0η0ur,σr*=σr/(λ+2μ)t*=c02η0t,τ0*=c02τ0η0,c02=(λ+2μ)/ρη0=ρcE/κ,θ*=βθ/(λ+2μ)         (4)

为简便起见,去掉各量右上方的“*”,将式(4)代入式(1)～式(3),可得:

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∂2ur∂t2=∂e∂r−∂θ∂r∇2θ=(1+τ0∂α∂tα)(∂θ∂t+ε∂e∂t)σr=e−2β21urr−θσθ=(1−2β21)e+2β21urr−θ         (5){∂2ur∂t2=∂e∂r-∂θ∂r∇2θ=(1+τ0∂α∂tα)(∂θ∂t+ε∂e∂t)σr=e-2β12urr-θσθ=(1-2β12)e+2β12urr-θ         (5)

式中,ε=β2T0ρ2cEc20ε=β2Τ0ρ2cEc02为热弹性耦合系数。

将Riemann-Liouville拉普拉斯变换式

L[Inf(t)]=1snL[f(t)],n>0         (6)L[Ιnf(t)]=1snL[f(t)],n>0         (6)

应用于式(5),可得:

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪s2u¯r=∂e−∂r−∂θ−∂r∇2θ−=(s+τ0sα)θ−+ε(s+τ0sα)e−σ−r=e−−2β21u¯rr−θ−σ−θ=(1−2β21)e−+2β21u¯rr−θ−         (7){s2u¯r=∂e-∂r-∂θ-∂r∇2θ-=(s+τ0sα)θ-+ε(s+τ0sα)e-σ-r=e--2β12u¯rr-θ-σ-θ=(1-2β12)e-+2β12u¯rr-θ-         (7)

对式(7)第一式运用微分算子∇=∂∂r+1r∇=∂∂r+1r,可得

(∇2−s2)e−=∇2θ−         (8)(∇2-s2)e-=∇2θ-         (8)

将式(8)代入式(7)第二式,得到:

∇4e−−[(s+τ0sα)(ε+1)+s2]∇2e−+s2(s+τ0sα)e−=0         (9)∇4e--[(s+τ0sα)(ε+1)+s2]∇2e-+s2(s+τ0sα)e-=0         (9)

同理,可得:

∇4θ−−[(s+τ0sα)(ε+1)+s2]∇2θ−+s2(s+τ0sα)θ−=0         (10)∇4θ--[(s+τ0sα)(ε+1)+s2]∇2θ-+s2(s+τ0sα)θ-=0         (10)

于是,式(9)和式(10)分别可写成

(∇4−d21∇2+d2)e−=0         (11)(∇4−d21∇2+d2)θ−=0         (12)(∇4-d12∇2+d2)e-=0         (11)(∇4-d12∇2+d2)θ-=0         (12)

其中,

{d1=(s+τ0sα)(ε+1)+s2d2=s2(s+τ0sα)         (13){d1=(s+τ0sα)(ε+1)+s2d2=s2(s+τ0sα)         (13)

进而将式(11)和式(12)分解为

(∇2−k21)(∇2−k22)e−=0         (14)(∇2−k21)(∇2−k22)θ−=0         (15)(∇2-k12)(∇2-k22)e-=0         (14)(∇2-k12)(∇2-k22)θ-=0         (15)

其中,k21,2=d1±d21−4d2√2k1,22=d1±d12-4d22

根据算子分解理论有

e−=e−1+e−2         (16)e-=e-1+e-2         (16)

其中,e−1,e−2Ze-1,e-2Ζ满足下式

(∇2−k21)e−1=0         (17)(∇2−k22)e−2=0         (18)(∇2-k12)e-1=0         (17)(∇2-k22)e-2=0         (18)

式(17)和式(18)可解得

e−1=A1K0(k1r)+B1I0(k1r)         (19)e−2=A2K0(k2r)+B2I0(k2r)         (20)e-1=A1Κ0(k1r)+B1Ι0(k1r)         (19)e-2=A2Κ0(k2r)+B2Ι0(k2r)         (20)

式中:Kn(·)和In(·)分别为n阶第二类和第一类Bessel函数;A1,A2,B1,B2为待定系数。

由边界条件,当r→∞Z时,应力和位移都趋向于零,结合虚宗量Bessel函数性质,有:

{Kn(kr)∼(π2kr)12e−krIn(kr)∼(2πkr)−12ekr,r→∞         (21){Κn(kr)∼(π2kr)12e-krΙn(kr)∼(2πkr)-12ekr,r→∞         (21)

因此,B1=B2=0

于是,可得:

e−=A1K0(k1r)+A2K0(k2r)         (22)θ−=D1K0(k1r)+D2K0(k2r)         (23)e-=A1Κ0(k1r)+A2Κ0(k2r)         (22)θ-=D1Κ0(k1r)+D2Κ0(k2r)         (23)

式中,D1,D2为待定系数。

利用Bessel函数的基本性质,由式(22)可得径向位移的表达式为

u¯r=−A1k1K1(k1r)−A2k2K1(k2r)         (24)u¯r=-A1k1Κ1(k1r)-A2k2Κ1(k2r)         (24)

因此,由式(7)第三式和第四式,可得径向应力与环向应力为

σ−r=2rβ21∑i=12AikiK1(kir)+∑i=12(Ai−Di)K0(kir)         (25)σ−θ=∑i=12[(1−2β21)Ai−Di]K0(kir)−2rβ21∑i=12AikiK1(kir)         (26)σ-r=2rβ12∑i=12AikiΚ1(kir)+∑i=12(Ai-Di)Κ0(kir)         (25)σ-θ=∑i=12[(1-2β12)Ai-Di]Κ0(kir)-2rβ12∑i=12AikiΚ1(kir)         (26)

将式(23)和式(24)代入式(7)第一式,可得如下关系

Di=k2i−s2k2iAi         (27)Di=ki2-s2ki2Ai         (27)

将式(27)代入式(25)和式(26),可得:

σ−r=∑i=12[2rβ21K1(kir)ki+s2K0(kir)k2i]Ai         (28)σ−θ=∑i=12{[(1−2β21)−k2i−s2k2i]K0(kir)−2rβ21K1(kir)ki}Ai         (29)σ-r=∑i=12[2rβ12Κ1(kir)ki+s2Κ0(kir)ki2]Ai         (28)σ-θ=∑i=12{[(1-2β12)-ki2-s2ki2]Κ0(kir)-2rβ12Κ1(kir)ki}Ai         (29)

2、边界条件及求解

由内壁表面处作用一随时间t变化的径向均布爆炸荷载且受到热冲击作用,可知:

σ−r=−f−(s),θ−=T−(s) r=R         (30)σ-r=-f-(s),θ-=Τ-(s) r=R         (30)

其中,式f−(s),T−(s)f-(s),Τ-(s)为f(t),T(t)的拉普拉斯变换表达式

⎧⎩⎨⎪⎪⎪⎪f−(s)=4q0(11+(0.5+s)2−11+(1+s)2)T−(s)=θ01s+0.5         (31){f-(s)=4q0(11+(0.5+s)2-11+(1+s)2)Τ-(s)=θ01s+0.5         (31)

将式(23)、(27)和式(28)代入式(30),可得:

⎧⎩⎨⎪⎪⎪⎪⎪⎪A2=−a1T−(s)+f−(s)b1a1b2−a2b1A1=−a2T−(s)+f−(s)b2a1b2−a2b1         (32){A2=-a1Τ-(s)+f-(s)b1a1b2-a2b1A1=-a2Τ-(s)+f-(s)b2a1b2-a2b1         (32)

其中,

a1=2Rβ21K1(k1R)k1+s2K0(k1R)k21         (33)a2=2Rβ21K1(k2R)k2+s2K0(k2R)k22         (34)b1=k21−s2k21K0(k1R)         (35)b2=k22−s2k22K0(k2R)         (36)a1=2Rβ12Κ1(k1R)k1+s2Κ0(k1R)k12         (33)a2=2Rβ12Κ1(k2R)k2+s2Κ0(k2R)k22         (34)b1=k12-s2k12Κ0(k1R)         (35)b2=k22-s2k22Κ0(k2R)         (36)

由此,可得热、力冲击作用下具有圆形隧洞的周围土体瞬态动力响应。

3、问题的退化解

3.1 与文献[17]对比

文献[17]研究了无限介质中圆形孔洞的瞬态响应,并未考虑热力耦合效应,即忽略温度的影响,则土体的运动方程和本构关系退化为

(λ+2μ)∂e∂r=ρ∂2ur∂t2         (37)σr=2μ∂ur∂r+λe−βθ         (38)σθ=2μurr+λe−βθ         (39)(λ+2μ)∂e∂r=ρ∂2ur∂t2         (37)σr=2μ∂ur∂r+λe-βθ         (38)σθ=2μurr+λe-βθ         (39)

3.2热弹性模型对比

当α=1时,上述热传导方程式(2)退化为Lord-Shulman(L-S)广义热弹性模型,则有:

κ∇2θ=(1+τ0∂∂t)(ρcE∂θ∂t+T0β∂e∂t)         (40)κ∇2θ=(1+τ0∂∂t)(ρcE∂θ∂t+Τ0β∂e∂t)         (40)

当热松弛时间τ0为零时,退化为经典热弹性模型,则有

κ∇2θ=ρcE∂θ∂t+T0β∂e∂t         (41)κ∇2θ=ρcE∂θ∂t+Τ0β∂e∂t         (41)

4、数值结果分析

为了得到时域中温度增量、径向位移、径向应力和环向应力的解答,直接对其进行Laplace逆变换较为困难。为此,本文采用如下Crump反演法计算[18]。

设函数F(s)为函数F(t)的Laplace变换,则Laplace逆变换的Crump反演算法为

F(t)≈eatT∗{12F(a)+∑k=1∞[Re[F(a+kπiT∗)]coskπtT∗−ImF(t)≈eatΤ*{12F(a)+∑k=1∞[Re[F(a+kπiΤ*)]coskπtΤ*-Ιm[F(a+kπiT∗)]sinkπtT∗]}         (42)[F(a+kπiΤ*)]sinkπtΤ*]}         (42)

若|F(t)|t/2。

在算例中,根据土体材料的特性,取文献[7]的部分参数,具体如下:

λ=7.76×107Pa,G=μ=3.86×107Pa,ρ=1800kg/m3,τ0=0.02s,T0=293K,αt=1.78×10-5℃-1,CE=2000m2/s2/℃,κ=3.8W。

4.1 计算与验证

如图1所示,本文的解与文献[17]的计算结果进行了对比分析。文献[17]未考虑温度与应力的相互作用,研究了力源作用下无限弹性介质中圆形孔洞的动力响应。本文忽略温度的影响,可退化为文献[17]的计算结果。可见,两者的结果几乎一致,验证了本文计算结果的正确性。

4.2 热源作用

图2～图5表示热源作用下分数阶参数α对径向位移ur、温度增量θ、径向应力σr和环向应力σθ的影响。分别取α=0.25和α=1时,分析圆形隧道内壁处(r=1)径向位移ur、温度增量θ、径向应力σr和环向应力σθ随无量纲时间t的变化规律。可见,随着分数阶参数α的增加,径向位移ur的峰值逐渐减小,但随着时间t的增加,衰减放慢(图2)。分数阶参数α对温度增量的影响显著,参数α越大,则峰值越大(图3)。图4和图5表明在圆形隧道内壁处(r=1)。分数阶参数α的变化值对径向应力σr无影响且都为零,而α=1对应的环向应力σθ值明显大于α=0.25对应的值。当r≈3时为临界点,径向应力的值有突变,r>3时,随着分数阶参数α的增加,径向应力σr逐渐减小。由此可见,热源作用下分数阶参数α对径向位移ur影响较小,而对温度增量θ、径向应力σr和环向应力σr影响较大。

图1与文献[17]计算结果对比

图2分数阶参数α对径向位移ur的影响

图3分数阶参数α对温度增量θ的影响

图4分数阶参数α对径向应力σr的影响

图5分数阶参数α对环向应力σθ的影响(热源)

4.3 力源作用

图6和图7表示力源荷载作用下剪切模量μ和分数阶参数α分别对径向位移ur的影响。由图6可见,随着剪切模量μ的增加,径向位移ur的峰值逐渐减小。这是由于剪切模量越大,土体越硬导致的。由图7可见,在力源作用下,分数阶参数α对响应幅值无任何影响。这是因为此时无热源作用,在沿半径未发生温差,从而未形成温度应力而产生变形。

图6剪切模量μ对径向位移ur的影响

图7分数阶参数α对径向位移ur的影响(力源)

4.4 热力耦合作用

图8～图11表示热力耦合作用下分数阶参数α对径向位移ur、温度增量θ、径向应力σr和环向应力σθ的影响。分别取α=0.25和α=1时,可见,在热力耦合作用下分数阶参数α对径向位移ur的影响较小。这是由于热力的叠加效应引起的。随着分数阶参数α的增加,温度增量θ逐渐增大。分数阶参数α对径向应力σr和环向应力σθ的影响较小。

图8热力耦合作用下分数阶参数α对径向位移ur的影响

图9热力共同作用下温度增量θ随半径r的衰减规律

图10热力共同作用下径向应力σr随半径r的衰减规律

图11热力共同作用下环向应力σθ随半径r的衰减规律

5、结论

本文基于分数阶广义热弹性理论,研究了热、力、热力共同冲击作用下无限大土体中圆形隧洞的热弹性耦合瞬态动力响应。利用Laplace变换及其逆变换,得到了土体的动力响应,得到了如下结论:

(1)热源作用下分数阶参数对径向位移的影响较小,对温度增量、径向应力和环向应力影响较大。而在热力共同作用下,分数阶参数对土体的温度增量影响显著,而对径向位移、径向应力和环向应力的影响较小。

(2)随着剪切模量的增加,位移的峰值逐渐减小。

(3)热源作用下分数阶参数对响应幅值的影响大于热力共同作用下分数阶参数对响应幅值的影响。力源作用下分数阶参数对径向位移无影响。

参考文献：

[1]张小勇,龚顺风.隧道内爆炸作用下衬砌结构损伤机理及抗爆性能研究[J].振动与冲击,2013,32(22):193-199.

[2]赵跃堂,冠伟晓,储程,等.管片衬砌结构在接触爆炸作用下的毁伤特性分析[J].振动与冲击,2018,37(11):95-100.

[5]孙飞,刘五祥,张伦伟.功能梯度半空间体的热弹性问题研究[J].力学季刊,2017,38(3):537-544.

[9]郑荣跃,刘干斌,唐国金.考虑扩散效应圆形隧洞热弹性耦合动力响应研究[J].国防科技大学学报,2008,30(3):27-31.

[17]杨桂通,张善元.弹性土动力学[M].北京:中国铁道出版社,1988.

闻敏杰,徐金明,熊厚仁.基于分数阶热弹性理论的圆形隧洞热弹性耦合动力响应[J].振动与冲击,2020,39(21):89-94+110.

基金:国家自然科学基金(41472254).