基于相对湿度场和变模量本构的岩土膨胀分析模型

2023-10-12程树范曾亚武

工程科学与技术 2023年5期

程树范，曾亚武,2*，叶阳，高睿,2，李涵

(1.武汉大学土木建筑工程学院，湖北武汉 430072；2.武汉大学岩土与结构工程安全湖北省重点试验室，湖北武汉 430072；3.中国地质大学（武汉）工程学院，湖北武汉 430074)

中国红层地区膨胀土、膨胀岩的分布极为广泛，吸水后产生的膨胀变形对工程结构的安全十分不利，特别是对于边坡、路基和隧道等交通工程，由岩体膨胀导致的地质灾害较多[1-3]。

岩土材料的吸水膨胀是一个十分复杂的流固耦合过程，与矿物成分、应力状态及浸水条件等因素均有关[4-6]。目前，从宏观角度直接建立的膨胀模型都是基于一些数学上的假设[7]，如：经典的1维Huder-Amberg膨胀模型假设膨胀率与膨胀压力的对数呈线性关系[8]；Gysel[9]提出的膨胀模型则假设膨胀力与吸水率服从指数关系。这一类模型属于半经验模型，模型参数的确定一般需要数据拟合，缺乏理论基础。缪协兴等[10-11]借鉴了热弹性力学的分析方法，提出湿度应力理论，并给出了控制方程组。白冰等[12-13]给出了湿度应力场的理论证明，认为该理论在力学模型上是完备的；同时，也指出湿度应力场虽然是一种合理的理论模型，但需要根据实际情况进行必要的修正和完善。为此，在湿度应力理论框架下，研究人员开展了大量拓展。如：陈有亮等[14]推导了考虑膨胀应力和剪胀效应的深埋隧道弹塑性解；陈钒等[15]建立了考虑时间效应的膨胀模型；任松等[16]得到了考虑时间效应的隧道弹-膨胀解析模型。然而，为方便理论推导，现有研究多假设固体相的变形是弹性或者弹塑性的，这一假设忽略了土水或水岩作用下岩土材料的软化效应，是一个明显的不足。且从矿物成分的角度来看，膨胀岩、膨胀土多含有蒙脱石和伊利石等亲水性黏土矿物，一般也有较高的软化系数，遇水后变形模量会持续下降[17-19]。对于有侧限或围压的情况，完全软化状态的岩土体变形模量一般会趋于定值[20]。吴建勋等[21]研究了围岩软化对隧道围岩自稳时间的影响，认为围岩软化作用不能忽略。综上所述，在湿度应力的理论框架内，如果能够充分考虑变形模量随湿度的变化，就能够更好地为膨胀岩地区的地下工程提供指导。

本文以压力固结的红层黏土（后文简称固结土）和红层黏土岩（后文简称黏土岩）为研究对象，首先，开展受荷膨胀试验，通过侧限膨胀率，确定变形模量随含水率的变化规律；然后，基于归一化的相对湿度建立了变模量本构，给出模型参数的确定方法，并通过室内试验和相应数值模拟的对比，验证了本文膨胀分析模型的可靠性；最后，给出了湿度平均扩散系数的标定方法，以圆形洞室遇水为算例，进行理论计算和数值模拟，讨论了本文模拟方法的优越性。

1 湿度应力场的改进

1.1 湿度应力场的控制方程

湿度应力理论是对复杂水-岩土体物理膨胀问题的一种简化研究方法，借鉴了热弹性问题的经典理论。湿度应力场的控制方程包含湿度场的扩散方程，应力场的平衡微分方程、几何方程，以及材料的本构方程[10-12]。

岩土体中的水分迁移过程十分复杂，与土或岩石的渗透性、饱和度、孔隙结构乃至应力状态等因素都密切相关，湿度场的引入是对水分迁移过程的一种简化。湿度是含水率（或含水量）的表征。根据质量守恒，单位体积含水量的增加应与水分的净流入量相等；而根据达西定律，渗流速度与湿度梯度成正比，那么湿度场的控制方程实际是一个无源扩散方程，各向同性条件下可以表示为：

式中：W为湿度场，对应的场变量为含水率w，是坐标x、y、z和时间t的函数；κ为湿度扩散率，在各向同性条件下为常数。

应力场的平衡方程、几何方程为：

式中，εij为应变分量，σij为应力分量，ui和uj为位移分量，xi和xj为坐标分量，fi为体力分量。

湿度应力场与固体力学场的区别主要表现在本构模型方面，除荷载外，湿度的变化也会引起材料的体变，无约束状态时[10]：

式中：α为膨胀系数；δij为克罗内克记号，仅当i=j时取1，其余取0。

有约束状态时，湿度场的变化会产生附加应力，附加应力与应变间满足广义虎克定律，总的本构模型就可以描述为[10]：

