Analytical solutions for 1-D thermo-consolidation process of semi-infinite saturated soil using similarity transformation method
-
摘要: 针对半无限空间饱和土在表面外荷载及变温作用影响下的一维热固结过程,考虑两类传热及渗流边界条件,给出了温度及超静孔压耦合演化过程的基本方程。引入组合变量对耦合方程进行了解耦,应用相似变换法获得了组合变量在半整数次幂函数边界条件下的通解,并基于该通解构建了任意幂级数形式荷载及边界条件下温度与超静孔压的解析解答。通过文献进行了对比验证,随后应用本解答计算分析了半无限空间饱和土在表面正弦温度及热流边界作用下土中温度及超静孔压的耦合响应特性,结果表明:土体变形功项对传热方程贡献很小,对温度及超静孔压演化的影响基本可以忽略;固结方程中温度化率项系数为正值时,温度变化会引起相反趋势的超静孔压变化,而该系数为负值时,温度变化则会引起相同趋势的超静孔压变化。Abstract: For a 1-D thermo-consolidation process of semi-infinite saturated soil induced by external loads and changing temperature at the surface, basic equations are proposed to describe the coupled evolvement of temperature and excessive porewater pressure, and two types of heat transfer and seepage boundary conditions are considered. The combination variables are introduced to decouple the governing equations, and their general solutions under the boundary conditions in the form of half-integer power functions are developed by applying the similarity transformation method. After expanding the external loads and the boundary conditions into power series, the analytical solutions for the temperature and the excessive porewater pressure can be obtained using these general solutions directly. After their verification by certain solution in the literature, our solutions are then employed to calculate and analyze the coupled response characteristics of soil temperature and excessive porewater pressure for semi-infinite saturated soil under the effects of sinusoidal temperature and heat flux at the soil surface. The results show that, the term of soil deformation work contributes very little to the heat transfer equation, and its impact on the evolution of temperature and excessive porewater pressure can basically be ignored. When the coefficient for the term of temperature changing rate in the consolidation equation is positive, temperature changes will cause changes in the excessive porewater pressure with the opposite trend, and when the coefficient is negative, temperature changes will cause changes in the excessive porewater pressure with the same trend.
-
0. 引言
饱和岩土体非等温条件下的热-水-力耦合问题在地热能开发利用、核废料处置、供热管道设计等领域均会遇到。非等温条件下超静孔压演化过程与传统太沙基固结存在显著差异,Paaswel[1]提出了热固结概念用于描述非等温条件下岩土体的固结过程,目前相关研究已成为岩土涉热工程领域的重要热点。
许多学者围绕岩土体热固结过程的理论模型及计算开展了大量研究[2-11]。Zhou等[3]考虑不可逆过程效应、固结水流对流、孔隙率变化影响等诸多复杂因素,建立了饱和多孔介质热-水-力耦合过程的理论模型,推导了模型简化条件下的解析解及复杂条件下的有限元解。Mctigue[4]忽略固结对传热影响而仅考虑温度对固结作用,建立了饱和岩土体单向耦合热固结模型;针对半无限区域表面恒定温度及热流情形,通过将温度解答代入固结方程求解,推导了热固结过程解析解。Blond等[5]、白冰[6]、吴瑞潜等[7]也针对单向耦合热固结模型,采用将已有温度解代入固结方程的方式,分别推导了半无限空间、有限空间岩土体热固结过程的解析解答。当考虑固结对传热影响后,热固结模型为双向耦合,此时模型解析求解更为困难。针对双向耦合热固结模型,Bai[8]采用Laplace变换及Gauss- Legendre数值反变换对半无限空间土体在周期性温度波动下的热固结问题进行了求解;王路君等[9]通过Laplace-Hankel变换,结合数值反变换与精细积分技术,获得了衰变热源作用下层状饱和多孔介质热固结问题解答;钮家军等[10]通过势函数与分离变量法,推导了饱和单层土体在三类边界条件下热固结过程的三角级数形式解答;随后,钮家军等[11]又通过正(余)弦变换及其逆变换,建立了半无限空间土体在两类边界条件下热固结过程的精确积分解。
半无限空间土体热固结理论对于季节性地表边界作用下土体热力响应研究十分重要。目前,双向耦合情形下半无限空间土体热固结问题的求解基本都采用积分变换及其数值反变换方式,解答相对复杂,不便于应用,迫切需要新的解析方法给出形式简洁的解答。文献[12,13]针对半无限空间传热问题建立了一种相似变换方法,该方法通过引入相似变量,将偏微分方程转化为常微分方程,通过求解常微分方程获得原偏微分方程的通解,该方法特别适用于半无限空间相关问题的解析求解。因此本文针对半无限空间饱和土在外荷载及表面渗流、传热边界作用下的一维热固结过程,给出考虑双向耦合作用的热固结模型基本方程,通过引入组合变量并应用相似变换方法,建立该模型解析解答,并对热固结过程中的耦合效应进行计算分析。
1. 半无限空间饱和土热固结模型
图 1所示为半无限空间饱和土在外荷载(应力)及表面渗流、传热边界作用下一维热固结过程示意图。图 1中上表面外荷载随时间变化记为Fs(t),渗流、传热边界条件后续给出;初始时刻土体超静孔压为均匀分布Fs(0),温度为均匀分布Tini,本节给出该一维热固结过程基本控制方程。
对于固结部分的控制方程,主要简化及假设如下:①忽略土颗粒和水的压缩性,仅考虑其热膨胀性;②热渗效应仅在土体渗透性极小时才会有明显影响[3],本文不考虑;③土中孔压适用太沙基有效应力原理,固结水流服从达西定律。
在上述简化及假设基础上,可以推导获得固结部分的控制方程为[3, 7, 13]
∂P∂t−∂Fs(t)∂t+A∂T∂t=Cg∂2P∂z2。 (1) 式中:P为超静孔压;T为温度;z为位置坐标;t为时间;Cg为固结系数,Cg = kgEs/γw, kg为渗透系数,γw为水的重度,Es为土体压缩模量。
式(1)中温度变化率项的系数A及其他相关参数为
A=β−Esˉa, (2) β=E1−2vasm, (3) Es=1−v1+vE1−2v, (4) ˉa=naw+(1−n)as。 (5) 式中:β为土体热应力系数;Es为土体压缩模量;ˉa为土骨架和水的热膨胀系数平均值,E为土体变形模量;ν为泊松比;n为孔隙率;aw为水的热膨胀系数;as为土骨架热膨胀系数;asm为土体热膨胀系数,一般可以取asm=ˉa。
对于传热部分方程,主要简化及假设如下:①根据Onsager倒易关系,不考虑孔压梯度引起的传热;②过程中温增不大[2-3],热力学温度近似有Tini+∆T ≈ Tini;③忽略固结水流引起的热对流作用。
在上述简化及假设基础上,可以推导获得传热部分的控制方程为[2-3]
Cv∂T∂t+Tiniβ∂εzz∂t=λ∂2T∂z2。 (6) 式中:T为温度;Cv为土体体积热容;λ为土体导热系数;左侧第二项为土体变形功项。
为了消除式(6)中的εzz,引入土体一维线性热弹性本构[3]:
σ′zz=Esεzz−β(T−Tini)。 (7) 式中:σ′zz为土体z方向有效应力,应力、应变均以拉为正,而孔压则以压为正。
利用式(7)及有效应力原理可知:
∂εzz∂t=1Es∂P∂t+βEs∂T∂t−1Es∂Fs(t)∂t。 (8) 将式(8)代入式(6)后变形得到
(Cv+Tiniβ2Es)∂T∂t+TiniβEs[∂P∂t−∂Fs(t)∂t]=λ∂2T∂z2。 (9) 式(1),(9)构成半无限空间饱和土热固结过程温度、超静孔压耦合演化的基本方程。
为了将式(1),(9)齐次化,引入新变量:
ˉP(z,t)=P(z,t)−Fs(t)。 (10) 将新变量代入式(1),(9)中得到
∂ˉP∂t+A∂T∂t=Cg∂2ˉP∂z2, (11) ˉCv∂T∂t+B∂ˉP∂t=λ∂2T∂z2。 (12) 式中:
¯Cv=Cv+Bβ ,B=TiniβEs 。} (13) 对比B的定义式与式(6)中的土体变形功项可知,该系数是描述土体变形功对传热方程影响的一个重要参数。
对于变量ˉP与T,初始条件可以表示为
¯P(z,0)=0 ,T(z,0)=Tini 。} (14) 表面z = 0处的渗流及传热边界条件主要考虑两类,分别为给定超静孔压及温度的边界Ⅰ,以及给定水流、热流的边界Ⅱ,具体形式为
P(0,t)=f1(t) ,T(0,t)=f2(t) ,} (15) ∂P∂z|z=0=g1(t) ,∂T∂z|z=0=g2(t) 。} (16) 式中:f1(t),f2(t),g1(t),g2(t)均为关于时间的任意函数。
对于变量ˉP与T,上述两类边界可以表示为
¯P(0,t)=f1(t)−Fs(t),T(0,t)=f2(t),} (17) ∂¯P∂z|z=0=g1(t) ,∂T∂z|z=0=g2(t) 。} (18) 2. 热固结模型解析求解
2.1 组合变量解耦
定义如下形式组合变量:
Ω=ˉP+bT。 (19) 利用该组合变量,将式(11),(12)中ˉP消除后得到
∂Ω∂t+(A−b)∂T∂t=Cg∂2Ω∂z2−bCg∂2T∂z2, (20) B∂Ω∂t+(ˉCv−bB)∂T∂t=λ∂2T∂z2。 (21) 式(20)乘以λ加上式(21)乘以bCg得到:
(λ+bCgB)∂Ω∂t+[bCg(ˉCv−bB)+λ(A−b)]∂T∂t=λCg∂2Ω∂z2。 (22) 为了消除式(22)中的温度变化率项以实现解耦,b需要满足:
BCgb2+(λ−CgˉCv)b−λA=0。 (23) 式(23)是关于b的二次方程,它有两个形式的解(b1取+):
b1,2=−(λ−CgˉCv)±√(λ−CgˉCv)2+4λABCg2BCg。 (24) 式(24)给出的两个解可能是两个实数、两个共轭复数,也可能退化为同一个实数,下面先考虑b1≠b2的情形,此时式(19)定义了如下两个不同的组合变量:
Ωi=¯P+biT (i=1,2)。 (25) 两个组合变量的控制方程为
∂Ωi∂t=Ci∂2Ωi∂z2,Ci=λCgλ+biCgB(i=1,2)。} (26) 两个组合变量的初始条件为
Ωi(z,0)=¯P(z,0)+biTini=biTini(i=1,2)。 (27) 定义如下新变量消除初始条件
¯Ωi=Ωi−biTini(i=1,2)。 (28) 于是,变量ˉΩi的控制方程及初始条件为
∂¯Ωi∂t=Ci∂2¯Ωi∂z2(i=1,2), (29) ˉΩi(z,0)=0。 (30) 根据式(17),(18),ˉΩi的两类边界条件分别为
ˉΩi(0,t)=f1(t)−Fs(t)+bif2(t)−biTini, (31) ∂ˉΩi∂z|z=0=g1(t)+big2(t)。 (32) 2.2 相似变换通解
对于半无限空间中形如式(29)的控制方程,文献[12, 13]建立的相似变换法可以较方便的应用。定义如下形式的相似变量:
ψi=¯Ωi(z,t)tn/2,ηi=z2√Cit(i= 1,2)。} (33) 式中:n为任意非负整数。
将相似变量代入式(29)后推导可得常微分方程为[12]
ψ″i+2ηiψ′i−2nψi=0 (i=1,2)。 (34) 常微分方程式(34)的通解为[12]
ψi(ηi)=Ai,ninerfc(ηi)+Bi,ninerfc(−ηi)(i=1,2)。 (35) 式中:Ai, n,Bi, n为任意常数,inerfc(⋅)为高斯误差函数的累次积分函数,
i−1erfc(ε)=2√πexp(−ε2),i0erfc(ε)=erfc(ε),} (36) inerfc(ε)=∫∞εin−1erfc(x)dx。 (37) 通过相似变换通解式(35),可以获得偏微分方程式(29)如下形式的通解:
ˉΩi(z,t)=∑∞n=0tn/2[Ai,ninerfc(z2√Cit)+Bi,ninerfc(−z2√Cit)] (i=1,2)。 (38) 2.3 模型解析解答
根据inerfc(⋅)的定义式可知,通解式(38)中每个级数的第一项均满足初始时刻为0的条件,因此可以取变量ˉΩi为
¯Ωi(z,t)=∑∞n=0Ai,ntn/2inerfc(z2√Cit)(i=1,2)。 (39) 将Fs(t),f1(t),f2(t),g1(t),g2(t)均展开为时间的幂级数,并代入式(31),(32),可以得到变量ˉΩi所要满足的两类幂级数形式的边界条件为
f1(t)−Fs(t)+bif2(t)−biTini=∑mn=0ci,ntn(i=1,2), (40) g1(t)+big2(t)=∑mn=0di,ntn(i=1,2)。 (41) 式中:ci, n,di, n为幂级数的系数。
函数inerfc(·)在0处的值为
inerfc(0)=12nΓ(n2+1)。 (42) 式中:Γ(·)为伽马函数。
对式(39)及其梯度在z = 0处取值,并与式(40)、(41)进行幂级数系数对比,可以确定边界条件Ⅰ即式(40)对应的解答为
ˉΩi(z,t)=∑mn=0ci,nΓ(n+1)(4t)ni2nerfc(z2√Cit)(i=1,2)。 (43) 边界条件Ⅱ即式(41)对应的解答为
ˉΩi(z,t)=∑mn=0(−di,n)Γ(n+1)√Ci(4t)(2n+1)/2⋅i2n+1erfc(z2√Cit)(i=1,2)。 (44) 根据式(28),可以确定组合变量Ωi(z,t),随后可得温度及超静孔压解答:
T(z,t)=Ω1(z,t)−Ω2(z,t)b1−b2, (45) P(z,t)=b1Ω2(z,t)−b2Ω1(z,t)b1−b2+Fs(t)。 (46) 式(45),(46)适用于b1≠b2的情形,也包括b1,b2为两个共轭复数时。实际上,当b1,b2为两个共轭复数,式(29)中系数C1,C2也是两个共轭复数,函数inerfc(⋅)会保持复数的共轭性,ˉΩ1(z,t),ˉΩ2(z,t)也是两共轭复数,因此式(45),(46)得到的温度和超静孔压均为实数。
当(λ−CgˉCv)2+4λABCg=0时(此时A<0),b1和b2为相同实数统一记为b0,通过上面的推导只能得到一个组合变量Ω0的解答。利用该已有的组合变量,式(11),(12)可以变形为
(1−Ab0)∂ˉP∂t=Cg∂2ˉP∂z2−Ab0∂Ω0∂t, (47) (ˉCv−Bb0)∂T∂t=λ∂2T∂z2−B∂Ω0∂t。 (48) 式(47),(48)中∂ˉP/∂t与∂T/∂t的系数必有一个为正值,不妨设式(48)中为正,则式(48)为半无限空间中含热源的一维热传导问题,许多方法可以应用,这里直接采用格林函数法给出解答[14]。
对于边界条件Ⅰ,温度解答为
T(z,t)=Tini∫∞0G1(z,t|z′,0)dz′−∫t0∫∞0ε∂Ω0(z′,τ)∂τG1(z,t|z′,τ)dz′dτ+ξ∫t0f2(τ)∂G1(z,t|0,τ)∂zdτ,ξ=λ¯Cv−Bb0,ε=B¯Cv−Bb0。} (49) 式中:格林函数为[14]
G1(z,t|z′,τ)=1[4πξ(t−τ)]1/2{exp[−(z−z′)24ξ(t−τ)]−exp[−(z+z′)24ξ(t−τ)]}。 (50) 对于边界条件Ⅱ,温度解答为
T(z,t)=Tini∫∞0G2(z,t|z′,0)dz′−∫t0∫∞0ε∂Ω0(z′,τ)∂τG2(z,t|z′,τ)dz′dτ−ξ∫t0g2(τ)G2(z,t|0,τ)dτ。 (51) 式中:格林函数为[14]
G2(z,t|z′,τ)=1[4πξ(t−τ)]1/2{exp[−(z−z′)24ξ(t−τ)]+exp[−(z+z′)24ξ(t−τ)]}。 (52) 随后,超静孔压的解答可以通过组合变量定义式(19)得到。
2.4 恒定边界条件特例(b1≠b2)
(1)恒定边界条件Ⅰ
考虑外荷载为瞬时荷载Pc,土体表面超静孔压f1(t)=0,表面温度f2(t)=Ts。根据式(40)可知,展开后的幂级数仅有1项,系数为
ci,0=bi(Ts−Tini)−Pc(i=1,2)。 (53) 利用2.3节中解答可得
¯Ωi(z,t)=ci,0erfc(z2√Cit)(i=1,2)。 (54) 于是温度及超静孔压解答为
T(z,t)=Tini+c1,0erfc(z2√C1t)−c2,0erfc(z2√C2t)b1−b2, (55) P(z,t)=b1c2,0erfc(z2√C2t)−b2c1,0erfc(z2√C1t)b1−b2+Pc。 (56) (2)恒定边界条件Ⅱ
考虑外荷载为瞬时荷载Pc,土体表面水流边界g1(t)=qw,表面热流边界g2(t)=qh。根据式(41),展开后的幂级数仅有1项,系数为
di,0=qw+biqh(i=1,2)。 (57) 利用2.3节的解答,可以得到
¯Ωi(z,t)=(−di,0)√Ci(4t)1/2ierfc(z2√Cit)(i=1,2)。 (58) 于是,温度及超静孔压解答为
T(z,t)=Tini−(4t)1/2d1,0√C1ierfc(z2√C1t)−d2,0√C2ierfc(z2√C2t)b1−b2, (59) P(z,t)=Pc+(4t)1/2b2d1,0√C1ierfc(z2√C1t)−b1d2,0√C2ierfc(z2√C2t)b1−b2。 (60) 3. 单向耦合模型解析解
前面传热部分分析考虑了土体变形功项,因而式(12)中包含了超静孔压变化率,模型是双向耦合的。实际热固结计算中,忽略固结过程对传热影响的单向耦合模型也应用较多[5-7],本节考虑单向耦合情形。
单向耦合模型中传热控制方程由式(12)退化为
∂T∂t=α∂2T∂z2。 (61) 式中:热扩散率α=λ/Cv。
该情形相当于双向耦合模型中B退化为0,于是关于b的方程式(23)退化为一次方程,仅有一组解答为
bc=Aαα−Cg。 (62) 此时相当于有T以及Ωb=ˉP+bcT两个组合变量,引入:
ˉΩb=Ωb−bcTini, (63) ˉT=T−Tini。 (64) 变量ˉΩb与ˉT所满足的控制方程为
∂ˉΩb∂t=Cg∂2ˉΩb∂z2, (65) ∂ˉT∂t=α∂2ˉT∂z2。 (66) 变量ˉΩb与ˉT初始条件均为0,所需满足的边界条件Ⅰ及边界条件Ⅱ分别为
¯Ωb(0,t)=f1(t)−Fs(t)+bcf2(t)−[bcTini−Fs(0)],¯T(0,t)=f2(t)−Tini,} (67) ∂¯Ωb∂z|z=0=g1(t)+bcg2(t),∂¯T∂z|z=0=g2(t)。} (68) 式(67),(68)展开为幂级数形式如下
¯Ωb(0,t)=∑mn=0¯c1,ntn,¯T(0,t)=∑mn=0¯c2,ntn,} (69) ∂¯Ωb∂z|z=0=∑mn=0¯d1,ntn,∂¯T∂z|z=0=∑mn=0¯d2,ntn。} (70) 利用相似变换通解,经过与第2节类似的推导,可以得到单向耦合模型解析解答。对于边界条件Ⅰ,温度及超静孔压解答为
T(z,t)=∑mn=0ˉc2,nΓ(n+1)(4t)ni2nerfc(z2√αt)+Tini, (71) P(z,t)=∑mn=0ˉc1,nΓ(n+1)(4t)ni2nerfc(z2√Cgt)−bc∑mn=0ˉc2,nΓ(n+1)(4t)ni2nerfc(z2√αt)+Fs(t)−Fs(0)。 (72) 对于边界条件Ⅱ,温度及超静孔压解答为
T(z,t)=∑mn=0(−1)ˉd2,nΓ(n+1)√α(4t)(2n+1)/2⋅i2n+1erfc(z2√αt)+Tini, (73) P(z,t)=∑mn=0(−ˉd1,n)Γ(n+1)√Cg(4t)(2n+1)/2⋅i2n+1erfc(z2√Cgt)+bc∑mn=0ˉd2,nΓ(n+1)√α(4t)(2n+1)/2⋅i2n+1erfc(z2√αt)+Fs(t)−Fs(0)。 (74) 4. 算例分析
4.1 对比验证
McTigue[4]给出了半无限空间岩土体热固结过程单向耦合模型的解析解,该解答形式简单,本节利用其对2.3节中解析解作对比验证。对比算例中土体初始温度为283 K,表面温度为308 K,无外荷载作用,表面超静孔压为0。土体导热系数为1.3 W/(m·K),比热为3000 J/(kg·K),密度为1300 kg/m3;土体变形模量E为5×106 Pa,泊松比v为0.3,固结系数Cg为10-8 m2/s;土体热膨胀系数按asm=ˉa确定,孔隙率为0.3,骨架热膨胀系数为2.5×10-5 K-1,孔隙水的热膨胀系数为2×10-4 K-1。
通过2.3节解答及McTigue[4]解答对该算例进行了计算,为了使2.3节双向耦合模型接近McTigue所采用的单向耦合模型以实现对比,式(12)中无量纲参数B取了很小的值0.01。解析解编程采用Matlab进行,其中inerfc(⋅)函数可以方便的调用mfun(‘erfc’, n, z)命令。图 2,3分别给出了两组解析解在不同时刻温度及超静孔压分布的对比情况,两组解答计算结果的一致性验证了本文解析解。
4.2 耦合效应分析
本文双向耦合热固结模型控制方程式(11),(12)中,温度变化对固结部分的影响主要通过温度变化率的系数A作用,而超静孔压变化对传热部分的影响则主要通过超静孔压变化率的系数B作用,本节将研究A,B参数对热固结过程中温度、超静孔压耦合发展情况的影响。
算例中半无限空间土体初始温度为283 K,无外荷载作用,考虑式(15),(16)所描述的两类边界条件,其中的具体函数为
f1(t)=0 ,f2(t)=283+20sin(2π8760⋅t) ,} (75) g1(t)=0,g2(t)=10sin(2π8760t)。} (76) 式中:t为时间(h)。
土体基本物性参数均与4.1节相同,而采用不同参数计算时会另作说明。经试算,式(75),(76)边界条件展开为幂级数10项即可收敛,随后可以应用2.3节中解答对温度场、超静孔压场进行计算。
图 4,5给出了边界条件Ⅰ下不同深度处温度、超静孔压随时间变化的曲线。温度从初始的283 K开始呈现升高—降低—升高的典型波动趋势;随着深度的增加,温度波动的幅值逐渐减小,而出现峰值的时间有所推后,这是温度波传播过程中典型的衰减和延迟性质。超静孔压的变化则是从初始的0开始,呈现降低—升高—降低的波动趋势,总体与温度的波动相反;图中4个深度超静孔压波动辐值最高出现在z=0.5 m,而不是更接近地表的z=0.2 m,这是由于地表超静孔压为0,太接近地表反而限制了超静孔压的波动。对比B为0.01,100两种情形发现,B值的变化对温度、超静孔压的发展基本没有影响。
图 6,7给出了边界条件Ⅱ下不同深度处温度、超静孔压随时间变化的曲线。由于式(76)给出的表面热流是先从地表向外流动的,温度变化主要呈降低—升高—降低的波动趋势,随着深度的增加,温度波动仍体现衰减和延迟特征。超静孔压则呈升高—降低—升高的波动趋势,同样与温度的波动相反;随着深度的增加,超静孔压的波动也出现了衰减和延迟,这是不透水边界与图 4,5中透水边界的区别。B=100时的结果与B=0.01时基本无差异(未作曲线);而B=1000时与B=0.01时结果相比有一些可见的差异,温度及超静孔压的波动幅度均有所减小。
图 4~7中选取不同B是为了研究式(12)中土体变形功对传热及固结过程的影响;若直接根据4.1节参数及式(13),可以确定B=4.07×10-2。从计算结果看,B值的变化对温度、超静孔压发展的影响较小,这一结论也可以从式(12)左侧两项的量级分析中获得。式(12)左侧温度变化率项的系数为Cv+Bβ, 其中Cv的大小是3.9×106 J/(m3·K),而Bβ的大小即使在B=100时也仅为0.09×106 J/(m3·K),是远小于Cv的;同时,孔压变化率项的系数B (即使取100)也远小于Cv,这就造成了图 4~7中B的变化对计算结果的影响较小。根据式(3),(4),(13),B不会超过0.3,B,Bβ相对于Cv都是很小的;因此,除非超静孔压变化率很高(102 Pa/s量级),土体变形功项对温度、超静孔压的影响基本都可以忽略。
文献[13]对参数A进行过分析,指出按式(2)及asm=ˉa取值时(上面算例即是如此),A始终为正值,而也有文献未按此确定asm及A,得到的A可能为负值[3]。下面具体计算分析A的正、负对温度、超静孔压耦合演化过程的影响。
算例参数也均与4.1节相同,通过设置不同A来分析温度变化率项对热固结过程的影响。图 8给出了两种边界情况下不同深度处温度随时间变化的曲线,图中温度变化趋势与图 4,6中的结果是一致的;参考文献[13]中总结的A的取值,选择了A为500,-500 Pa/K两种情形,图中结果表明,A的变化对两种情形下温度场的演化基本没有影响。
图 9给出了边界Ⅰ情形下不同深度超静孔压随时间的变化情况。在A=500 Pa/K时,超静孔压随时间的变化总体与温度随时间的变化呈相反趋势,这与图 4~7中的结果是一致的。然而,在A=-500 Pa/K时,超静孔压随时间的变化总体与温度变化呈相同的趋势;以z=1.0 m位置为例,t<107 s,超静孔压及温度总体呈升高趋势,而在107 s<t<2.6×107 s,超静孔压与温度总体呈降低趋势,最后在t>2.6×107 s,超静孔压及温度转变为升高趋势。
图 10给出了边界Ⅱ情形下不同深度超静孔压随时间的变化情况。与图 8中温度变化对比可知,在A= 500 Pa/K时,超静孔压的变化总体与温度变化呈相反趋势,而在A=-500 Pa/K时,超静孔压的变化总体与温度变化呈相同趋势,这与图 9中的结果是一致的。
A的正负对超静孔压与温度耦合演化的影响特征可以通过对式(11)的分析进一步阐明。对比式(11)与传统固结方程可知,温度变化项的影响相当于给固结过程增加了一个源项-A∂T/∂t,温度变化即通过该源项引起相应的超静孔压变化。当A>0时,∂P/∂t与∂T/∂t前面的符号相反,温度变化引起的超静孔压变化与之呈相反趋势;当A<0时, ∂P/∂t与∂T/∂t前面的符号相同,温度变化引起的超静孔压变化与之呈相同趋势,这便是图 8~10中所展示的结果。
从式(11)也可以看出,若将A = 500 Pa/K情形下该式的解答记为P(z, t),由于A的变化对温度解答基本没有影响,A=-500 Pa/K情形下该式的解答则为-P(z, t),即A = -500 Pa/K情形下的超静孔压与A = 500 Pa/K情形下的超静孔压是关于直线P = 0对称的,这与图 9,10中展示的结果一致。
5. 结论
针对半无限空间土体在外荷载及表面两类渗流、传热边界条件下的热固结问题,给出了描述温度、超静孔压双向耦合演化过程的基本方程,并基于相似变换方法建立了幂级数形式荷载及渗流、传热边界条件下温度与超静孔压的解析解答。多数情形下,解答可以直接用inerfc(⋅)函数的代数组合形式给出,应用十分方便。
(1)温度波动从地表向深处传播有明显的衰减和延迟特性,而超静孔压传播则随渗流边界条件不同而有所差异。地表为不透水边界时,超静孔压波动随深度增加也具有衰减和延迟特性;地表为透水边界时则有所不同,超静孔压波动幅度最大值出现在距地表一定深度处。
(2)传热控制方程中土体变形功项的贡献很小,模型计算基本可以忽略土体变形功对温度、超静孔压演化过程的影响。
(3)固结控制方程中温度变化率项对温度、超静孔压的耦合演化有显著影响,其特征主要决定于该项的系数A。当系数A为正时,温度变化将引起相反趋势的超静孔压变化;当系数A为负时,温度变化将引起相同趋势的超静孔压变化。
-
-
[1] PAASWELL R E. Temperature effects on clay soil consolidation[J]. Journal of the Soil Mechanics and Foundations Division, 1967, 93(3): 9-22. doi: 10.1061/JSFEAQ.0000982
[2] BIOT M A. Thermoelasticity and irreversible thermodynamics[J]. Journal of Applied Physics, 1956, 27(3): 240-253. doi: 10.1063/1.1722351
[3] ZHOU Y, RAJAPAKSE R K N D, GRAHAM J. A coupled thermoporoelastic model with thermo-osmosis and thermal-filtration[J]. International Journal of Solids and Structures, 1998, 35(34/35): 4659-4683.
[4] MCTIGUE D F. Thermoelastic response of fluid-saturated porous rock[J]. Journal of Geophysical Research: Solid Earth, 1986, 91(B9): 9533-9542. doi: 10.1029/JB091iB09p09533
[5] BLOND E, SCHMITT N, HILD F. Response of saturated porous media to cyclic thermal loading[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27(11): 883-904. doi: 10.1002/nag.301
[6] 白冰. 循环温度荷载作用下饱和多孔介质热-水-力耦合响应[J]. 工程力学, 2007, 24(5): 87-92. doi: 10.3969/j.issn.1000-4750.2007.05.015 BAI Bing. Thermo-hydro-mechanical responses of saturated porous media under cyclic thermal loading[J]. Engineering Mechanics, 2007, 24(5): 87-92. (in Chinese) doi: 10.3969/j.issn.1000-4750.2007.05.015
[7] 吴瑞潜, 谢康和, 程永锋. 变荷载下饱和土一维热固结解析理论[J]. 浙江大学学报(工学版), 2009, 43(8): 1532-1537. WU Ruiqian, XIE Kanghe, CHENG Yongfeng. Analytical theory for one-dimensional thermal consolidation of saturated soil under time-dependent loading[J]. Journal of Zhejiang University (Engineering Science), 2009, 43(8): 1532-1537. (in Chinese)
[8] BAI B. Fluctuation responses of saturated porous media subjected to cyclic thermal loading[J]. Computers and Geotechnics, 2006, 33(8): 396-403. doi: 10.1016/j.compgeo.2006.08.005
[9] 王路君, 艾智勇. 衰变热源作用下饱和多孔介质热固结问题的扩展精细积分法[J]. 力学学报, 2017, 49(2): 324-334. WANG Lujun, AI Zhiyong. Epim for thermal consolidation problems of saturated porous media subjected to a decaying heat source[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(2): 324-334. (in Chinese)
[10] 钮家军, 凌道盛, 王秀凯, 等. 饱和单层土体一维热固结精确解[J]. 岩土工程学报, 2019, 41(9): 1715-1723. doi: 10.11779/CJGE201909016 NIU Jiajun, LING Daosheng, WANG Xiukai, et al. Exact solutions for one-dimensional thermal consolidation of single-layer saturated soil[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(9): 1715-1723. (in Chinese) doi: 10.11779/CJGE201909016
[11] 钮家军, 王秀凯, 凌道盛, 等. 半无限空间饱和土一维热固结精确积分解[J]. 岩土工程学报, 2021, 43(8): 1426-1433. doi: 10.11779/CJGE202108007 NIU Jiajun, WANG Xiukai, LING Daosheng, et al. Exact integral solutions for one-dimensional thermal consolidation of semi-infinite saturated soils[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(8): 1426-1433. (in Chinese) doi: 10.11779/CJGE202108007
[12] ZHOU Y, ZHANG L Y, XU C, et al. Analytical solution for classical one-dimensional thaw consolidation model considering unfrozen water effect and time-varying load[J]. Computers and Geotechnics, 2020, 122: 103513. doi: 10.1016/j.compgeo.2020.103513
[13] 周扬, 武子寒, 许程, 等. 高温下饱和冻土一维融化热固结模型及解答[J]. 岩土工程学报, 2021, 43(12): 2190-2199. doi: 10.11779/CJGE202112005 ZHOU Yang, WU Zihan, XU Cheng, et al. One-dimensional thaw thermo-consolidation model for saturated frozen soil under high temperature and its solution [J]. Chinese Journal of Geotechnical Engineering, 2021, 43(12): 2190-2199. (in Chinese) doi: 10.11779/CJGE202112005
[14] COLE K D, BECK J V, HAJI-SHEIKH A, et al. Heat Conduction Using Green's Functions[M]. 2nd ed. Boca Raton: CRC Press, 2010.
-
其他相关附件