0.   引言

多孔介质的瞬态响应问题是动力学的重要分支，在地震工程、海洋工程和岩土工程等领域受到广泛的关注，具有重要的理论价值和实际应用。

自Biot^[1-2]提出饱和多孔介质波动理论以来，众多学者通过求解Biot波动方程，获得一维、半空间以及层状饱和多孔介质在动力荷载下时域或频域内的基本解或者动力响应结果，研究波的传播和弥散机理以及多孔介质的振动和应力孔压等响应情况。在基本解方面，Zheng等^[3]将波动问题简化为由Helmholtz控制方程的4个标量势确定，给出柱坐标系下多孔弹性半空间的通解。Sahebkar等^[4]采用标量势函数和积分变换方法给出点荷载和流体源作用下横观各向同性多孔弹性介质基本解。Wang等^[5]采用傅里叶变换给出频域内三维各向异性饱和半空间的格林函数。对于成层饱和多孔介质响应问题，Lu等^[6]通过Hankel变换和Fourier变换提出点荷载或流体点源作用下分层饱和半空间的透射和反射矩阵法（TRM）。Ba等^[7]分析了多层横观各向同性多孔弹性半空间的三维动态响应。对于移动荷载下多孔介质响应问题，Cai等^[8]考虑路面不平顺，采用傅里叶变换分析了移动卡车荷载下饱和多孔弹性半空间的动力响应。Li等^[9]将轨道视为Timoshenko梁，研究了移动荷载下成层横观各向同性饱和路基的动力响应。在一维瞬态响应方面，Schanz等^[10]采用Laplace变换给出无黏性耦合情况下的解析解。Shan等^[11]采用分离变量法和傅里叶正余弦变换法给出单层饱和多孔介质在三类典型边界条件下的级数解以及半无限问题的积分解。凌道盛等^[12]给出时域内顶部和底部均为透水边界的单层饱和多孔介质半解析解。Schanz等^[13]采用Kelvin-Voigt模型替换Biot的线弹性固体骨架模型，获得Laplace域内一维黏弹性饱和多孔介质的解析解。Pisanò等^[14]采用改进的分步Taylor- Galerkin有限元法获得黏弹性数值解。Shan等^[15]采用傅里叶变换法获得时域内饱和黏弹性多孔介质一维瞬态响应解。Shan等^[16]采用特征函数法给出上覆饱和流体时的饱和多孔介质一维瞬态响应半解析解。

可以看出，目前关于Biot方程的解析研究主要通过积分变换的方式求解，已有诸多解答。然而，由于多孔介质的惯性、黏滞和机械耦合作用比较复杂，多数问题都进行了简化，多数解答均是通过时间和空间离散的数值方法获得，频域的解析解答需经过逆变换转为时域，而逆变换过程如Laplace逆变换往往存在数值不稳定问题^[17]。多数精确解和解析解集中在一维问题中，即使是一维问题，也仅能给出满足一定边界条件和材料参数要求的解答，需要寻求适应更多边界条件的高效解答。此外，上述研究中边界透水条件均沿用传统的太沙基类边界，采用孔压或液相位移描述边界透水情况，然而，这并不能反映边界实际复杂的透水情况，如分析域边界处存在上覆砂垫层等工况，边界将变为半透水情况。众所周知，边界情况对响应结果具有较大的影响，如何针对复杂的实际工况，对边界进行修正是动力问题另一个需要进一步研究的工作。

对于一维问题的研究，一是可以验证相关数值解法的正确性和精度，更重要的是可以深入理解多孔介质的动力特性。本文以单层饱和多孔介质为研究对象，对其在任意竖向荷载作用下的一维瞬态响应问题展开研究，提出具有普遍意义的一般性透水边界条件，考虑多孔介质各相间的惯性、黏滞和机械耦合作用，结合特征函数法和常数变易法，在时域内直接求解，并系统分析了不同透水条件下单层饱和多孔介质的瞬态响应特点。

