EPIM FOR THERMAL CONSOLIDATION PROBLEMS OF SATURATED POROUS MEDIA SUBJECTED TO A DECAYING HEAT SOURCE
-
摘要: 热源作用下饱和多孔介质热固结效应是土木及能源工程领域的一个重要课题.由于问题的复杂性，已有的研究大多将介质假定为均匀各向同性，且将热源假定为恒定强度.实际工程中，天然饱和多孔介质常表现出明显的分层特性，热源强度也存在衰变性，为此本工作采用扩展精细积分法对衰变热源作用下层状饱和多孔介质的热固结问题进行研究.借助于积分变换，将饱和多孔介质热固结问题的偏微分方程转化为变换域内的常微分方程；然后对饱和多孔介质微层元进行合并消元，并结合边界条件，推导出衰变热源作用下层状饱和多孔介质热固结问题在积分变换域内的扩展精细积分解；对所得解答进行相应的数值积分逆变换，可获得所求温度、超静孔压及竖向位移在物理域内的解答.基于上述求解过程，编制相应的计算程序进行数值计算，通过与已有文献对比，验证本文扩展精细积分法在求解层状饱和多孔介质热固结问题中的适应性和正确性；最后通过几组算例，分析热源衰变周期、热源埋深及介质的成层性对热固结效应的影响.结果表明：热源衰变周期对温度和超静孔压的峰值、以及达到峰值的时间均有明显影响，衰变周期越长，二者峰值均越大，且达到峰值所需时间越长；热源埋深对超静孔压及竖向位移变化影响显著，深埋热源作用时热源两侧竖向位移呈对称分布，而浅埋热源两侧则无此现象；饱和多孔介质的分层特性对热固结效应影响明显.Abstract: The thermal consolidation of saturated porous media subjected to a heat source is an important subject in civil engineering and energy engineering. For the complexity of the problem, the porous media are usually treated as homogeneous isotropic media and the heat source is assumed to be a heat source with constant strength in the existing studies. In engineering practice, natural saturated porous media usually show obvious layered characteristics and the heat source is decaying with time. In this case, the extended precise integration method (EPIM) is presented in this study to investigate the thermal consolidation problems of layered saturated porous media subjected to a decaying heat source. The partial differential equations are reduced to ordinary ones by means of the integral transform techniques. Combining the adjacent layer elements and considering the boundary conditions, the EPIM solutions in the transformed domain of the problems are deduced. With the aid of corresponding numerical integral inversion, the temperatures, excess pore pressures and vertical displacements in the physical domain are obtained. A numerical example with the corresponding calculation program is performed to compare with the existing results, which confirm the applicability and validity of the presented method in dealing with the thermal consolidation problems of layered saturated porous media. Finally, numerical examples are carried out to analyse the influence of the heat source's half-life and buried depth, as well as the stratification of medium on the thermal consolidation behaviour. Numerical results show that:the decay period of heat sources has significant influence on the peak values and peak time of temperature and excess pore pressure, the longer the decay period, the greater the peak values and the longer the peak time of temperature and excess pore pressure; burial depths have obvious influence on the variations of excess pore pressure and vertical displacement, the evolutions of vertical displacements against time on both side of the deeply buried heat source are symmetrical, while there is no such phenomenon for the shallow heat source; stratification characteristics of the saturated porous media shows prominent effects on the thermal consolidation.
-
引言
热源作用下饱和多孔介质热固结问题是土木、能源、环境等领域的热点课题之一.饱和多孔介质中，热源作用会引起孔隙水压力的产生及消散，而孔隙水压力的变化会对固体骨架变形及强度产生重要影响.因此，热固结问题是渗流场、温度场和应力应变场相互耦合作用的结果，该课题在放射性核废料地下处置、地热资源开发、能量桩等工程中具有重要的研究和应用价值[1-5].在Biot固结理论[6]的基础上，Biot[7]建立了饱和多孔介质热-水-力耦合问题的基本控制方程，并对方程中的耦合系数进行了明确的定义和完整的物理解释，较早开展了饱和多孔介质热固结响应的研究.由于Biot[7]所建立的热-水-力耦合理论是位移、孔压及温度变化耦合的复杂偏微分方程，故在求解上遇到了一定的困难.因此，国内外诸多学者对该理论进行简化和修正，并在求解方法上进行了诸多探索.Booker和Savvidou[8]，Savvidou和Booker[9]分析了饱和岩土介质热固结问题的作用机理，并给出了球形和点热源作用下热固结问题的解析解.Mctigue[10]对热源作用于半无限多孔饱和介质表面的热-水-力耦合问题进行了研究，并得到了其解析解.Bai和Abousleiman[11]基于热-水-力耦合理论，探讨了完全耦合、部分耦合及完全非耦合三种理论的适用条件，并给出了三种情况下的解析解答.白冰等[4, 12]基于饱和多孔介质热-水-力耦合控制方程，借助于积分变换对热弹性固结模型进行求解，并给出温度、超静孔压和位移演化过程的解析式.郑荣跃等[13]研究了半无限地基在内置点热力源作用下的响应问题，并给出热力源作用下应力、位移、孔隙水压力的解.吴瑞潜等[14]建立了变载荷作用下饱和土体一维热固结问题的解析解.Lu和Lin[15]对衰变热源作用下多孔弹性半空间的热固结问题进行了研究，并给出了解析解.Selvadurai和Suvorov[16]对固体骨架为Hooke弹性或弹塑性体的饱和介质进行研究，并分析了不同条件下的热-水-力耦合响应.
天然介质经过长期的沉积过程，往往表现出明显的分层特性.为更真实地描述热固结过程，部分学者开展了层状介质热固结问题的研究.Giraud等[17]推导了两层介质一维热固结问题的半解析解答.白冰[18]采用解析方法，对变温载荷作用时的双层半无限饱和介质的热固结问题进行了研究.Ai和Wang[5]采用解析层元法求解了层状饱和介质轴对称热固结问题.目前针对层状弹性体系，常用的求解方法有传递矩阵法[19-21]、有限层法[22-23]、刚度矩阵法或解析层元法[5, 24-26]等.对于层状弹性体及饱和地基等本构方程较为简单的情况，采用上述方法往往可推导出其显式解析解.但对于较复杂的求解模型，如横观各向同性体及热-水-力耦合问题，其控制方程较复杂，极难采用上述方法推导出其显式解析解.另一方面，目前主流的数值计算方法，如有限元法等，在求解热固结问题时需耗费大量时间和内存，且需较强的经验性.鉴于此，本文拟基于具有稳定性好、计算效率与精度高特点的精细积分法[27-29]，并结合积分变换，建立衰变热源作用下层状饱和多孔介质热固结问题的扩展精细积分解.与文献[28]所给出的精细积分法相比，本文方法不仅可以求解表面受荷问题，而且还可以求解层状介质内部任意位置的受荷问题；另外，本文方法的整个求解过程均在积分变换域内进行，是精细积分法[27-28]与积分变换法相结合的一个新尝试[30-31].由于本文方法是精细积分法[28]应用的一个扩展，为此本文将其称之为"扩展精细积分法".首先推导出热固结问题在Laplace-Hankel变换域内的常微分矩阵方程；然后对介质微层元进行合并消元，得到该问题在变换域内的扩展精细积分解；最后对该解进行相应的积分逆变换，可得其在物理域内的解答.
1. 热固结问题的常微分矩阵方程
轴对称条件下，不考虑体力时，饱和多孔介质热固结问题的平衡方程为
$${\frac{{\partial {\sigma _r}}}{{\partial r}} + \frac{{\partial {\sigma _{rz}}}}{{\partial z}} + \frac{{{\sigma _r}-{\sigma _\varphi }}}{r} = 0}$$ (1a) $${\frac{{\partial {\sigma _z}}}{{\partial z}} + \frac{{\partial {\sigma _{rz}}}}{{\partial r}} + \frac{{{\sigma _{rz}}}}{r} = 0}$$ (1b) 式中，$\sigma _r$, $\sigma _\varphi$, $\sigma _z$分别为r, $\varphi$, z方向的正应力，$\sigma _{rz}$为$r - z$面上的剪应力.
结合广义热弹性Hooke定律和有效应力原理，可得到用位移和温度增量表示的本构关系
$${{\sigma _r} = \lambda {\varepsilon _v} + 2G\frac{{\partial {u_r}}}{{\partial r}}-\beta \theta-\sigma }$$ (2a) $${{\sigma _\varphi } = \lambda {\varepsilon _v} + 2G\frac{{{u_r}}}{r}-\beta \theta-\sigma }$$ (2b) $${{\sigma _z} = \lambda {\varepsilon _v} + 2G\frac{{\partial {u_z}}}{{\partial z}}-\beta \theta-\sigma }$$ (2c) $${{\sigma _{rz}} = G\left( {\frac{{\partial {u_r}}}{{\partial z}} + \frac{{\partial {u_z}}}{{\partial r}}} \right)}$$ (2d) 式中，$\varepsilon _v = \dfrac{\partial u_r}{\partial r} + \dfrac{u_r }{r} + \dfrac{\partial u_z }{\partial z}$，其中，$u_r$和$u_z$分别表示r和z方向的位移；$\lambda =2G\mu \eta$为Lamé常数，$\eta = 1 /{(1 - 2\mu })$, $G = E/{2(1 + \mu )}$为剪切模量，E和$\mu$分别表示杨氏模量和泊松比，$\beta = 2G\alpha \eta (1 + \mu)$，$\theta$和$\sigma$分别为温度增量和超静孔隙压力，$\alpha$为固体骨架的线膨胀系数.
热量在介质内的传导符合Fourier导热定律，则0至t时间内，z方向的热流量$Q_\theta$可定义为
$$Q_\theta = - \int_0^t {K\frac{\partial \theta}{\partial z}{\rm d}t}$$ (3) 式中，K为热传导系数.
根据Fourier定律与能量守恒方程，热传导方程可表示为[9]
$$\frac{\partial \theta }{\partial t} = \kappa {{\nabla}}^2\theta$$ (4) 式中，$\kappa = K /m$，$m = \rho C$，$\rho$和C分别表示密度和比热容，${\nabla}^2 = \dfrac{\partial^2}{\partial r^2} + \dfrac{1}{r}\dfrac{\partial }{\partial r} +\dfrac{\partial ^2}{\partial z^2}$.
据Darcy定律，初始时刻至t时刻竖向流量为
$$Q_f = - \int_0^t{\frac{c}{\gamma _w }\frac{\partial \sigma }{\partial z}{\rm d}t}$$ (5) 式中，c为渗透系数.
结合Darcy定律和渗流连续条件，渗流连续方程可表示为
$$\frac{\partial\varepsilon _v }{\partial t} + \alpha _u \frac{\partial \theta}{\partial t} = \frac{c}{\gamma _w }{\nabla}^2\sigma$$ (6) 式中，$\alpha _u = 3\alpha _{\rm s} (n - 1) - 3n\alpha _{\rm f}$，$\alpha _{\rm s}$和$\alpha _{\rm f}$分别为固体颗粒和孔隙水的线膨胀系数，n表示孔隙率.
由于直接求解上述偏微分方程较困难，本文借助于积分变换方法，将偏微分方程转化为易于求解的常微分方程.函数$f(r, z, t)$关于时间t的Laplace变换及其逆变换为[32]
$${\tilde f(r, z, s) = \int_0^{ + \infty } {f(r, z, t){{\rm {e}}^{-st}}{\rm {d}}t} }$$ (7a) $${f(r, z, t) = \frac{1}{{2\pi {\rm {i}}}}\int_{a-{\rm {i}}\infty }^{a + {\rm {i}}\infty } {\tilde f(r, z, s){{\rm {e}}^{st}}{\rm {d}}s} }$$ (7b) 式中，s表示关于t的Laplace变换参数，${\rm i} = \sqrt { - 1}$.
函数$\tilde{f}(r, z, s)$关于坐标r的m阶Hankel变换及其逆变换定义为[33]
$${{{\bar f}^m}(\xi, z, s) = \int_0^{ + \infty } {\tilde f(r, z, s){{\rm{J}}_m}(\xi r)r{\rm{d}}r} }$$ (8a) $${\tilde f(r, z, s) = \int_0^{ + \infty } {{{\bar f}^m}(\xi, z, s){{\rm{J}}_m}(\xi r)\xi {\rm{d}}\xi } }$$ (8b) 式中，$\xi$是关于r的Hankel变换参数，${\rm J}_m (\xi r)$为m阶Bessel函数.
对式 (1)~式 (6) 进行关于时间t和坐标r的Laplace及Hankel变换，所得常微分方程表示为矩阵形式
$$\frac{\mbox{d}}{\mbox{d}z}\left[{{\begin{array}{*{20}c}{{\varGamma}} \\{\varLambda}\\\end{array} }} \right] = \left[{{\begin{array}{*{20}c} {{{\varPhi}}_1 }&{{{\varPhi}} _2 } \\ {{{\varPhi}}_3 }&{{{\varPhi}} _4 } \\\end{array} }} \right] \cdot \left[{{\begin{array}{*{20}c} {{\varGamma}} \\ {{\varLambda}}\\\end{array} }} \right]$$ (9) 式中，${\varGamma}= \left[{\bar {\sigma }_{rz}^1, \bar {\sigma}_z^0, \bar {\sigma }^0, \bar {\theta }^0} \right]^{\rm T}$为变换域内的广义应力向量，${\varLambda} = \left[{\bar{u}_r^1, \bar {u}_z^0, \bar {Q}_f ^0, \bar {Q}_\theta ^0}\right]^{\rm T}$为广义位移向量，上标1和0表示对该量进行1阶和0阶Hankel变换
$${\varPhi} _1 \! = \!\left[\!\!{{\begin{array}{*{20}c} 0 &{b_2 \xi } &{b_3 \xi } &{b_6 \xi } \\ {-\xi } &0 &0 &0 \\ 0 &0 &0 &0 \\ 0 &0 &0 &0 \\\end{array} }} \!\!\right], {\varPhi} _2 \!\mbox{ = }\!\left[\!\!{{\begin{array}{*{20}c} {b_5 \xi ^2} &0 &0 &0 \\ 0 &0 &0 &0 \\ 0 &0 &{f_1 } &0 \\ 0 &0 &0 &{f_2 } \\\end{array} }} \!\!\right] \\ {\varPhi} _3 \!\mbox{ = }\!\left[\!\!{{\begin{array}{*{20}c} {b_7 } &0 &0 &0 \\ 0 &{b_1 } &{-b_1 } &{b_4 } \\ 0 &{-b_1 } &{g_1 } &{g_2 } \\ 0 &0 &0 &{g_3 } \\\end{array} }} \!\!\right], {\varPhi} _4 \!\mbox{ = }\!\left[\!\!{{\begin{array}{*{20}c} 0 &\xi &0 &0 \\ {-b_2 \xi } &0 &0 &0 \\ {b_3 \xi } &0 &0 &0 \\ 0 &0 &0 &0 \\\end{array} }} \!\!\right]$$ 其中，$b_1 = \dfrac{1 - 2\mu }{2G(1 - \mu )}, ~ b_2= \dfrac{\mu }{1 - \mu }, ~ b_3 = - \dfrac{1 - 2\mu}{1 - \mu }b, ~ b_4 = \dfrac{1 + \mu }{1 - \mu }\alpha,$$ ~ b_5 =\dfrac{2G}{1 - \mu }, ~ b_6 = - 2Gb_4,$$ ~ b_7 =\dfrac{1}{G}, ~ f_1 = - {s\gamma _w }/c, ~f_2 = - s /K, $$~ g_1 = -b_1 b^2 - {c\xi ^2} /{s\gamma _w }, ~ g_2 = 3\alpha _u - b_4, ~ g_3= - C\rho - {K\xi ^2} /s.$
2. 层状体系的精细积分法
对于两点边值问题的精细积分法，文献[28]给出了严格的理论推导.对于z向任意微层元$\left[{z_a, z_b }\right]$，其上下界面状态向量如图 1.
对于式 (9) 所描述的两端边值问题，其两端状态向量存在如下关系[28]
$${\varGamma} _b = { F}{\varGamma}_a - { G}{\varLambda}_b$$ (10a) $${\varLambda}_a = { Q}{\varGamma}_a +{ E}{\varLambda}_b$$ (10b) 式中，${ F}$, ${ G}$, ${ Q}$, ${E}$为待求关系矩阵，它们建立了微层元两端位移和应力向量间的关系，均为关于$z_a$，$z_b$的函数.若层元厚度$\ell = \Delta z\mbox{ = }z_b - z_a$非常小，则关系矩阵${ F}$, ${ G}$, ${ Q}$, ${E}$可进行Taylor级数展开，此时其表达式中仅包含层元厚度$\ell$及矩阵${\varPhi}_i (i = 1, 2, 3, 4)$，而矩阵${\varPhi} _i$只与介质的材料参数有关，即层元厚度$\ell$和介质材料参数一经确定，则关系矩阵${ F}$, ${ G}$, ${Q}$, ${ E}$可推导求出，具体可参考文献[34].
式 (10) 建立了任意微层元上下界面状态变量之间的关系.对于任意两相邻微层元1和2，如图 2所示.
根据交界面处变量的协调条件，可以推导出$z_a$和$z_c$界面上状态变量间的关系式[34]
$${\varGamma}_c = { F}_3 {\varGamma} _a -{ G}_3 {\varLambda}_c$$ (11a) $${\varLambda}_a = { Q}_3 {\varGamma} _a + { E}_3 {\varLambda}_c$$ (11b) 式中，${ F}_3$, ${ G}_3$, ${ Q}_3$, ${ E}_3$为合并后形成的微层元3的关系矩阵，其具体表达式为
$${{F_3} = {F_2}{{(I + {G_1}{Q_2})}^{-1}}{F_1}}$$ (12a) $${{G_3} = {G_2} + {F_2}{{(G_1^{-1} + {Q_2})}^{-1}}{E_2}}$$ (12b) $${{Q_3} = {Q_1} + {E_1}{{(Q_2^{-1} + {G_1})}^{-1}}{F_1}}$$ (12c) $${{E_3} = {E_1}{{(I + {Q_2}{G_1})}^{-1}}{E_2}}$$ (12d) 由层元合并消元次序无关定理[28]知，式 (12) 可看作微层元消元合并的一个递归表达式.若将微层元3与其相邻的微层元4继续合并形成微层元5，并将微层元序号3, 4, 5依次用1, 2, 3替换，则微层元5上下界面处状态向量间关系仍可用式 (11) 和式 (12) 表述.可见，对于多层介质中的任一层元，将其划分为多个微层元，再进行微层元间的合并操作，最终可建立该层上下表面状态变量间的关系.
由以上推导可知，只要$\ell$足够小，则数值计算结果的误差是微小的.但当$\ell$过度小时，计算时可能会因计算机存储精度问题而导致有效位数丢失.为避免该类问题，递归式 (12) 需采用如下表达式
$$\left.\begin{array}{lll} { F}_3^\# = ({ F}_2^\# - { G}_1 {Q}_2 / 2){ J} +{ F}_2^\# {{JF}}_1^\# +\\\quad { J}({ F}_1^\# - { G}_1{ Q}_2 / 2)\\{ G}_3 = { G}_2 + ({ I} + { F}_2^\# )({ G}_1^{ - 1}+ { Q}_2 )^{ - 1}({ I} + { E}_2^\# )\\{ Q}_3 = { Q}_1 + ({ I} + { E}_1^\# )({ Q}_2^{ - 1}+ { G}_1 )^{ - 1}({ I} + { F}_1^\# )\\{ E}_3^\# = ({ E}_1^\# - { Q}_2 { G}_1 / 2){ K} +{ E}_1^\# {{KE}}_2^\#+ \\\quad { K}({ E}_2^\# - { Q}_2{ G}_1 / 2)\\{ F}_i = { I} + { F}_i^\#, \quad { E}_i = { I} +{ E}_i^\#, (i =1, 2, 3)\end{array}\right\}$$ (13) 式中，${ J} =({ I} + { G}_1 { Q}_2 )^{ - 1}$，${ K} = ({ I} +{ Q}_2 { G}_1 )^{ - 1}$.式 (12) 与式 (13) 实际上是相同的，只是将原来的${ F}$, ${E}$分别表示为${ I}$与${ F}^\#$, ${ E}^\#$之和的形式，这是由于$\ell$很小时，${ F}^\#$, ${ E}^\#$与${ I}$相比为极小的矩阵.若直接用${ F}$和${E}$进行计算，${ F}^\#$和${ E}^\#$可能会因计算机存储精度问题被消去而得不到精确结果.故递归操作时仅涉及${ F}^\#$和${ E}^\#$运算，以避免发生精度丢失的问题.
3. 内部受荷问题的扩展精细积分法
如图 3所示，对于内部作用有热源的多层介质体系，根据热源作用深度$H_F$和计算点深度$H_C$，可将多层体系分成三个部分，为便于表述，将其称之为"层块".三个层块分别为：层块1$\left[{a, b}\right]$，介质体系表面与热源作用 (或计算点) 深度间的层块；层块2$\left[{b, c} \right]$，热源作用深度与计算点间的层块；层块3$\left[{c, d}\right]$，计算点 (或热源作用) 深度与介质体系底面间的层块.各层块可能是单个自然层，也可以是某个自然层的一部分，或是由多个自然层组合而成.如图中虚线所示，各层块内部微层元采用式 (12) 或式 (13) 进行消元凝聚，得到各自的关系矩阵，分别记为$\hat{{ F}}_i$, $\hat {{ G}}_i$, $\hat {{ Q}}_i$, $\hat{{ E}}_i (i = 1, 2, 3)$.根据$H_F$和$H_C$的关系，可分不同情况进行讨论.
当$H_F > H_C$时，各层块的关系矩阵方程如下：
层块1
$$\left.\begin{array}{lll}{\varGamma}_b = \hat {{ F}}_1 {\varGamma} _a - \hat{{ G}}_1 {\varLambda}_b \\{\varLambda} _a = \hat {{ Q}}_1 {\varGamma}_a + \hat{{ E}}_1 {\varLambda}_b\end{array}\right\}$$ (14) 层块2
$$\left.\begin{array}{lll}{\varGamma}_c^ - =\hat {{ F}}_2 {\varGamma}_b - \hat {{ G}}_2{\varLambda}_c^ - \\{\varLambda}_b = \hat {{ Q}}_2 {\varGamma}_b + \hat {{E}}_2 {\varLambda}_c^ -\end{array}\right\}$$ (15) 层块3
$$\left.\begin{array}{lll}{\varGamma}_d =\hat {{ F}}_3 {\varGamma}_c^ + - \hat {{ G}}_3{\varLambda}_d \\{\varLambda}_c^ + = \hat {{ Q}}_3 { \varGamma}_c^ + +\hat {{ E}}_3 {\varLambda}_d\end{array}\right\}$$ (16) 式中，${\varGamma}_c^ -$和${\varLambda}_c^ -$为层块2在$H_F$处的广义应力和位移，${\varGamma}_c^ +$和${\varLambda} _c^+$为层块3在$H_F$处的广义应力和位移.
界面c处可能作用有广义载荷${ P}_\varGamma$或广义位错${P}_\varLambda$，此时有如下关系
$$\left.\begin{array}{lll}{\varGamma} _c^ + ={\varGamma} _c^ - - { P}_{\varGamma}\\{\varLambda} _c^ + = {\varLambda} _c^ - - { P}_\varLambda\end{array}\right\}$$ (17) 结合式 (15)~式 (17)，消去c界面的状态向量可得
$$\left.\begin{array}{lll}{\varGamma}_d = \hat {{ F}}_{23} {\varGamma}_b- \hat {{ G}}_{23} {\varLambda} _d - \hat {{ H}}_3{ P}_\varGamma - \hat {{ L}}_3 { P}_\varLambda\\{\varLambda}_b = \hat {{ Q}}_{23} {\varGamma} _b + \hat{{ E}}_{23} {\varLambda}_d - \hat {{ H}}_4 {P}_\varGamma + \hat {{ L}}_4 { P}_\varLambda\end{array}\right\}$$ (18) 式中
$$\left.\begin{array}{llll}\hat {{ F}}_{23}= \hat {{ F}}_3 ({ I} + \hat {{ G}}_2 \hat {{ Q}}_3)^{ - 1}\hat {{F}}_2\\\hat {{ G}}_{23} = \hat {{ G}}_3 + \hat {{ F}}_3 (\hat{{ G}}_2^{ - 1} + \hat {{ Q}}_3 )^{ - 1}\hat {{ E}}_3\\\hat {{ E}}_{23} = \hat {{ E}}_2 ({ I} + \hat {{ Q}}_3\hat {{ G}}_2 )^{ - 1}\hat {{ E}}_3\\\hat {{ Q}}_{23} = \hat {{ Q}}_2 + \hat {{ E}}_2 (\hat{{ G}}_2 + \hat {{ Q}}_3^{ - 1} )^{ - 1}\hat {{ F}}_2\\\hat {{ H}}_3 = \hat {{ F}}_3 ({ I} + \hat {{ G}}_2\hat {{ Q}}_3 )^{ - 1}, \quad \hat {{ H}}_4 = \hat {{L}}_4 \hat {{ Q}}_3\\\hat {{ L}}_3 = \hat {{ H}}_3 \hat {{ G}}_2, \quad \hat{{ L}}_4 = \hat {{ E}}_2 ({ I} + \hat {{ Q}}_3 \hat{{ G}}_2 )^{ - 1}\end{array}\right\}$$ (19) 由式 (14) 和式 (18)，可得计算点界面b处广义应力和位移的表达式
$$\left.\begin{array}{lll}{\varGamma}_b = \hat {{ F}}_1 {\varGamma}_a -\hat {{ G}}_1 {\varLambda}_b\\{\varLambda}_b = ({ I} + \hat {{ Q}}_{23} \hat {{G}}_1 )^{ - 1}(\hat {{ Q}}_{23} \hat {{ F}}_1 {\varGamma}_a - \hat {{ E}}_{23} {\varLambda}_d - \\\quad \hat {{ H}}_4 { P}_\varGamma +\hat {{ L}}_4 {P}_\varLambda )\end{array}\right\}$$ (21) 当$H_F < H_C$时，采用类似的推导方法，可建立计算点界面c处应力和位移的表达式
$$\left.\begin{array}{ll}{\varGamma}_c = ({ I} + \hat {{ G}}_{12} \hat{{ Q}}_3 )^{ - 1}(\hat {{ F}}_{12} {\varGamma} _a - \hat{ G}_{12} \hat { E}_3 {\varLambda}_d - \\\quad\hat {{ H}}_1{ P}_\varGamma - \hat {{ L}}_1 { P}_\varLambda )\\{\varLambda} _c = \hat {{ Q}}_3 {\varGamma} _c + \hat{ E}_3 {\varLambda}_d\end{array}\right\}$$ (22) 式中，关系矩阵$\hat {{F}}_{12}$, $\hat {{ G}}_{12}$, $\hat {{ Q}}_{12}$, $\hat{{ E}}_{12}$, $\hat{{ H}}_1$, $\hat{ H}_2$, $\hat {{L}}_1$, $\hat {{ L}}_2$与式 (18) 中$\hat {{ F}}_{23}$, $\hat {{ G}}_{23}$, $\hat {{ Q}}_{23}$, $\hat {{E}}_{23}$, $\hat {{ H}}_3$, $\hat {{ H}}_4$, $\hat {{L}}_3$, $\hat {{ L}}_4$对应，此时只需将其下标2和3逐次用1和2替换即可.
上述推导建立了广义载荷作用于层状介质体系内部时的一般性解答，对于$H_F= H_C$, $H_C = 0$或$H_F = 0$等工况，属于以上解答的特例.计算时只需根据相应工况，将三层块中不存在层块的关系矩阵$\hat {{F}}$, $\hat {{ G}}$, $\hat {{ Q}}$, $\hat {E}$分别用I, 0, 0, I代替，而后采用上述方法，可实现层状体系特殊受荷时的解答.
当层状体系内部分别作用有恒定总强度为Q0的点热源和半径为的圆形热源时，有
$$\left.\begin{array}{lll} {P}_\varGamma = \left[{0, 0, 0, {Q_0 } /{2\pi s^2}} \right]^{\rm T}\\{ P}_\varGamma = \left[{0, 0, 0, {Q_0 {\rm J}_1 (\xi r_0 )}/{\pi r_0 \xi s^2}} \right]^{\rm T}\end{array}\right\}$$ (23) 若内部作用的热源为衰变热源时，假定其初始强度为Q0，其强度随时间衰减$Q = Q_0 {\rm e}^{ - \gamma t}$，其中$\gamma ={\ln 2} /{t'}$，t'为热源的半衰期，此时有
$$\left.\begin{array}{lll} \label{eq47} { P}_\varGamma = \left[{0, 0, 0, {Q_0 } /{2\pi s(s +\gamma )}} \right]^{\rm T}\\{ P}_\varGamma = \left[{0, 0, 0, {Q_0 {\rm J}_1 (\xi r_0 )}/{\pi r_0 \xi s(s + \gamma )}} \right]^{\rm T}\end{array}\right\}$$ (24) 结合式 (21)~式 (24)，可得到所求问题在积分变换域内的解答.
4. 数值计算与分析
第3节得到了问题在变换域内的解答，对其进行相应的积分逆变换，可求得其物理域内的最终解.本文Laplace数值逆变换采用FT[35]法
$$\begin{array}{lll} g\left( {t, M} \right) = \dfrac{\phi }{M}\Big\{ {\dfrac{1}{2}\tilde{g}\left( \phi \right)\exp \left( {\phi t} \right) + } \\\quad {\sum\limits_{\chi = 1}^{M - 1} {{\rm Re}\left[{\exp\left( {ts\left( {\zeta _\chi } \right)} \right)\tilde {g}\left({s\left( {\zeta _\chi } \right)} \right)\left( {1 + {\rm i}\vartheta \left( {\zeta _\chi } \right)} \right)} \right]} }\Big\} \end{array}$$ (25) 式中
$$\left.\begin{array}{ll}\vartheta \left( \zeta \right) = \zeta +\left( {\zeta \cot \zeta - 1} \right)\cot \zeta\\s\left( \zeta \right) = \phi \zeta \left( {{\rm i} + \cot \zeta }\right), - \pi < \zeta < \pi\\\phi = {2M} /{ {5t}}\\\zeta _\chi = {\chi \pi } /M\end{array}\right\}$$ (26) 其中$g(t)$为$\tilde{g}(s)$在物理域内相对应的函数.
由式 (25) 可知，$g(t)$的求解表达式中仅含有一个自由的参数M，即累加求和项数.为控制$g(t)$求解计算中的舍入误差，参考文献[36]的建议，对精度要求需定义M，即所需计算精度的有效位数.
因此，FT法的求解可概述为：确定积分变换式，并指定计算参数t和M的值.首先设定与计算精度有关的参数M(参考文献[30, 36]的建议，本文M值取为10)，随后根据式 (25) 和式 (26) 计算$g(t, M)$的值.参数M对计算精度及效率的影响分析可参考文献[30, 36].
本文Hankel数值逆变换借鉴文献[37]的方法实现.考虑到Bessel函数是震荡衰减函数，Gauss积分点或正或负，为保证数值计算的稳定性，计算时采用零点分段的方法.首先，选取Bessel函数两相邻零点为一积分区段，将原Hankel逆变换的半无限积分区间划分为多个积分区段；然后，每个区段采用Gauss-Legendre法进行数值积分计算；最后，将每个区段的积分结果进行叠加可得到最终的积分值.此法在理论上可表达为：
对于函数$\bar {g}^m\left( \varphi \right)$的第i个积分区间$\left[{\xi _i, \xi _{i + 1} } \right]$，若令$\varphi = {\left( {\xi _i +\xi _{i + 1} } \right)} /2 + {\left( {\xi _{i + 1} - \xi _i }\right)\varpi } / 2$，则函数可转化为积分区间为$\left[{-1, 1}\right]$的积分
$$\int_{\xi _i }^{\xi _{i + 1} } {\bar {g}^m\left(\varphi \right){\rm d}\varphi } = \frac{\xi _{i + 1} - \xi _i}{2}\int_{ - 1}^1 {\bar {g}^m\left( \varpi \right)} {\rm d}\varpi$$ (27) 则积分区间$\left[{\xi _i, \xi _{i + 1} }\right]$的Gauss-Legendre数值积分式如下
$$\int_{\xi _i }^{\xi _{i +1} } {\bar {g}^m\left( \varphi \right){\rm d}\varphi } = \frac{\xi_{i + 1} - \xi _i }{2}\sum\limits_j {\omega _j } \bar {g}^m\left({\varphi _j } \right)$$ (28) 将N个子积分区段的结果进行叠加，可得整个区间的积分值
$$\begin{array}{ll} g\left( r \right) \approx H^{ - 1}\left\{ {\bar {g}^m\left( \xi\right), \left[{\xi _0, \xi _N } \right]} \right\} \approx \\\quad \sum\limits_{i = 0}^{N - 1} {\left( {\frac{\xi _{i + 1} -\xi _i }{2}\sum\limits_j {\omega _j \bar {g}^m\left( {\varphi_{ij} } \right){\rm J}_m \left( {\varphi _{ij} r} \right)\psi_{ij} } } \right)} \end{array}$$ (29) 式中，$H^{ - 1}\left\{ {\bar {g}^m\left( \xi \right), \left[{\xi_0, \xi _N } \right]} \right\}$表示对$\bar {g}^m\left( \xi\right)$进行积分限为$\left[{\xi _0, \xi _N }\right]$的Hankel逆变换；$g\left( r \right)$为$\bar {g}^m\left( \xi\right)$对应物理域内的值；其中$\xi _i \left( {i = 0\sim N}\right)$为m阶Bessel函数的第i个零点对应的区间限值，且$\xi _0 =0$，N取足够大以近似表示无穷大积分区间，$\varphi _{ij} = {\left({\xi _{i + 1} + \xi _i } \right)}/2 + {\left( {\xi _{i + 1} - \xi_i } \right)\varpi _j } /2$；$\varpi _j$和$\omega _j$分别为Gauss-Legendre积分法中第j个Gauss积分点和Gauss系数.参数对逆变换的影响分析详见文献[36].
4.1 验证
为验证本文方法及程序的正确性，本文对点热源作用下饱和多孔介质热固结问题进行计算，并将结果与文献[9]进行对比，如图 4所示.点热源强度为Q0且保持不变，热源埋置于介质的深部.无量纲因子定义见图 4，由图可知，超静孔压随着时间的发展先逐渐增大，达到一峰值后而逐渐消散.另外，对比图中结果可知，本文结果与文献[9]的成果在各时刻均显示出较高的吻合度，从而证明本文方法在求解该类问题中的适用性.
![]()
图 4 点热源作用时超静孔压随时间变化Figure 4. Variation of excess pore pressure with time due to a point heat source4.2 热源衰变周期的影响分析
图 5(a)和图 5(b)展示了单层介质热固结问题中热源衰变周期对温度和超静孔压的影响.计算点位于z轴，与点热源间距为b.泊松比$\mu =0.3$保持不变，参数设置如图 5(a)所示，无量纲因子与4.1节设置相同.由图 5(a)可知，衰变周期对温度峰值及达到峰值所需时间影响明显；衰变周期越大，温度峰值越大，且达到峰值所需时间越长；当热源为无衰变热源时 ($t_0 = \infty )$，温度随时间的发展而增大，最终趋于稳定.由图 5(b)可知，衰变周期越长，超静孔压峰值越大且达到峰值所需时间越长；衰变周期较短时会有负超静孔压出现，这与热源的迅速冷却有关.
![]()
图 5 衰变周期对热固结效应的影响Figure 5. Effect of the heat source's half-life on thermal consolidation4.3 热源埋深的影响分析
本节通过两组算例研究热源埋深对多孔饱和介质热固结响应的影响.图 6(a)和图 6(b)分别展示了点热源埋深h对超静孔压及竖向位移的影响，假定介质厚度为H，且$H = 10h_0$，h0为常数，泊松比$\mu =0.25$，无量纲因子设置仍与4.2节相同.热源埋深考虑三种情况，分别为$h= h_0$，$h=2h_0$和$h=5h_0$.由图 6(a)可知，随着计算点与热源点间距离的增大，超静孔压峰值逐渐减小，且达到峰值时间逐渐延长.
![]()
图 6 热源埋深对超静孔压及竖向位移的影响Figure 6. Effect of the heat source's buried depth on the variation of excess pore pressure and vertical displacement从图 6(b)可以发现，热源埋置较深时，表面竖向位移峰值较大，但热源埋置深度对竖向位移的稳定值影响不大.
接下来分析深埋和浅埋热源作用下竖向位移随时间的变化规律，具体如图 7.热源为圆形热源，其半径和总强度分别为r0和Q0，深埋和浅埋深度分别为$1~000r_0$和$10r_0$，无量纲因子见图 7(a).由图知，深埋热源上下两侧等距离处的位移变化曲线呈对称分布，并分别向两侧发生膨胀变形.当$\tau <1$时，浅埋热源两侧位移变化基本呈对称分布，这是由于在较短时间内，热量仅传导至周围较小的区域；随着时间的发展，热量继续向更深处传递，由于为浅埋热源，热量传递至介质表面而使其上部介质膨胀量不再产生变化，热源下部由于为半空间而导致竖向位移继续变化，下侧介质不断膨胀使热源及计算点位移逐渐向上发展，从而引起图 7(b)中现象的产生.
![]()
图 7 不同埋深条件下竖向位移随时间的变化曲线Figure 7. Evolutions of vertical displacement against time with different conditions of buried depths4.4 介质成层性的影响分析
实际工程中的天然介质往往呈层状分布，本节以四层饱和多孔半空间为例分析层状特性对热固结响应的影响.如图 8(a)所示，1~3层为有限深度层，第4层为半空间，选5种不同工况进行分析，如表 1所示.
![]()
图 8 不同时刻多层介质超静孔压沿深度分布Figure 8. Distribution of excess pore pressure along depth of multilayered medium at different time表 1 多层饱和半空间计算参数Table 1. Parameters of the multilayered half-space
其他参数有如下关系：$\mu _i = 0.25(i = 1, 2, 3, 4)$，$\Delta h_1 :\Delta h_2 :\Delta h_3 :h = 5:3:2:10$，$G_1 :G_2:G_3 :G_4 = 8:5:2:1$，h为点热源的埋置深度，其强度为Q0，其他无量纲参数设置见图 8(d).
图 8展示了不同时刻z轴上不同深度计算点的超静孔压分布.由图可知，在较早时刻 ($\tau = 0.02, 0.05)$，超静孔压主要产生于热源周围区域，且距离热源越近，其值越大；由于在不同介质层分界面处，介质参数出现突变，导致分布曲线在分界面处出现明显的折点；随着时间的延长，超静孔压的区域范围逐渐扩大，且分界面处的曲线逐渐趋于平滑.
5. 结论
本文在Laplace-Hankel变换域内推导出饱和多孔介质热固结问题的常微分控制方程，并借助于层状体系任意深度作用载荷时的精细积分法和数值逆变换技术，得到了衰变热源作用下层状饱和多孔介质热固结问题的解答.通过与已有文献对比，验证了本文方法在求解热固结问题时的正确性.最后，通过算例分析了热源衰变周期、热源埋深及介质成层性对热固结效应的影响，结果表明：
(1) 热源衰变周期对温度和超静孔压的峰值、达到峰值的时间均有明显影响；衰变周期越大，二者峰值均越大，且达到峰值所需时间越长.
(2) 热源埋深对超静孔压及竖向位移变化影响显著；深埋热源作用时热源两侧竖向位移呈对称分布，浅埋热源两侧位移曲线则无此现象.
(3) 本文方法能有效求解多层介质的热固结问题，介质的分层特性对热固结效应影响明显.
-
![]()
图 4 点热源作用时超静孔压随时间变化
Figure 4. Variation of excess pore pressure with time due to a point heat source
![]()
图 5 衰变周期对热固结效应的影响
Figure 5. Effect of the heat source's half-life on thermal consolidation
![]()
图 6 热源埋深对超静孔压及竖向位移的影响
Figure 6. Effect of the heat source's buried depth on the variation of excess pore pressure and vertical displacement
![]()
图 7 不同埋深条件下竖向位移随时间的变化曲线
Figure 7. Evolutions of vertical displacement against time with different conditions of buried depths
![]()
图 8 不同时刻多层介质超静孔压沿深度分布
Figure 8. Distribution of excess pore pressure along depth of multilayered medium at different time
-
[1] Delage P, Sultan N, Cui YJ. On the thermal consolidation of Boom clay. Canadian Geotechnical Journal, 2000, 37(2):343-354 doi: 10.1139/t99-105
[2] 王铁行, 李宁, 谢定义.土体水热力耦合问题研究意义、现状及建议.岩土力学, 2005, 26(3):488-493 http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX20050300U.htm Wang Tiehang, Li Ning, Xie Dingyi. Necessity and means in research on soil coupled heatmoisture-stress issues. Rock and Soil Mechanics, 2005, 26(3):488-493(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX20050300U.htm
[3] 蒋中明, Dashnor H.核废料贮存库围岩体热响应耦合场研究.岩土工程学报, 2006, 28(8):953-956 http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200608004.htm Jiang Zhongming, Dashnor H. Studies on coupled field of thermal response in rock mass of nuclear waste repository. Chinese Journal of Geotechnical Engineering, 2006, 28(8):953-956(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200608004.htm
[4] Bai B, Guo LJ, Han S. Pore pressure and consolidation of saturated silty clay induced by progressively heating/cooling. Mechanics of Materials, 2014, 75:84-94 doi: 10.1016/j.mechmat.2014.04.005
[5] Ai ZY, Wang LJ. Axisymmetric thermal consolidation of multilayered porous thermoelastic media due to a heat source. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(17):1912-1931 doi: 10.1002/nag.v39.17
[6] Biot MA. General theory of three-dimensional consolidation. Journal of Applied Physics, 1941, 12(2):155-164 doi: 10.1063/1.1712886
[7] Biot MA. Thermoelasticity and irreversible thermodynamics. Journal of Applied Physics, 1956, 27(3):240-253 doi: 10.1063/1.1722351
[8] Booker JR, Savvidou C. Consolidation around a spherical heat source. International Journal of Solids and Structures, 1984, 20(11-12):1079-1090 doi: 10.1016/0020-7683(84)90091-X
[9] Savvidou C, Booker JR. Consolidation around a heat source buried deep in a porous thermoelastic medium with anisotropic flow properties. International Journal for Numerical and Analytical Methods in Geomechanics, 1989, 13(1):75-90 doi: 10.1002/(ISSN)1096-9853
[10] Mctigue DF. Thermoelastic response of fluid-saturated porous rock. Journal of Geophysical Research Atmospheres, 1986, 91(B9):9533-9542 doi: 10.1029/JB091iB09p09533
[11] Bai M, Abousleiman Y. Thermoporoelastic coupling with application to consolidation. International Journal for Numerical and Analytical Methods in Geomechanics, 1997, 21(2):121-132 doi: 10.1002/(ISSN)1096-9853
[12] 白冰.岩土介质非稳态热固结耦合问题的热源函数法.力学学报, 2004, 36(4):427-434 http://lxxb.cstam.org.cn/CN/abstract/abstract141271.shtml Bai Bing. Heat source function method for coupling analyses of thermal consolidation in saturated soil. Acta Mechanica Sinica, 2004, 36(4):427-434(in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract141271.shtml
[13] 郑荣跃, 刘干斌, 梧松.半空间饱和土内置点载荷作用下的热弹性波动.力学学报, 2008, 40(3):413-420 http://lxxb.cstam.org.cn/CN/abstract/abstract141668.shtml Zhen Rongyue, Liu Ganbin, Wu Song. Coupling thermo-hydro-mechanical dynamic response of saturated soil subjected to internal excitation. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(3):413-420(in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract141668.shtml
[14] 吴瑞潜, 谢康和, 程永锋.变载荷下饱和土一维热固结解析理论.浙江大学学报 (工学版), 2009, 43(8):1532-1537 http://www.cnki.com.cn/Article/CJFDTOTAL-ZDZC200908034.htm Wu Ruiqian, Xie Kanghe, Cheng Yongfeng. Analytical theory for one-dimensional thermal consolidation of saturated soil under time-dependent loading. Journal of Zhejiang University (Engineering Science), 2009, 43(8):1532-1537(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-ZDZC200908034.htm
[15] Lu JCC, Lin F. Thermal consolidation of a poroelastic full space subjected to a decaying point heat source//Proceedings of the 2nd International ISCM Symposium and the 12nd International EPMESC Conference, 2010:407-412
[16] Selvadurai APS, Suvorov AP. Thermo-poromechanics of a fluidfilled cavity in a fluid-saturated geomaterial//Proceedings of the royal society a mathematical physical and engineering sciences, 2014, 470(2163):20130634
[17] Giraud A, Homand F, Rousset G. Thermoelastic and thermoplastic response of a double-layer porous space containing a decaying heat source. International Journal for Numerical and Analytical Methods in Geomechanics, 1998, 22(2):133-149 doi: 10.1002/(ISSN)1096-9853
[18] 白冰.变温度载荷作用下半无限成层饱和介质的热固结分析.应用数学和力学, 2006, 27(11):1341-1348 http://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200611010.htm Bai Bing. Thermal consolidation of layered porous half-space to variable thermal loading. Applied Mathematics and Mechanics, 2006, 27(11):1341-1348(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YYSX200611010.htm
[19] Yue ZQ, Yin JH. Backward transfer-matrix method for elastic analysis of layered solids with imperfect bonding. Journal of Elasticity, 1998, 50(2):109-128 doi: 10.1023/A:1007421014760
[20] 艾智勇, 吴超.渗透各向异性可压缩地基固结的平面应变分析.力学学报, 2009, 41(5):801-807 http://lxxb.cstam.org.cn/CN/abstract/abstract141853.shtml Ai Zhiyong, Wu Chao. Analysis on plane strain consolidation of a multi-layered soil with anisotropic permeability and compressbility constituents. Chinese Journal of Theoretical and Applied Mechanics, 2009, 41(5):801-807(in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract141853.shtml
[21] 赵宇昕, 陈少林.关于传递矩阵法分析饱和成层介质响应问题的讨论.力学学报, 2016, 48(5):1145-1158 http://lxxb.cstam.org.cn/CN/abstract/abstract146012.shtml Zhao Yunxin, Chen Shaolin. Discussion on the matrix propagator method to analyze the response of saturated layered media. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5):1145-1158(in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract146012.shtml
[22] Booker JR, Small JC. Finite layer analysis of consolidation. Ⅰ. International Journal for Numerical and Analytical Methods in Geomechanics, 1982, 6(2):151-171 doi: 10.1002/(ISSN)1096-9853
[23] 宰金珉, 梅国雄.有限层法求解三维比奥固结问题.岩土工程学报, 2002, 24(1):31-33 http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200201006.htm Zai Jinmin, Mei Guoxiong. Finite layer analysis of three dimensional Biot consolidation. Chinese Journal of Geotechnical Engineering, 2002, 24(1):31-33(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200201006.htm
[24] 钟阳, 耿立涛.多层弹性平面问题解的精确刚度矩阵法.岩土力学, 2008, 29(10):2829-2832 http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200810046.htm Zhong Yang, Geng Litao. Explicit solution of multiplayer elastic plane by exact stiffness matrix method. Rock and Soil Mechanics, 2008, 29(10):2829-2832(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200810046.htm
[25] 艾智勇, 曹国军, 成怡冲.平面应变Biot固结的解析层元.力学学报, 2012, 44(2):401-407 http://lxxb.cstam.org.cn/CN/abstract/abstract143147.shtml Ai Zhiyong, Cao Guojun, Cheng Yichong. Analytical layer-element of plane strain Biot's consolidation. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(2):401-407(in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract143147.shtml
[26] 艾智勇, 王路君, 曾凯.稳定温度场下层状路面体系的解析层元解.同济大学学报 (自然科学版), 2014, 42(11):1665-1669 http://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ201411006.htm Ai Zhiyong, Wang Lujun, Zeng Kai. Analytical layer-element solution for layered pavement in stable temperature field. Journal of Tongji University (Natural Science), 2014, 42(11):1665-1669(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ201411006.htm
[27] 钟万勰.结构动力方程的精细时程积分法.大连理工大学学报, 1994, 34(2):131-136 http://www.cnki.com.cn/Article/CJFDTOTAL-DLLG199402003.htm Zhong, Wanxie. On precise time-integration method for structural dynamics. Journal of Dalian University of Technology, 1994, 34(2):131-136(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-DLLG199402003.htm
[28] 钟万勰.弹性力学求解新体系.大连:大连理工大学出版社, 1995 Zhong Wanxie. A New Systematic Methodology for Theory of Elasticity. Dalian:Dalian University of Technology Press, 1995(in Chinese)
[29] 韩泽军, 林皋, 李建波.二维层状地基格林函数的求解.土木工程学报, 2015, 48(10):99-107 http://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201510014.htm Han Zejun, Lin Gao, Li Jianbo. The solution of Green's functions for two-dimensional layered ground. China Civil Engineering Journal, 2015, 48(10):99-107(in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201510014.htm
[30] Ai ZY, Cheng YC. Extended precise integration method for consolidation of transversely isotropic poroelastic layered media. Computers & Mathematics with Applications, 2014, 68(12):1806-1818 https://www.researchgate.net/publication/268695249_Extended_precise_integration_method_for_consolidation_of_transversely_isotropic_poroelastic_layered_media
[31] Wang LJ, Ai ZY. Plane strain and three-dimensional analyses for thermo-mechanical behavior of multilayered transversely isotropic materials. International Journal of Mechanical Sciences, 2015, 103:199-211 doi: 10.1016/j.ijmecsci.2015.09.006
[32] Talbot A. The accurate numerical inversion of Laplace transforms. Journal of Institute of Mathematics and Its Application, 1979, 23(1):97-120 doi: 10.1093/imamat/23.1.97
[33] Sneddon IN. The Use of Integral Transform. New York:McGrawHill, 1972
[34] Zhong WX, Lin JH, Gao Q. The precise computation for wave propagation in stratified materials. International Journal for Numerical Methods in Engineering, 2004, 60(1):11-25 doi: 10.1002/(ISSN)1097-0207
[35] Bailey DH, Swarztrauber PN. A fast method for the numerical evaluation of continuous Fourier and Laplace transforms. SIAM Journal on Scientific Computing, 1994, 15(5):1105-1110 doi: 10.1137/0915067
[36] Abate J, Valko PP. Multi-precision Laplace transform inversion. International Journal for Numerical Methods in Engineering, 2004, 60(5):979-993 doi: 10.1002/(ISSN)1097-0207
[37] Ai ZY, Yue ZQ, Tham LG, et al. Extended Sneddon and Muki solutions for multilayered elastic materials. International Journal of Engineering Science, 2002, 40(13):1453-1483 doi: 10.1016/S0020-7225(02)00022-8
下载:
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
