CN106202980B - 一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法 - Google Patents

Abstract

本发明公开了一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法，采用流‑固耦合模块并编制相关FISH语言程序，实现非饱和渗流过程中非饱和区基质吸力、非饱和渗透系数以及土体强度参数随含水率的变化关系；然后基于“湿度应力场理论”，根据渗流连续性微分方程和热传导方程的相似性，推导出热力学参数和渗流参数的等效转换关系，并采用热力学模块，实现膨胀土在增湿条件下的膨胀变形过程。本发明能够有效且全面地实现膨胀土在降雨入渗条件下增湿膨胀的全过程数值计算，从而为膨胀土工程灾害高风险地区的工程设计、施工提供有效的科学指导手段，同时有力地推动计算土力学在实际工程建设和地质灾害防治等技术领域的应用。

Description

一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法

技术领域

本发明属于计算土力学和工程灾害防治技术领域，涉及一种数值计算方法，尤其涉及一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法。

背景技术

全世界分布有膨胀土的国家有40多个，而相对其他国家，我国的膨胀土问题尤为突出。首先，膨胀土在我国分布范围很广，基本覆盖了西南云贵高原到华北平原之间的大片土地。其次，膨胀土种类较多，不同地区膨胀土的形成原因及工程特性不尽相同，具有显著的地域性。随着我国基础设施建设进程加速，许多工程建设将不可避免地穿越膨胀土地区。由于我国现有膨胀土地区建筑技术规范还不成熟和完善，膨胀土地区工程设计施工多依据以往工程经验。然而膨胀土具有显著的地域差异性，盲目套用其他地区膨胀土设计施工经验和参数将导致许多严重工程问题的产生，造成巨大的经济损失和人员伤亡。

众多膨胀土工程实例表明，降雨入渗是诱发膨胀土工程灾害事故发生的主要元凶。然而，由于其涉及复杂的非饱和渗流和膨胀变形过程，采用理论求解时，为计算方便，常对其实际情况进行过度简化，导致计算结果精度较低，且耗时耗力。虽然目前已经探索了相关膨胀土增湿膨胀变形的数值模拟方法，但仍存在一些不足，例如，基于“湿度应力场”理论，采用温度场模拟降雨增湿时，其实际上只适用于饱和渗流分析，不能考虑降雨过程中雨水重力以及非饱和区基质吸力和渗透系数变化的影响，因此不能真实反映膨胀土在降雨入渗条件下的非饱和渗流过程。同时在考虑膨胀土膨胀效应时，通常是通过在单元上施加垂直于临空面的外力来模拟膨胀力，但实际中膨胀变形是向四周发展的，膨胀力并不能单纯以外力方式施加。

综上所述，目前针对膨胀土增湿膨胀数值模拟研究中，均存在一定的问题，尚不能合理且全面地模拟膨胀土在增湿条件下的非饱和渗流和膨胀变形过程，计算结果与实际情况存在较大误差，不能很好地指导膨胀土地区工程的设计和施工。

发明内容

本发明的目的在于解决现有数值计算中不能真实反映膨胀土在增湿条件下膨胀变形的技术问题，进而提供了一种合理且全面考虑降雨入渗条件下非饱和渗流和膨胀变形过程的数值计算方法。

为了实现上述目的，本发明采用如下技术方案：

一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法，其特征在于：所述方法包括以下步骤：

1)建立膨胀土工程力学计算模型；采用FLAC^3D程序建立三维几何模型并划分计算网格，划分有计算网格的三维几何模型包括模型单元以及模型节点，所述膨胀土工程力学计算模型是根据具体膨胀土工程的实际情况建立的；

2)对步骤1)所建立的膨胀土工程力学模型进行初始平衡地应力计算；

3)结束步骤2)的计算后进行非饱和渗流计算：

4)结束步骤3)的计算后进行热传导计算；

5)对数值模拟结果进行分析。

作为优选，本发明所采用的步骤3)的具体实现方式是：

3.1)对模型单元设置相应的渗流参数，并根据具体膨胀土工程的实际情况设置初始饱和度和孔隙水压力，膨胀土工程力学模型的上边界设置为流量边界以模拟降雨入渗过程，并设置渗流计算时间；