1.   控制方程和定解模型

本文以一般性透水边界条件下单层饱和多孔介质瞬态响应问题为例，建立如图 1所示模型，H为单层饱和介质厚度。模型顶部作用任意竖向荷载，底部为固定位移，边界处为任意的透水条件。

图  1  不同透水边界条件下单层饱和多孔介质模型

Figure  1.  1D model for saturated porous media under various permeable boundary conditions

下载: 全尺寸图片 幻灯片

1.1   基本方程

采用Simon等^[18]的表述方法，Biot饱和多孔介质整体平衡方程和流体平衡方程基本方程为

渗流连续方程和固体骨架本构方程为

式中：u_z为固体骨架的绝对位移；w_z为流体的相对位移，${w_z} = n({\bar w_z} - {u_z})$，${\bar w_z}$为流体的绝对位移，n为材料孔隙率；${\sigma _z}$为介质总应力，以拉为正；P为孔隙流体压力，以压为正；$\rho$为饱和多孔介质质量密度，$\rho = {\text{(1}} - n{\text{)}}{\rho _\text{s }} + n{\rho _\text{f }}$，${\rho _\text{f }}$为流体质量密度，${\rho _{\text{s}}}$为固体颗粒的质量密度；$\bar \rho = {\rho _\text{f }}/n$；下标t为时间；${k_\text{f }}$为动力渗透系数，与Darcy定律的渗透系数和流体重度${\gamma _\text{f }}$间满足${k_{\text{f}}} = k/{\gamma _\text{f }}$；$\lambda$，μ为固体骨架的Lame常数；$\alpha$，Q为与固相和液相压缩性相关的系数，且$\alpha = 1 -$ ${K_\text{b}}/{K_\text{s}}$，$1/Q = (\alpha - n)/{K_\text{s}} + n/{K_\text{f}}$，${K_\text{b}} = \lambda + 2\mu /3$，${K_\text{s}}$，${K_\text{f }}$和${K_\text{b}}$分别为固体颗粒，液相和固体骨架的体积模量。当${K_\text{s}} \to \infty$时，固体颗粒不可压缩；当${K_\text{f}} \to \infty$时，流体不可压缩。

为方便起见，将式（1）~（4）量纲归一化，并写出矩阵形式：

式中：$\boldsymbol{u} = \left( {\begin{array}{*{20}{c}} u \\ w \end{array}} \right)$，$\boldsymbol{K} = \left( {\begin{array}{*{20}{c}} 1&{\alpha \kappa } \\ {\alpha \kappa }&\kappa \end{array}} \right)$，$\boldsymbol{M} = \left( {\begin{array}{*{20}{c}} 1&\beta \\ \beta &\gamma \end{array}} \right)$，$\boldsymbol{C} = \left( {\begin{array}{*{20}{c}} 0&0 \\ 0&{2\eta } \end{array}} \right)$，$\tau = \frac{{\sqrt {m/\rho } }}{H}t$。其中$\zeta = \frac{z}{H}$，$u = \frac{{{u_z}}}{H}$，$w = \frac{{{w_z}}}{H}$，$\beta = \frac{{{\rho _\text{f }}}}{\rho }$，$\kappa = \frac{Q}{m}$，$\gamma = \frac{{\overline \rho }}{\rho }$，$\eta = \frac{H}{{2{k_\text{f }}\sqrt {m\rho } }}$，$m = \lambda + 2\mu + {\alpha ^2}Q$。

1.2   定解条件

在固结问题中，沿用太沙基边界条件，用孔压表示透水和不透水，当孔压采用第三类边界条件时，则可以表示半透水或任意透水条件，

式(6b)中${b_0}$和${c_0}$不同时为0。当${a_0} = 0$或${b}_{0}=0，{c}_{0}\ne 0$时，式（6）退化为太沙基不透水条件，当${a_0} \to \infty$或${b}_{0}\ne 0，{c}_{0}=0$时，式（6）退化为太沙基透水条件。其他条件，则式（6）可以不同程度的半透水条件。