式中， Θ为应力张量的第1不变量，E、 ν为弹性模量和泊松比。

湿度应力场理论的数学实质是在初边界条件下求解由式（1）～（3）和式（5）组成的偏微分方程组。

1.2 基于相对湿度的扩散方程

湿度场与温度场具有相似性，但也存在一些实质性的区别，如：湿度必然有一个上限指标，即饱和含水率。当岩土体达到饱和后，不会自发地继续吸水，湿度也不会继续增加。同时，岩土体会固结排水，湿度也不仅仅是空间位置和时间的函数，还是应力状态的函数，严格地说，在湿度应力理论中，湿度场和应力场是双向耦合的。为简化耦合模式，本文对含水率行归一化处理，引入相对湿度的概念，将湿度场转化为相对湿度场来描述水分在岩土体的运移规律，相对湿度ϖ的定义为：

式中，w为含水率，wsat为饱和含水率。对于相对湿度，采用体积含水率和质量含水率是等效的。

相对湿度场的控制方程仍可简化为一个无源的扩散方程，即：

式中，K为相对湿度的扩散率，是求解瞬态问题的基础，理论上可通过非饱和土的渗流试验进行标定[22-23]。

1.3 考虑软化效应的本构方程

文献[8]在提出湿度应力场时就强调了岩石软化问题的重要性，即在考虑岩土体吸水膨胀的同时，不应忽略其弹性模量的降低。因此，式（5）中的弹性模量应当是含水率的函数，进而也是相对湿度的函数，而非一个固定的材料参数，即：

那么，式（5）也不再是线性的偏微分方程，应力的求解难度也将显著增加。

综上所述，基于相对湿度并考虑软化效应的湿度应力场瞬态分析模型控制方程为：

由式（9）可知，相对湿度扩散率K、膨胀系数α的取值和弹性模量E与相对湿度ϖ间的函数关系E(ϖ)都需要通过试验确定。对于稳态问题，不考虑具体的湿度场变化过程，不涉及湿度扩散方程，那么可以不对相对湿度扩散率K进行准确标定[10-11,15]。

2 湿度应力场稳态参数的确定

以采自湖南株洲的红层黏土（固结土，压实系数大于0.95，编号A）和采自甘肃天水的红层黏土岩（编号B）为研究对象，开展物理特性试验，结果见表1。

表1 土样的基础物理参数Tab.1 Basic physical parameters of samples

由表1可知，固结土和黏土岩的基础物理参数较为接近，且均处于硬塑状态，具有较强的可比性。

2.1 膨胀力及侧限膨胀试验

膨胀试验通过WG型三联高压固结仪完成，试验装置测得土样A、B的膨胀压力时程曲线如图1所示。

图1 室内试验得到的土样膨胀压力时程曲线Fig.1 History curves of expansive pressure in laboratory test

由图1可知：固结土的吸水速度较快，在30 min时膨胀压力已达到稳定；而黏土岩的膨胀压力在45 min时才达到稳定。分别选取试验中位值112.5 kРa和225.0 kРa作为固结土和黏土岩的膨胀应力进行后续分析。

在侧限膨胀试验中，以膨胀应力为参考值，固结土组的法向压力分别设定为0、12.5、25.0、50.0、75.0和100.0 kР；黏土岩组的法向压力分别设定为0、25、50、100、150和200 kРa，各进行3次平行试验，结果见表2和3。

表2 固结土侧限膨胀率Tab.2 Confined expansive rates of consolidated clay

表3 黏土岩侧限膨胀率Tab.3 Confined expansive rates of clay rock

2.2 稳态参数的标定方法

假设膨胀变形结束后试样已处于饱和状态，引入相对湿度后，膨胀系数α与无荷载作用时的侧限膨胀率δep0为同一参数的不同表述。膨胀系数可以通过侧限膨胀试验直接得到，固结土和黏土岩的膨胀系数分别为0.205和0.345。

基于文献[24]的修正Huder-Amberg模型，对表2中的数据进行拟合，得到侧限膨胀率ε与法向压力间的函数关系如式（10）所示，曲线如图2所示。

图2 侧限膨胀率-法向压力曲线Fig.2 Confined expansive ratio-normal stress curves

式中：εb为膨胀参数；σ0为膨胀应力；σn为施加的法向压力；σe为侧限膨胀试验的预加荷载，可取为σe=10 kРa。

由图2可知，试验值与理论曲线的一致性较好，可以对不同法向压力对应的侧限膨胀率进行预测。

饱和含水率可以通过初始孔隙率和侧限膨胀率按式（11）计算：

式中，e0为初始孔隙率，δep为侧限膨胀率。

根据式（11）可以在法向压力-饱和含水率坐标系内得到一系列散点，通过线性拟合建立两者的函数关系，如图3所示。由于在相同的法向压力条件下黏土岩具有更大的膨胀变形，因此也对应更大的饱和含水率。