3.2)在每一渗流计算时步中，提取模型节点饱和度S_r，并判断模型节点饱和度S_r是否等于1.0；若是，则直接进行步骤3.3)；若否，则依据通过试验得出的土-水特征曲线拟合公式计算该模型节点的基质吸力值，并将模型节点的基质吸力值赋值给节点的孔隙水压力后进行步骤3.3)；

3.3)首先通过反距离加权插值法，将模型节点饱和度转化为模型单元饱和度，然后根据渗透系数和强度参数与饱和度之间的关系计算相应饱和度下模型单元的渗透系数以及强度参数值，并将相应饱和度下模型单元的渗透系数以及强度参数值赋给模型单元；

3.4)检查地表节点的孔隙水压力pˊ，若孔隙水压力大于0，将该地表节点的流量边界修改为压力边界，固定该地表节点的孔隙水压力为0；若地表节点的孔隙水压力小于或等于零，则直接进行步骤3.5)，所述地表节点是膨胀土工程力学模型上边界节点；

3.5)判断是否达到非饱和渗流计算的终止条件，若是，则结束计算过程并保存结果文件；若否，重复步骤3.2)-3.4)直至达到非饱和渗流计算的终止条件；所述非饱和渗流计算的终止条件是渗流计算时间。

作为优选，本发明所采用的步骤4)的具体实现方式是：

4.1)基于湿度应力场理论，根据渗流连续性微分方程和热传导微分方程的相似性，建立渗流参数与热力学参数的等效对应关系；

其中，所述渗流连续性微分方程的表达形式是：

所述热传导微分方程的表达形式是：

式中：

k_x、k_y以及k_z分别为x、y以及z这三个方向的渗透系数；

h_m是基质吸力水头；

C_w是比水容重；

λ_x、λ_y以及λ_z分别为x、y以及z这三个方向的热传导系数；

T是温度；

ρ是介质密度；

C_v是介质的比热容；

t是时间；

比较渗流连续性微分方程以及热传导微分方程可知，渗流参数与热力学参数存在等效对应关系，其中，渗透系数k_i对应热传导系数λ_i；基质吸力水头h_m对应温度T；比水容重C_w对应比热容ρC_v；根据以上对应关系，确定出热传导计算中模型单元和模型节点的热力学参数；

4.2)提取步骤3)计算结束后保存的结果文件，给模型单元赋予上述按等效原理计算出的热力学参数；

4.3)设定当模型单元含水率从残余含水率θ_r增加到饱和含水率θ_s时对应的温度为100℃，各模型节点对应的等效温度T_p采用线性插值法计算，计算公式如下：

4.4)固定降雨入渗边界处模型节点的温度，并进行热传导计算，在每一热力计算时步中，判断各模型节点温度是否全部达到其等效温度T_p，若是，则结束计算过程；若否，则重复进行步骤4.3)直至达到等效温度T_p后结束计算。

作为优选，本发明所采用的步骤2)中是FLAC^3D程序采用显式有限差分计算方法进行迭代结算，在迭代过程中通过监控最大不平衡力比率，检测计算是否达到收敛，一旦最大不平衡力比率小于预先设置的值，则计算过程终止。

作为优选，本发明所采用的步骤5)中的计算结果分析包括但不限于不同降雨历时下围岩含水率分布、孔隙水压力分布、围岩应力以及工程结构的受力和变形。

与现有技术相比，本发明的优点是：

(1)通过本发明的数值计算方法，能够全面模拟膨胀土在降雨入渗条件下的增湿膨胀变形全过程，克服了现有数值计算方法的不足和缺陷，从而能够精确计算膨胀土在增湿膨胀变形过程中对工程结构受力和变形的影响，揭示膨胀土工程灾害孕育机理和演化规律，为膨胀土地区相关工程的设计和施工提供参考和依据。

(2)本发明的数值计算方法具有简单、灵活以及通用性强的特点，方法便于实现，既能保证计算精度，又能节约时间成本。

附图说明

图1是本发明膨胀土增湿膨胀数值模拟方法流程图；