在动力问题中，对于透水条件同样采用太沙基边界，即孔压$p = 0$表示，而不透水边界的表示方法则有两种表述方式，一类是同样沿用太沙基边界条件：$\partial p/$ $\partial z = 0$^[19]，另一类则是采用液相位移$w = 0$表示^[11]。两者分别从孔压和位移的角度描述边界的不透水情况，在一定程度上都是合理的。然而，两者并不完全等价，研究表明，在动力问题中，如果采用$\partial p/\partial z = 0$表示不透水边界，往往会使解答不适定，即解不唯一，如果想要获得唯一解，需要在边界处补充$w = f(\tau )$的条件^[19-20]。因此，在动力问题中，较多采用$w = 0$的边界条件表示不透水。

以往动力研究中较少分析半透水形式，实际上，边界条件尤其是透水条件，对动力响应结果同样有较大影响，极端的透水条件并不能准确反映结构和路基的振动情况，因此有必要引入半透水边界。

本文采用以下边界条件：

式中：两个参数${a_i}$，${b_i}$为不同时为0的已知实数，$i = 1,2$。式（7b），（7d）反映边界的透水情况，当${a_{1or2}} = 0$时，为不透水边界，当${b_{1or2}} = 0$时，为透水边界，均不为0时，$- {a_i}/{b_i}$的大小表示不同程度的半透水边界，反映模型边界处附近渗透系数和多孔介质层厚度的影响。

量纲归一化的初始条件表示为

式中：${\boldsymbol{g }}_{1}(\zeta )={({g}_{11}(\zeta )，{g}_{12}(\zeta ))}^{\text{T}}$；${\boldsymbol{g }_2}(\zeta ) =$$({g}_{21}(\zeta )，$ ${g_{22}}(\zeta ){)^\text{T}}$。${g_{11}}(\zeta )$，${g_{12}}(\zeta )$，${g_{21}}(\zeta )$，${g_{22}}(\zeta )$均为任意函数。

2.   半解析解推导

2.1   边界齐次化

将式（5）的解分为两部分：

式中：${\boldsymbol{u }^\text{d}} = {({u^\text{d}},{w^\text{d}})^\text{T}}$；${\boldsymbol{u }^\text{s}} = {({u^\text{s}},{w^\text{s}})^\text{T}}$。

${\boldsymbol{u }^\text{s}}$采用如下表达式满足边界条件，${u^\text{s}}$，${w^\text{s}}$可分别采用形式为

式（10），（11）两两组合共有4种可尝试齐次化方式。

将式（10），（11）代入边界条件（7）可得

式中：$i=1，\text{ }2，\text{ }3，\text{ }4$。共可获得4个系数矩阵，首先判断系数矩阵行列式是否为0，不为0时，表明可采用这种形式进行齐次化，之后等式两边左乘系数矩阵的逆矩阵，即可给出$a$，$b$，$c$，$d$的值，继而给出${\boldsymbol{u }^\text{s}}$的值。

将式（9）~（11）代入式（5），（7），（8），可得${\boldsymbol{u }^\text{d}}$满足的定解条件：

2.2   特征值和特征函数

首先求解式（13）对应的无阻尼特征值问题。

假定方程具有如下分离变量形式的解：

式中：$\boldsymbol{U } = {(\varphi (\zeta ),\psi (\zeta ))^\text{T}}$；ω≥0。

代入无阻尼特征方程和边界条件（14）得

当$\omega > 0$时，式（17）具有如下形式解：

式中：$\chi$为待定常数。

将式（19）代入式（17）得

式（20）有非零解，系数矩阵的行列式等于0，由此可得

式中：${h_3} = \kappa - {\alpha ^2}{\kappa ^2}$；${h_2} = \kappa + \gamma - 2\alpha \kappa \beta$；${h_1} = \gamma - {\beta ^2}$。