图3 法向压力-饱和含水率曲线Fig.3 Normal stress-saturated water content curves

根据表2的数据，将不同法向压力状态视为饱和试样侧限排水压缩时的中间状态，按式（12）反演得到不同含水率阶段时的弹性模量：

式中，Δσn为侧限膨胀试验中法向压力的增量，Δδep为对应的侧限膨胀率增量。

对图3中的饱和含水率指标进行归一化处理，可以计算得到不同相对湿度下的压缩模量，如图4所示。侧限膨胀试验的相对湿度变化范围较小，因此本文还补充了塑限、缩限和干燥状态下的侧限压缩试验。

图4 弹性模量-相对含水率曲线Fig.4 Elastic modulus-relative humidity curves

假设弹性模量与相对湿度呈线性关系，那么基于图4中的数据即可建立变模量本构如式（13）所示：

式中，E1为饱和弹性模量，E0为饱和过程弹性模量的减小量，E0+E1为干燥弹性模量。

由图4可知：相对湿度相同时，黏土岩的压缩模量大于固结土；但随着含水率的增加，其压缩模量衰减速度也更快，在饱和时的压缩模量基本相同。

综合图3和4可知，黏土岩的膨胀性和软化性均高于固结土。

3 模型可靠性验证

由于相对湿度场与温度场的控制方程具有相同的形式，在数值模拟过程中可以通过在0～1范围内变化的温度场来等效模拟，并定义弹性模量随温度变化实现变模量本构。

3.1 验证性模型

通过膨胀应力和侧限膨胀率对本文模型的可靠性进行验证，建立的几何模型如图5所示。

图5 数值模拟中的几何模型Fig.5 Geometric models in numerical simulation

如图5所示，侧限膨胀模型包括试样和试验装置（环刀、垫片）两部分，其中，试样采用温变弹塑性材料进行模拟，试验装置采用刚体材料进行模拟。模型的环刀部分节点均为固定节点，试样底部施加法向位移约束。试样和环刀间不共节点，采用无摩擦接触；垫片底面与试样顶面也不共节点，采用点-面双向接触。在模拟侧限膨胀试验时，在垫片上表面施加均布压力，并对垫片法向位移数据进行监测和输出；在模拟膨胀应力试验时，对垫片顶部节点进行固定约束，监测试样与垫片间的接触力并输出。

初始时刻，将模型的相对湿度设置为0；定义上、下表面为恒定湿度边界，相对湿度设置为1；模型采用显式方法按瞬态问题进行求解，求解器为LS-solver（R11.0）。由于膨胀应力和侧限膨胀率均为稳态指标，因此并不需要准确的相对湿度扩散率Κ。为加速求解，计算时将相对湿度扩散率设置为10 m2/s。数值模拟得到无荷载的侧限膨胀率时程曲线如图6所示，膨胀压力时程曲线如图7所示。由于采用了放大的相对湿度扩散率，图6和7中横轴T对应数值计算时的虚拟时长，而不是实际时间t。

图6 无荷载侧限膨胀率时程曲线Fig.6 History curves of confined expansive rates

图7 膨胀压力时程曲线Fig.7 History curves of expansive pressure

如图6和7可知：当虚拟时长达到0.14 s时，膨胀应力模型的相对湿度达到稳定；而在膨胀率模型中，黏土岩和固结土相对湿度达到稳定的虚拟历时分别为0.16和0.19 s，略长于膨胀应力模型。这是由于模型采用了相同的相对扩散率，而膨胀率模型中体积膨胀增加了湿度扩散的距离。数值模拟得到的膨胀应力和侧限膨胀率与室内试验结果在数值上基本一致，也验证了本文方法在求解稳态问题时的可靠性。

3.2 相对湿度扩散率的标定

根据湿度场控制方程（7），增大湿度扩散率后，虚拟时长T和实际的渗透时间t之间的变换是线性的，且缩小的比例因子与湿度扩散率的放大比例相同。将图6的数值曲线与室内试验结果对比，通过相对湿度场达到稳定的历时来间接地标定相对湿度扩散率Κ。以固结土为例，膨胀压力稳定历时为30 min，而模拟中对应的虚拟历时为0.14 s，相对湿度扩散率放大比例为12 857，那么实际湿度扩散率应当为预设值的1/12 857，即7.8×10-4m2/s；同理，可以计算得到黏土岩的相对湿度扩散率为6.0×10-4m2/s。严格地说，本方法标定的相对湿度场湿度扩散率Κ是湿度场的平均扩散率，而非瞬态值。

3.3 模型的验证

除无荷载侧限膨胀试验外，本文进行了有法向压力的侧限膨胀数值模拟，得到的受荷侧限膨胀率时程曲线如图8所示，其中横坐标为已转换为实际时间t。