图2是隧道平面力学模型网格划分图；

图3是采用流-固耦合模块模拟非饱和渗流过程的数值计算流程图；

图4是采用热力学模块模拟膨胀变形过程的数值计算流程图；

图5是隧道初期支护监测点布置示意图；

图6是初期支护弯矩随降雨历时变化曲线；

图7是降雨过程中初期支护变形转化示意图。

具体实施方式

下面结合附图和具体的实施例对本发明的膨胀土在降雨入渗条件下增湿膨胀数值模拟方法做进一步的详细说明：

膨胀土在降雨入渗条件下增湿膨胀数值模拟方法流程图见附图1，本发明以某浅埋膨胀土隧道工程为例，对其在降雨入渗条件下支护结构的动力响应进行模拟分析，具体按照以下步骤实施：

1、建立浅埋膨胀土隧道力学计算模型；

某山岭隧道为浅埋膨胀土隧道，单洞净宽10.5m，净高5.0m，洞口段埋深8.4m。为简化计算，考虑围岩处于平面应变状态，取模型厚度1m。为减小边界效应的影响，模型宽度取100m，隧道底部至模型下边界35m，埋深8.4m，如附图2所示。对浅埋膨胀土隧道力学计算模型四周施加法向约束，并将底部固定，浅埋膨胀土隧道力学计算模型四周及底部均设置为不透水边界，顶部设置为降雨边界。围岩采用Mohr-Coulomb模型，初期支护采用Shell单元，厚度为26cm，其弹性模量由钢拱架、钢筋网片和喷射混凝土按刚度等效原理进行折算，围岩及支护结构的力学参数列于表1，其中土体饱和渗透系数k_sat取为10^-4cm/s，非饱和区渗透系数是关于饱和度的函数，可采用以下经验公式进行计算。

k(s)＝k_sats²(3-2s) (1)

表1围岩及支护结构的力学参数

	密度/(kg/m<sup>3</sup>)	弹性模量/GPa	泊松比
				围岩	2000	0.08	0.35
锚杆	—	210	0.25
				初支	2500	23	0.20

采用滤纸法测定现场膨胀土样的基质吸力，并以Van Genuchten(V-G)提出的土-水特征曲线模型进行参数拟合，V-G模型表达式及相关拟合参数如下：

式中：θ为含水率，θ_s为饱和含水率，θ_r为残余含水率，ψ为基质吸力。拟合参数a＝0.0116，n＝2.2610，m＝0.2841，θ_r＝9.02％，θ_s＝31.07％。

通过室内直剪试验测得现场膨胀土样的黏聚力c和内摩擦角随含水率θ变化的拟合关系式为：

2、进行浅埋膨胀土隧道力学计算模型初始平衡状态求解计算；

FLAC^3D程序采用显示有限差分计算方法进行迭代结算，在迭代过程中通过监控最大不平衡力比率，检测计算是否达到收敛，一旦最大不平衡力比率小于预先设置的值，则计算过程终止。

3、进行非饱和渗流计算；

采用FLAC^3D内置的流-固耦合模块进行非饱和渗流计算的流程图如附图3所示，具体过程表述如下：

①对模型单元设置相应的渗流参数，并根据实际情况设置初始饱和度和孔隙水压力，并将模型上边界设置为降雨入渗边界，本实例中假设降雨强度较大，地表在短时间内便达到饱和状态，进入压力入渗阶段，因此设置地表孔隙水压力为零，并设置渗流计算时间为72h。

②在每一渗流计算时步中，提取模型节点饱和度S_r，若模型节点饱和度不等于1.0，则依据土-水特征曲线拟合公式求得该模型节点的基质吸力值，并将该模型节点的基质吸力值赋值给相应节点的孔隙水压力(负孔隙水压)，若饱和度等于1.0，则直接进行下一步计算；

③由于强度参数以及渗透系数是单元变量，而饱和度是节点变量，因此首先通过反距离加权插值法，将模型节点饱和度转化为模型单元饱和度，然后计算相应饱和度下模型单元的渗透系数以及强度参数值，并将其值赋给模型单元。