则可获得方程（21）相应的4个特征值：

式中：

对应的特征矢量${\boldsymbol{J }_i}$（$i = 1,2,3,4$）为

则$U$的表达式可写为

将式（25）代入边界条件（18）并整理可得

式中：

$\begin{array}{l}{c}_{41}={b}_{2}{J}_{12}\mathrm{cos}(\omega {\chi }_{1})+{a}_{2}\kappa \omega {\chi }_{1}{J}_{12}\mathrm{sin}(\omega {\chi }_{1})+{a}_{2}\alpha \kappa \omega {\chi }_{1}{J}_{11}\mathrm{sin}(\omega {\chi }_{1})；\\ {c}_{42}={b}_{2}{J}_{12}\mathrm{sin}(\omega {\chi }_{1})-{a}_{2}\kappa \omega {\chi }_{1}{J}_{12}\mathrm{cos}(\omega {\chi }_{1})-{a}_{2}\alpha \kappa \omega {\chi }_{1}{J}_{11}\mathrm{cos}(\omega {\chi }_{1})；\\ {c}_{43}={b}_{2}{J}_{22}\mathrm{cos}(\omega {\chi }_{2})+{a}_{2}\kappa \omega {\chi }_{2}{J}_{22}\mathrm{sin}(\omega {\chi }_{2})+{a}_{2}\alpha \kappa \omega {\chi }_{2}{J}_{21}\mathrm{sin}(\omega {\chi }_{2})；\\ {c}_{44}={b}_{2}{J}_{22}\mathrm{sin}(\omega {\chi }_{2})-{a}_{2}\kappa \omega {\chi }_{2}{J}_{22}\mathrm{cos}(\omega {\chi }_{2})-{a}_{2}\alpha \kappa \omega {\chi }_{2}{J}_{21}\mathrm{cos}(\omega {\chi }_{2})。\end{array}$式（26）若有非零解，则系数矩阵行列式为0，由此可得到关于特征值$\omega$的方程：

式（27）为超越方程，由此可获得一系列特征值$\omega$，由小到大排列，并记为${\omega _k}$（，），其对应的特征矢量为${\boldsymbol{R }_k} = {({A_{11}},{A_{12}},{A_{13}},{A_{14}})^{\text{T}}}$，${A_{1j}}$为式（26）系数矩阵第一行对应的代数余子式，$j=1，\text{ }2，\text{ }3，\text{ }4$。

当$\omega = 0$时，方程（17）有解：

式中：${d_j}$为待求常数，$j=1，\text{ }2，\text{ }3，\text{ }4$。

将式（28）代入边界条件（18）得

当系数行列式不为0时，有零解，此时$\boldsymbol{U }$为平凡解。

当系数行列式等于0时，则有非零解，此时需要边界条件各参数满足：

对应的特征向量为

式中：${B_{1j}}$为式（29）系数矩阵第一行对应的代数余子式，$j = 1,2,3,4$。

令式（26）中$\omega = 0$，并乘以${\boldsymbol{R }_0}$，使其等于${\boldsymbol{U }_0}$，可获得合适的${\boldsymbol{R }_0}$值，从而将$\omega = 0$的模态并入到${\boldsymbol{U }_k}$中。

即特征值${\omega _{\text{k}}}$对应的特征函数：

假定${\omega _{\text{k}}}$和${\omega _{\text{p}}}$为不同特征值，对应的特征函数为${\boldsymbol{U }_{\text{k}}}(\zeta )$和${\boldsymbol{U }_{\text{p}}}(\zeta )$，分别代入方程（17）并相减，并由$\boldsymbol{M }$的对称性可得

采用分部积分技术结合边界条件，可证明第一项为0。则当${\omega _{\text{k}}} \ne {\omega _{\text{p}}}$时，

即特征函数满足正交性。