图8 有荷载的侧限膨胀率时程曲线Fig.8 History curves of loaded confined expansive rates

由图8可知，在有法向荷载作用时，数值模拟得到的稳定侧限膨胀率在数值上与室内试验的一致性仍较高，但模拟的侧限膨胀率时程曲线与试验结果存在一定差异，主要表现为膨胀变形稳定前实测的膨胀率大于模拟值。这是由于试验对象为非饱和的黏土（岩），其湿度扩散率与含水率密切相关，随含水率的增加，湿度扩散率减小[25]，而模拟过程中并未考虑这一问题。然而，由平均扩散率引起的差异较小，且会随时间的增加而持续减小，对模拟结果的影响并不明显。

4 算例

4.1 稳态理论解

文献[26]给出了用湿度应力理论求解圆形洞室遇水的算例，简化后的力学模型如图9（a）所示。图9（a）中：ρ为岩土体的密度；H为埋深；a为隧道半径；b为影响区范围，取b=5a；σr=a、σr=b分别为隧道壁和影响区边缘的应力边界条件，取σr=a=0，σr=b=ρgH。

图9 圆形洞室模型及稳态理论解Fig.9 Model of circular chambers with continuous internal water erosion and its theoretical solution

根据弹性力学小孔问题，考虑5倍半径区域的影响，并忽略重力梯度，采用均布的压应力ρgH代替，孔内侧为自由边界无荷载，并假设湿度场为线性分布，即：

式中，W0为洞室内边界湿度，a为洞室半径。

若采用相对湿度来描述，则为：

对应的线弹性理论解为[26]：

式中，σr、σθ为极坐标下的应力分量，约定以受拉为正，受压为负。

作为对比，本文采用相同的简化条件，并考虑由湿度场引起的软化效应，即，将式（13）代入式（16），得到对应的理论解为：

以固结土中半径a=2 m的圆形洞室为例，本文及文献[22]模型在5倍洞室半径的稳态理论解如图9（b）所示。对于式（16）的模型，为增加可比性，取弹性模量E=0.5[E1+(1/5)E0]，即本文模型变形模量变化范围内的中间值。

由图9可知，本文模型计算得到的洞室周边环向压应力分量较文献[22]的模型小，而在较远处环向应力分量则较大，这是由于本文模型考虑了围岩遇水软化引起的应力重分布。

4.2 瞬态分析

本文模型也可适用于瞬态问题的数值分析，瞬态分析模型可采用显-隐式方法求解，即：首先将相对湿度设置为0，施加应力边界条件σr=b=ρgH，通过隐式迭代法进行应力初始化；然后，定义内边界相对湿度为1，外边界相对湿度为0，并对应力初始化的模型进行显式的重启动求解。

围岩的材料参数在第2、3节中已经进行了标定，模拟中将渗透系数放大864 000倍进行求解，那么虚拟时长0.1 s即相当于真实时间1 d。仍取洞室半径a=2 m，埋深H=10 m。求解得到洞室在120 d内相对湿度的变化过程如图10所示。为便于比较，图10中给出了ϖ=0.167的参考线。

图10 数值模拟得到的相对湿度场变化规律Fig.10 Evolution of relative humidity field in numerical simulation

由图10可知：随着时间的增加，固结土和黏土岩内相对湿度场的变化均趋缓，在第60天至第120天阶段，黏土岩内相对湿度ϖ=0.167参考线移动距离仅为0.2 m；相同阶段，固结土内参考线的移动距离已小于0.1 m，但均未达到理论上的稳态。可见，岩土体吸水膨胀问题实际是一个持续的过程，相对湿度场在短时间内无法达到稳态，本算例作为瞬态问题分析更具有工程意义。

第3、15、30、60和120天时的σr和σθ随r的变化规律如图11所示。

图11 极坐标下应力分量的变化规律Fig.11 Evolution of stress components in the polar coordinate system

由于本文模型可以考虑围岩软化引起的应力重分布，在图11中，随着时间的增加，洞室周边的环向应力有所下降。对比固结土与黏土岩可以发现，黏土岩虽然具有更大的膨胀性，但由于其遇水软化明显，洞室周边岩体实际应力水平低于相同条件下的固结土洞室。且黏土岩内应力下降比例较高，在第3天至第60天阶段，环向最大压应力由385下降为331 kРa，下降比例为14.0%；而该阶段固结土内的环向最大压应力下降比例仅为9%，这是因为黏土岩软化引起的应力重分布更加明显。因此，进行膨胀岩（土）洞室稳定分析时，除围岩膨胀性外也应关注其软化性。由图11还可知，在围岩软化尚未稳定前，本文数值模拟得到的最大压应力已小于理论模型，这是由于式（17）的理论模型是基于小变形假设推导的，未考虑围岩大变形产生的卸载，计算应力偏大。