④检查地表节点的孔隙水压力pˊ，若孔隙水压力大于0，将该地表节点的流量边界修改为压力边界，固定其孔隙水压力为0。(由于本实例中将地表空隙水压力设置为零来模拟强降雨情况，因此可跳过此步骤)

⑤检查是否达到渗流计算终止条件(预先设定的渗流计算时间72h)，若未达到则重复步骤②～⑤，否则结束计算。

4、进行热传导计算；

采用FLAC^3D内置的热力学模块进行膨胀变形计算的流程图如附图4所示，具体过程表述如下：

①基于“湿度应力场理论”，根据渗流连续性微分方程和热传导微分方程的相似性，建立渗流参数与热力学参数的等效转换关系。其中，渗流连续性微分方程和热传导微分方程的表达形式如下：

式中：k_x、k_y、k_z分别为三个主方向的渗透系数，h_m为基质吸力水头，C_w为比水容重，λ_x、λ_y、λ_z分别为三个主方向的热传导系数，T为温度，ρ为介质密度，C_v为介质的比热容。

比较式(4)和式(5)可知，渗流参数与热力学参数存在等效关系，其中，渗透系数k_i对应热传导系数λ_i；基质吸力水头h_m对应温度T；比水容重C_w对应比热容ρC_v。根据以上对应关系，可确定出热传导计算中，模型单元和节点对应的热力学参数。

②提取非饱和渗流计算结果文件，给模型单元赋予上述按等效原理计算出的热力学参数(λ_i、C_v、α)。

③设定当模型单元含水率从残余含水率θ_r增加到饱和含水率θ_s时对应的温度为100℃，各模型节点对应的等效温度T_p采用线性插值法计算，计算公式如下：

④固定降雨入渗边界处模型节点的温度，并进行热传导计算，检查各模型节点温度是否全部达到其等效温度，如未达到则继续热传导计算，否则，结束计算。

5、对数值模拟结果进行分析。

本实施步骤主要包括：降雨入渗过程中围岩孔隙水压力和围岩应力分布变化规律分析以及支护结构受力和变形情况分析。

(1)降雨入渗过程中围岩孔隙水压力和围岩应力分布变化规律分析

依据上一步计算所得的不同降雨历时下围岩孔隙水压力以及围岩应力分布图，分析降雨入渗过程中围岩孔隙水压力和围岩应力分布变化规律。可以看出随降雨历时增加，浸润锋线向下部推进，降雨影响范围内围岩孔隙水压力升高(基质吸力减小)。降雨过程中围岩垂直应力变化不大，但水平应力显著增大，降雨前拱顶处水平应力为0.06MPa，降雨72h后拱顶水平应力增加至0.33MPa，且在隧道拱顶上方产生应力集中。

(2)支护结构受力和变形情况分析

取衬砌结构上5处进行监测，监测点布置如附图5所示，监测点处弯矩随降雨历时变化曲线如附图6所示，其中正值表示衬砌临空一侧受拉，负值表示衬砌围岩一侧受拉。从图中可看出，降雨增湿过程中，衬砌在监测点4(即拱脚)处弯矩绝对值最大。拱顶处弯矩随降雨历时增加逐渐减小至0，后反向增大；而拱腰及边墙处逐渐由负弯矩转变为正弯矩。这是由于在降雨初期，垂直荷载较大，支护结构在上覆荷载作用下发生竖向挤压变形，此时拱顶处支护临空一侧受拉，拱腰和边墙处支护围岩一侧受拉；降雨持续一段时间后，围岩增湿膨胀效应引起水平压力逐渐增大，支护结构逐渐由竖向挤压变形转化为水平向挤压变形，这使得拱顶、拱腰及边墙处衬砌弯矩发生正负值变化，降雨入渗过程中初期支护结构变形转化示意图如附图7所示。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照实施例对本发明进行了详细说明，凡对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种膨胀土在降雨入渗条件下增湿膨胀数值模拟方法，其特征在于：所述方法包括以下步骤：

1)建立膨胀土工程力学计算模型；采用FLAC^3D程序建立三维几何模型并划分计算网格，划分有计算网格的三维几何模型包括模型单元以及模型节点，所述膨胀土工程力学计算模型是根据具体膨胀土工程的实际情况建立的；