2.3   动力问题的解

根据常数变异法，方程（13）具有如下形式解：

将式（35）代入式（13）得

式（36）左乘$boldsymbol{U}_{\text{p}}^{\text{T}}$，并对从0到1积分，利用特征函数的正交性可得：

式中：${c_{{\text{p}}k}} = \frac{1}{{{G_{\text{p}}}}}\int_0^1 {\boldsymbol{U}_{\text{p}}^{\text{T}}} \boldsymbol{C}{\boldsymbol{U}_k}{d_\zeta }$；${s_{\text{p}}}(\tau ) = - \frac{1}{{{G_{\text{p}}}}}\int_0^1 {\boldsymbol{U}_{\text{p}}^{\text{T}}} (\boldsymbol{Mu}_{,\tau \tau }^s +$ $\boldsymbol{Cu}_{,\tau }^s - \boldsymbol{Ku}_{,\zeta \zeta }^s){d_\zeta }$。

将（35）代入初始条件（15）中得

式（38）左乘$\boldsymbol{U}_{\text{p}}^{\text{T}}\boldsymbol{M}$，并对从0到1积分，利用特征函数的正交性，可得${T_k}$的初始条件：

求解式（37），（39）构成的初值问题，即可获得方程（13）${\boldsymbol{u}^d}$的解，继而加上${\boldsymbol{u}^s}$，就可获得原问题（5）的解。

2.4   常微分方程精细时程积分解

常微分方程组（37）中包含无限个常系数微分方程，特征值${\omega _k}$按由小到大得顺序排列。实际求解过程中，根据精度要求取前N次，式（37）改写为矩阵形式：

式中：$\hat{\boldsymbol{T}}=\left\{\begin{array}{c} T_1 \\ T_2 \\ \cdots \\ T_N \end{array}\right\}$，$\hat{\boldsymbol{C}}=\left\{\begin{array}{llll} c_{11} & c_{12} & \cdots & c_{1 N} \\ c_{21} & c_{22} & \cdots & c_{2 N} \\ \cdots & \cdots & \cdots & \cdots \\ c_{N 1} & c_{N 2} & \cdots & c_{N N} \end{array}\right\}$，$\hat{\boldsymbol{K}}=\left[\begin{array}{llll} \omega_1^2 & & & \\ & \omega_2^2 & & \\ & & \cdots & \\ & & & \omega_N^2 \end{array}\right]$，$\hat{\boldsymbol{f}}=\left\{\begin{array}{c} f_1 \\ f_2 \\ \cdots \\ f_N \end{array}\right\}$。

对应的初始条件仍为式（39），采用精细时程积分法^[21]给出上述初值问题的解答。

3.   算例

3.1   方法正确性验证

3个算例用来验证本文模型的正确性，3个算例分别代表不同的透水条件，顶部作用阶跃荷载，如表 1所示，计算中采用的参数如表 2所示。

表  1  验证算例中的边界参数

Table  1.  Boundary conditions of verification examples

算例来源	边界条件	备注
Shan等[11]	$\begin{gathered} {a_1} = {b_2} = 1 \hfill \\ {a_2} = {b_1} = 0 \hfill \\ q(\tau ) = - 1 \hfill \\ \end{gathered}$	顶部透水，底部不透水
凌道盛等[12]	$\begin{gathered} {a_1} = {a_2} = 1 \hfill \\ {b_1} = {b_2} = 0 \hfill \\ q(\tau ) = - 1 \hfill \\ \end{gathered}$	顶部和底部均透水
单振东[22]	$\begin{gathered} {a_1} = {a_2} = 0 \hfill \\ {b_1} = {b_2} = 1 \hfill \\ q(\tau ) = - 1 \hfill \\ \end{gathered}$	顶部和底部不透水

下载: 导出CSV 

| 显示表格

表  2  材料参数^[10-11]

Table  2.  Parameters of saturated porous media