2)对步骤1)所建立的膨胀土工程力学模型进行初始平衡地应力计算；

3)结束步骤2)的计算后进行非饱和渗流计算：

3.1)对模型单元设置相应的渗流参数，并根据具体膨胀土工程的实际情况设置初始饱和度和孔隙水压力，膨胀土工程力学模型的上边界设置为流量边界以模拟降雨入渗过程，并设置渗流计算时间；

3.2)在每一渗流计算时步中，提取模型节点饱和度S_r，并判断模型节点饱和度S_r是否等于1.0；若是，则直接进行步骤3.3)；若否，则依据通过试验得出的土-水特征曲线拟合公式计算该模型节点的基质吸力值，并将模型节点的基质吸力值赋值给节点的孔隙水压力后进行步骤3.3)；

3.3)首先通过反距离加权插值法，将模型节点饱和度转化为模型单元饱和度，然后根据渗透系数和强度参数与饱和度之间的关系计算相应饱和度下模型单元的渗透系数以及强度参数值，并将相应饱和度下模型单元的渗透系数以及强度参数值赋给模型单元；

3.4)检查地表节点的孔隙水压力pˊ，若孔隙水压力大于0，将该地表节点的流量边界修改为压力边界，固定该地表节点的孔隙水压力为0；若地表节点的孔隙水压力小于或等于零，则直接进行步骤3.5)，所述地表节点是膨胀土工程力学模型上边界节点；

3.5)判断是否达到非饱和渗流计算的终止条件，若是，则结束计算过程并保存结果文件；若否，重复步骤3.2)-3.4)直至达到非饱和渗流计算的终止条件；所述非饱和渗流计算的终止条件是渗流计算时间；

4)结束步骤3)的计算后进行热传导计算：

4.1)基于湿度应力场理论，根据渗流连续性微分方程和热传导微分方程的相似性，建立渗流参数与热力学参数的等效对应关系；

其中，所述渗流连续性微分方程的表达形式是：

所述热传导微分方程的表达形式是：

式中：

k_x、k_y以及k_z分别为x、y以及z这三个方向的渗透系数；

h_m是基质吸力水头；

C_w是比水容重；

λ_x、λ_y以及λ_z分别为x、y以及z这三个方向的热传导系数；

T是温度；

ρ是介质密度；

C_v是介质的比热容；

t是时间；

比较渗流连续性微分方程以及热传导微分方程可知，渗流参数与热力学参数存在等效对应关系，其中，渗透系数k_i对应热传导系数λ_i；基质吸力水头h_m对应温度T；比水容重C_w对应比热容ρC_v；根据以上对应关系，确定出热传导计算中模型单元和模型节点的热力学参数；

4.2)提取步骤3)计算结束后保存的结果文件，给模型单元赋予上述按等效原理计算出的热力学参数；

4.3)设定当模型单元含水率从残余含水率θ_r增加到饱和含水率θ_s时对应的温度为100℃，各模型节点对应的等效温度T_p采用线性插值法计算，计算公式如下：

4.4)固定降雨入渗边界处模型节点的温度，并进行热传导计算，在每一热力计算时步中，判断各模型节点温度是否全部达到其等效温度T_p，若是，则结束计算过程；若否，则重复进行步骤4.3)直至达到等效温度T_p后结束计算；

5)对数值模拟结果进行分析。

2.根据权利要求1所述的膨胀土在降雨入渗条件下增湿膨胀数值模拟方法，其特征在于：所述步骤2)中是FLAC^3D程序采用显式有限差分计算方法进行迭代结算，在迭代过程中通过监控最大不平衡力比率，检测计算是否达到收敛，一旦最大不平衡力比率小于预先设置的值，则计算过程终止。

3.根据权利要求2所述的膨胀土在降雨入渗条件下增湿膨胀数值模拟方法，其特征在于：所述步骤5)中的计算结果分析包括但不限于不同降雨历时下围岩含水率分布、孔隙水压力分布、围岩应力以及工程结构的受力和变形。