参数	单位	值
Lame常数λ	MPa	144.7
剪切模量μ	MPa	98
固相体积模量Ks	GPa	11
动力渗透系数kf	m4/(N·s)	0.01
液相体积模量Kf	GPa	3.3
固相密度ρs	kg/m3	2700
液相密度ρf	kg/m3	1000
孔隙率n	—	0.48
厚度H	m	10

下载: 导出CSV 

| 显示表格

图 2，3给出了本文孔压和总应力的计算结果与文献中不同透水边界结果的对比，可以看出，本文结果和前人研究结果具有较好的一致性，从而证明了本文方法的正确性，与以往文献相比，本文的边界条件和求解方法更具有一般性。

图  2  流体孔压结果

Figure  2.  Comparison of fluid pressure results at ζ=0.5

下载: 全尺寸图片 幻灯片

图  3  总应力结果

Figure  3.  Comparison of total stress results at ζ=0.5

下载: 全尺寸图片 幻灯片

3.2   不同透水条件下响应

以底部固定不透水，顶部作用阶跃荷载，不同透水条件为例，研究饱和多孔介质在不同透水边界条件下的响应结果，各参数取值：${a}_{1}=-1，{b}_{2}=1，{a}_{2}=0，$ $q(\tau )=-1，$${b_1} = - 0.01, - 5, - 20, - 100, - 10000$。

从图 4~6可以看出，半透水条件下的瞬态响应结果介于透水和不透水之间，当${b_1}$为-0.01，-10000时，各分量响应解答已经与不透水条件和透水条件结果相近。透水条件下产生的位移要远大于不透水条件，在透水和不透水条件下，固相绝对位移和液相相对位移绝对值呈现先增大后趋于稳定的变化趋势，而在半透水条件下，固相绝对位移和液相相对位移绝对值均呈现先增加后减小的现象。

图  4  不同透水条件下孔压结果

Figure  4.  Dimensionless fluid pressure responses at ζ=0.5 under different boundary conditions

下载: 全尺寸图片 幻灯片

图  5  不同透水条件下总应力结果

Figure  5.  Dimensionless total stress responses at ζ=0.5 under different boundary conditions

下载: 全尺寸图片 幻灯片

图  6  不同透水条件下位移结果

Figure  6.  Dimensionless displacement responses at ζ=0.5 under different boundary conditions

下载: 全尺寸图片 幻灯片

3.3   半透水界面波动情况

图 7，8给出不同时刻，入射阶跃荷载在饱和多孔介质内部孔压的传播情况，计算工况和各参数取值均同3.2节，半透水情况时${b_1}$取值为-20。可以看出，入射波激发两种纵波，随着时间进行，两种波分开，快波和慢波分别约在0.9，4.0 s时到达底部，反射波和入射波产生相同相位的孔压，即若入射波激发正的孔压，反射波亦激发正的孔压。

图  7  不同深度处孔压

Figure  7.  Dimensionless fluid pressure responses under different depths

下载: 全尺寸图片 幻灯片

图  8  不同深度处孔压

Figure  8.  Dimensionless fluid pressure responses under different depths

下载: 全尺寸图片 幻灯片

4.   结论

本文针对单层饱和多孔介质一维瞬态响应问题，提出具有普遍意义的一般性透水边界条件，并在时域内直接给出一维问题瞬态响应半解析解，分析了不同透水条件下单层饱和多孔介质的瞬态响应特点。主要得出3点结论。

（1）本文建立的一般性透水边界条件可以较好地描述边界透水情况，可精确退化至现有边界条件。

（2）本文考虑惯性、黏滞和机械耦合作用，在时域内直接求解，给出的半解析解可较好克服Laplace逆变换存在的数值不稳定等问题，具有较好的通用性。

（3）半透水条件下，饱和多孔介质响应介于透水和不透水之间，在顶部半透水和底部固定位移的边界条件下，入射波和反射波产生相同相位的孔压。