CN105160138A - 隧道流固耦合稳定性分析方法 - Google Patents

Abstract

本发明提供一种隧道流固耦合稳定性分析方法。本发明方法，包括：根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，所述弹性阶段的状态指标包括剪切状态和拉伸状态；根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标，破坏临界点的等效塑性应变为极限等效塑性剪切应变或极限等效塑性拉伸应变；根据所述弹性阶段的状态指标、所述屈服阶段以及破坏阶段的状态指标确定所述材料单元耦合过程中渗透系数的变化值；根据所述渗透系数进行隧道的流固耦合稳定性分析。本发明提高了隧道流固耦合稳定性分析的准确率。

Description

隧道流固耦合稳定性分析方法

技术领域

本发明实施例涉及数字处理分析领域，尤其涉及一种隧道流固耦合稳定性分析方法。

背景技术

富水区隧道开挖是一个复杂的流固耦合过程，一直以来，地下隧道围岩的稳定性分析并没有十分科学合理的方法，很多停留在以某点位移、塑性区大小的经验值或者整体的安全系数作为判断稳定性的依据。然而这样的判别不能完全反映出岩石在实际复杂应力状态下的危险性。

富水区盾构隧道围岩的整个的渐进破坏过程包含了弹性、屈服、破坏各阶段，其中伴随着复杂的渗流-应力耦合作用。应力、损伤及破坏的不同单元状态会对渗透性产生不同的影响，而传统流固耦合分析忽略不同单元状态渗透性的变化，因此带来很大误差。

发明内容

本发明实施例提供一种隧道流固耦合稳定性分析方法，以克服现有技术中流固耦合分析存在误差，准确率低的问题。

本发明的技术方案，包括：

根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，所述弹性阶段的状态指标包括剪切状态和拉伸状态；

根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标；

根据所述弹性阶段的状态指标、所述屈服阶段以及破坏阶段的状态指标确定所述材料单元耦合过程中渗透系数的变化值；

根据所述渗透系数进行隧道的流固耦合稳定性分析。

进一步地，所述根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，包括：

将所述应力状态点进行受力分析；

根据应力状态下的所述应力状态点处于同一平面上的最稳定参考点沿罗德角方向到屈服面的距离与所述应力状态点沿相同罗德角方向到屈服面的距离确定所述状态指标为

其中，所述ZSI为状态指标，所述O′B为空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离，所述O′P为该点在相同罗德角方向上到屈服面的距离，所述τ'_π为子午面GBO′面上B点纵坐标，所述τ_π为子午面GBO′面上P点纵坐标，所述I₁为应力张量第一不变量，所述为内摩擦角，所述θ_σ为应力罗德角，所述J₂为偏应力张量第二不变量，所述c为粘聚力。

进一步地，所述根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标，包括：

根据等效塑性剪切应变和极限等效塑性剪切应变确定所述材料单元在剪切屈服段内距离破坏段的接近程度；

根据所述接近程度确定屈服阶段和破坏阶段的状态指标为

$ZSI=min(1-\frac{{\mathrm{\ϵ}}^{ps}}{{\stackrel{\‾}{\mathrm{\ϵ}}}^{ps}},1-\frac{{\mathrm{\ϵ}}^{pt}}{{\stackrel{\‾}{\mathrm{\ϵ}}}^{pt}})---\left(2\right)$

其中，所述ε^ps为等效塑性剪切应变，所述为极限等效塑性剪切应变，所述ε^pt为等效塑性拉伸应变，所述为极限等效塑性拉伸应变。

进一步地，所述渗透系数的变化值为

$k=\left\{\begin{array}{cc}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI\&GreaterEqual;1\\ {\mathrm{\&xi;k}}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& 0\&le;ZSI<1\\ {\mathrm{\&xi;}}^{\&prime;}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI<0\end{array}\right.---\left(3\right)$

其中，所述k为当前的渗透系数，所述k₀为初始渗透系数，所述ε_v为体应变，所述n₀为初始孔隙度，所述ξ为屈服阶段突跳系数，所述ξ'为破坏阶段突跳系数。

本发明在H-M耦合过程中通过针对岩石弹性、屈服、破坏全阶段整体的、统一的安全性或危险性量化评价的指标建立岩土材料从弹性到破坏全阶段的定量评价方法，以单元状态指标ZSI值表征岩体单元的安全或危险程度。并通过值域能够直接明确岩石所处状态的评价方法，使用负值进行破坏阶段的评价，以便统一流固耦合体评价体系，提高了流固耦合分析的准确率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明隧道流固耦合稳定性分析方法流程图；

图2为本发明应力空间上应力点状态示意图；

图3为本发明H-M耦合模拟仿真流程图；

图4为本发明开挖面附近孔隙水压、渗透系数等值线图与渗流矢量示意图；

图5为本发明与中心点不同距离单元ZSI和渗透系数变化曲线；

图6a和图6b为本发明数值模拟结果与实测值对比示意图；

图7a、图7b以及图7c为本发明开挖面的横纵断面ZSI等值线图；

图8为本发明开挖面破坏率变化曲线图；

图9为本发明支护压力比对开挖面破坏率的影响。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明隧道流固耦合稳定性分析方法流程图，如图1所示，本实施例方法，包括：

步骤101、根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，所述弹性阶段的状态指标包括剪切状态和拉伸状态；

步骤102、根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标；

步骤103、根据所述弹性阶段的状态指标、所述屈服阶段以及破坏阶段的状态指标确定所述材料单元耦合过程中渗透系数的变化值；

步骤104、根据所述渗透系数进行隧道的流固耦合稳定性分析。

进一步地，所述根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，包括：

将所述应力状态点进行受力分析；

根据应力状态下的所述应力状态点处于同一平面上的最稳定参考点沿罗德角方向到屈服面的距离与所述应力状态点沿相同罗德角方向到屈服面的距离确定所述状态指标为

其中，所述ZSI为状态指标，所述O′B为空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离，所述O′P为该点在相同罗德角方向上到屈服面的距离，所述τ'_π为子午面GBO′面上B点纵坐标，所述τ_π为子午面GBO′面上P点纵坐标，所述I₁为应力张量第一不变量，所述为内摩擦角，所述θ_σ为应力罗德角，所述J₂为偏应力张量第二不变量，所述c为粘聚力。

具体来说，岩土工程通常针对整体的安全性，根据安全系数进行评价。然而，岩土体的整体失稳往往是由于局部范围内岩体的失稳开始的，因此，定量评价局部岩土体单元的安全性与危险性对于分析整体稳定性有重要的参考价值。材料单元在弹性阶段单元的塑性应变为0，如图2所示，材料单元实际的应力状态点P计算状态指标。当σ₁≤0时，由于没有拉应力存在，按照单元处于弹性剪切状态来考虑。Mohr-Coulomb准则屈服面在主应力空间中是一个不规则的六角形截面的角锥体表面，在π平面内投影为不等角六边形屈服面。I₁为应力张量第一不变量，J₂为偏应力张量第二不变量。点状线为偏平面上σ₂′轴的垂线，与σ₂′轴的夹角等于应力罗德角θ_σ。O′点为等倾线上的点，亦为相对最安全的参考点，连接O′P并延长交EFG平面于B，在GBO′面上O点坐标为(0，0)，O′点坐标为(σ_π，0)，σ_π为π平面上的正应力分量。P点的坐标为(σ_π，τ_π)。B点坐标为(σ_π，τ_π′)。此时ZSI可表述为：空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离与该点在相同罗德角方向上沿到屈服面的距离之比。

此时ZSI可表述为：空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离与该点在相同罗德角方向上沿到屈服面的距离之比，因此，

${\mathrm{\&tau;}}_{\mathrm{\&pi;}}=\sqrt{2J_2}---\left(4\right)$

根据公式(2)、公式(3)以及公式(4)得到该材料单元弹性阶段的状态指标。

所述弹性阶段的状态指标包括剪切状态和拉伸状态：当应力点在屈服面上时，ZSI＝1，应力点在等倾线上时ZSI＝+∞，表示单元在该阶段所能达到的最安全状态。需要说明的是，ZSI＝+∞并不是指单元无限的安全，而是在弹性剪切阶段处于相对最安全的状态。σ₁≤0，对应剪切状态，σ₁>0，对应与拉伸状态。岩土材料的抗拉强度相对抗剪强度而言要小很多，因此必须考虑材料在拉伸状态下的安全性，此时ZSI＝σt/σ1。σ1趋近于0，受拉状态为相对最安全状态ZSI＝+∞，然而此时仍不能排除单元剪切破坏的可能，为更好的把握岩土体的危险性，ZSI值可取为剪切与拉伸同时考虑的情况下的较小值。弹性段的ZSI表达式为：

弹性阶段ZSI∈[1，+∞)，ZSI大于1时，单元处于弹性阶段，ZSI值越大，安全性越高，当ZSI＝1时，单元开始进入屈服阶段。

进一步地，所述根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标，包括：

根据ε^ps/确定所述材料单元在剪切屈服段内距离破坏段的接近程度；

根据所述接近程度确定屈服阶段和破坏阶段的状态指标为

$ZSI=min(1-\frac{{\mathrm{\&epsiv;}}^{ps}}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^{ps}},1-\frac{{\mathrm{\&epsiv;}}^{pt}}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^{pt}})---\left(6\right)$

其中，所述ε^ps为等效塑性剪切应变，所述为极限等效塑性剪切应变，所述ε^pt为等效塑性拉伸应变，所述为极限等效塑性拉伸应变。

具体来说，屈服阶段的岩体出现了塑性变形，等效塑性应变大于0。采用应用应变来建立屈服与破坏条件更能反映材料屈服和破坏的本质。该阶段将破坏临界点的等效塑性应变值为破坏的判据，ε^ps/表示材料在剪切屈服段内距离破坏段的接近程度，其相补参量1-ε^ps/作为剪切屈服状态的单元状态指标。这样处理的优点是可以将屈服阶段的值域与弹性阶段明显的区分开来，并且仍然保持值越小，安全性越低的趋势。如考虑拉伸情况则ZSI值取两种情况中的较小值，即公式(6)。此时ZSI∈[0，1)，ZSI越大单元相对越安全，与弹性段的单调性相一致。材料刚刚进入屈服段时，ε^ps＝0，ZSI＝1。到达等效塑性应变阈值时，ZSI＝0。

破坏阶段时，剪切与拉伸状态的ZSI表达式与屈服阶段ZSI表达式相同。ZSI∈(-∞，0)，值越小破坏程度越高。

进一步地，所述渗透系数的变化值为

$k=\left\{\begin{array}{cc}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI\&GreaterEqual;1\\ {\mathrm{\&xi;k}}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& 0\&le;ZSI<1\\ {\mathrm{\&xi;}}^{\&prime;}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI<0\end{array}\right.---\left(7\right)$

其中，所述ε^ps为等效塑性剪切应变，所述为极限等效塑性剪切应变，所述ε^pt为等效塑性拉伸应变，所述为极限等效塑性拉伸应变。

具体来说，FLAC3D在模拟渗流的过程中，默认单元的渗透系数是不变的，然而这与实际的渗流情况并不相符。岩石的渗透性是与应力应变状态密切相关的。岩石在进入屈服阶段之前，即弹性阶段时，渗透系数将保持在较低水平上，而一旦进入屈服阶段，岩石的渗透系数将远远大于弹性阶段的渗透系数，而岩石破裂后渗透系数存在突跳现象，目前也较难找到连续光滑的函数表达。由于试件越接近均质，渗透系数突跳点越接近峰值点，可以认为细观单元突跳点和峰值点一致，通过细观单元的非均匀赋值，反映宏观试件应力应变-渗透系数的非线性。关于渗透系数与应力、应变关系的方程有很多。推导方法包括经验公式、间接公式和理论模型等方法。但大多数方程参数较多，适应性不好。体应变能更好的反应单元屈服、软化和破坏过程中渗透系数的变化。也即，初始压密阶段，ZSI≥1，体应变为负值，渗透性减小但变化并不明显；当应力达到峰值时，进入屈服状态，0≤ZSI<1，产生大量的微裂纹，体应变迅速增加；当单元破坏后，ZSI<0，裂隙扩展、贯通，变形随应力迅速增长，裂隙贯通形成畅通的导水通道，渗透系数突跳。随着变形的进一步发展，破裂的凹凸部分被剪断或磨损，在围压作用下，破坏试件又出现一定程度的压密闭合。本实施例中材料单元在弹性阶段到破坏阶段渗透系数k是与体积应变ε_v相关的函数，结合单元状态指标ZSI基于Kozeny-Carman公式表达为公式(7)，其中孔隙度的演化方程为：

$n=\frac{n_0+{\mathrm{\&epsiv;}}_v}{1+{\mathrm{\&epsiv;}}_v}---\left(8\right)$

该方程不仅表达形式简单，参数明确，更容易在仿真软件中实现。由于破裂后的单元渗透系数与弹性或屈服状态相差若干个等级，屈服状态比弹性状态同样高出很多倍，增速比体应变突变速率要快，因此用两个阶段分别用突跳系数ξ、ξ'来表征。ξ、ξ'取决于岩体性质。

采用FLAC3D建立数值模型仿真验证本发明的渗透系数，模拟仿真的流程如图3所示，以大连地铁2号线202标段盾构区间段包括促进路始发井-春光街站为例，春光街站-香工街站，香工街站-沙河口火车站三区间，三盾构区间单线全长2.3km。采用双线圆型盾构隧道。2台土压平衡式盾构机对左右两线分别进行施工。盾构区间地质构造较为复杂，隧道围岩表现为上软下硬。以桩号为DK14+634～DK14+670段为研究对象，软土的特性及地质勘探资料如下：(1)素填土：稍湿-饱和，松散-稍密状态，主要由粘性土、碎石等组成，层厚0.40～5.00m；(2)粉质粘土：含10-20％左右石英岩碎石、角砾，湿，可塑状态，层厚2.70～10.00m；(3)全风化钙质板岩：风化节理裂隙极发育，层厚1.40～7.80m；(4)强风化钙质板岩：碎裂结构，揭露层厚0.50～10.20m。(5)中风化钙质板岩：层理和节理裂隙较发育，矿物主要为云母、石英、方解石。勘察期间稳定地下水位为地下2.6m～8.7m，年水位变化幅度为1～3m。

该数值模型共16920个单元，18955个节点。模型水平方向为30m，隧道长度为36m，模型高度为37m。侧向施加法向约束，底部固定约束。布置如图所示的地表沉降监测点，测点间距为7.2m。模拟计算采用如下假设：

(1)土体本身变形与时间无关；(2)渗流是通过开挖面的透水实现的，渗流模型为各向同性渗流模型；(3)围岩为各向同性、连续的弹塑性材料，服从Mohr-Coulomb屈服准则。

一、计算参数和开挖步骤：

由于地质条件复杂，对分析地层进行简化，通过现场勘查资料和部分室内试验获得围岩参数如表1。计算模型为应变软化模型，岩层进入塑性后参数折减，土层不进行折减。岩石残余强度的粘聚力和内摩擦角按照折减系数1.1取值，折减公式如下：

c_F＝c₀/F_trial(9)

其中，所述c_F为折减后粘聚力，所述为折减后内摩擦角，所述c₀为初始粘聚力，所述为初始内摩擦角，所述F_trial为折减系数，c_r为残余粘聚力，为残余内摩擦角。所述c₀与c_r，与之间随等效塑性剪切应变ε^ps的增加进行线性插值。表1为各围岩层及支护的参数取值。

表1

盾构循环进尺为2.4m，共15步开挖完成。管片采用结构单元模拟，材料为C50钢筋混凝土，外半径为3m，每环长度为1.2m，底部管片施加160kPa的施工荷载。地下水埋深3m，水面为自由边界，模型侧面及底部，隧洞四周为不透水边界。管片外围为注浆层，由于地层较为软弱，采用全断面径向注浆加固。注浆圈厚度为2m。根据现场勘查资料，风化岩石层综合渗透系数k_r＝2.3×10^-2m/d，土层综合渗透系数为k_s＝0.1m/d，FLAC^3D中渗透系数与达西定律中的渗透系数不同，需乘换算系数才能用于计算^[23]，因此岩层换算后的渗透系数k_r＝2.71×10^-11(m²/P_a·sec)，初始孔隙率n₀＝0.5，ξ＝5，ξ′＝138。土层渗透系数则为k_s＝1.2×10^-10(m²/P_a·sec)，本文施加的开挖面支护压力为梯形荷载，如图3，计算采用的支护压力比为0.7。

σ＝λ(σ_p+σ_s)(11)

其中，所述λ为支护压力比，σ_p为孔隙水压力，σ_s为水平静止土压力，σ_p和σ_s根据隧道埋深和水头高度计算。

二、计算结果与分析：

1、渗流特性结果分析

隧道掘进的过程中，围岩的力学行为不断发生变化，渗透性也便发生了改变，图4为开挖面附近渗透矢量图。由于开挖面出现了破坏区域(ZSI<0)和屈服区域(0≤ZSI<1)，渗透系数在开挖面底部和中心处发生突跳，比没有破坏的部位渗透系数要高两个数量级左右，破坏区周围的屈服区域渗透性也有所增加，从渗流矢量图中可以看出，进入破坏阶段的单元渗流速度远大于其他区域。也验证了应变-渗透系数FISH程序的有效性。

单元的状态与渗透系数是紧密相关的，在隧道开挖面中心点沿x方向，选取等间距的8个单元，提取每个单元的渗透系数与ZSI值如图5所示，沿着x方向，渗透系数随着距离的增加先减小后增大，而ZSI随着距离的增加先增大后减小，两者的变化趋势相反。前2个单元ZSI在0到1之间，处于屈服状态，渗透系数明显大于弹性压密状态的第3、4、5、6单元，而第7、8单元的屈服程度相对较大，趋近于破坏状态，因此渗透系数的渗透系数更大，由于该方向上没有ZSI小于0的单元，各单元的渗透系数没有出现更大的突跳。计算结果与渗透性方程的变化理论均是相符的。

2、位移结果分析

图6a为数值计算沉降，图6b为实测沉降。显示随着隧道的掘进，监测点均产生持续的沉降，施工完成后各监测点的沉降值相差不大，累计沉降值约为-11mm。随着开挖面逐渐远离，DB-1和DB-2的沉降速率呈现先增加后减缓的趋势，而DB-3、DB-4的沉降速率随着开挖面的临近，出现先期减缓，后期增加的趋势。所有的监测点与实测值趋势相同，由于模型的简化、地质条件的复杂性，监测点位置间隔及监测时间的间隔误差，数值略有差异，但整体趋势较为吻合，也验证了模拟的正确性。在模拟的过程中，监测点沉降在达到最终稳定之前会有小范围的上抬，表明施工对周围土体扰动造成的隧道瞬时隆起存在滞后性。

表2为4个监测点在不同的计算方法得出的最终沉降量，其中不考虑地下水与实测沉降值相差最大，其次为H-M耦合但不考虑渗透性变化的情况，与实际监测最为接近的是H-M耦合并且渗透性变化的情况。考虑H-M耦合作用的沉降值要远远高出单独力学计算得出的沉降值。考虑渗透性变化模式的最终沉降也要略大于渗透性不变化的模式，因此地下水的影响是不可忽略的，同时考虑渗透性的变化是较为符合实际的，在富水区盾构隧道的开挖模拟中应当值得注意。

表2

3、单元状态指标分析

根据本次盾构隧道ZSI计算情况可知，ZSI<0的区域(破坏区)均集中在开挖面及前方的待开挖土体上。破坏区体积沿隧道掘进方向呈三角形逐渐减小，对隧道的每个横断面进行剖分，开挖面破坏面积最大，因此开挖面上的破坏面积的大小对于分析稳定程度具有一定的代表性。对于为便于分析提出开挖面破坏率的概念，即：开挖面上单元状态指标ZSI<0的单元总面积之和与开挖面断面面积之比。

$\mathrm{\&eta;}=\frac{A_Z}{A_0}---\left(12\right)$

其中，所述η为开挖面破坏率，所述AZ为开挖面破坏区域面积，单位为m2，所述A0为开挖面区断面面积，单位为m2，该区段中风化钙质板岩的极限等效塑性剪切应变为极限等效塑性拉应变为图7a为开挖第3步后(η＝0.109)，图7b开挖第6步后(η＝0.157)，图7c开挖第9步后(η＝0.201)，为不同开挖步下，开挖面的横断面与纵断面ZSI等值线图。图7a、图7b以及图7c中虚线区域ZSI<0，即破坏区。值得注意的是：单元破坏为剪切破坏还是拉伸破坏无法直接从等值线上看出，但从每个单元的状态或破坏形式可以通过FISH语言遍历单元后输出查看。

开挖面的剪切破坏区域主要集中在开挖面的底部，而拉伸破坏区域出现在开挖面中心附近。屈服区域则分布在弹性区域与破坏区域中间。开挖面前方破坏区域影响范围约为3.6m，并且随着开挖距离的增加变化并不明显，只有底部剪切破坏区域向前略有扩展。而开挖面的破坏率则是呈增加趋势的，并且屈服区域面积也在逐渐扩大，仍处在弹性区域的面积逐渐减小，说明随着开挖面的推进，扰动增加，隧道开挖面的稳定性在减弱，如图8所示。

为提高第13步开挖面的稳定性，减小破坏率，根据分析原因，在开挖12步的基础上，调整第13步开挖面的支护压力比λ，从图9中可以看出，支护压力比的增加，使得开挖面的破坏率逐渐呈负指数减小，拟合曲线方差为0.997。当λ＝0.6时，开挖面几乎全部破坏。若将破坏率控制在0.2左右，则应增加支护压力比，并且不应小于0.8。说明根据实际情况适当调整支护压力比对控制开挖面破坏率具有重要的作用。

本发明的有益效果为：

(1)结出单元状态指标ZSI表征单元安全性的评价方法，以此为基础引入渗透系数随体应变变化方程，渗透系数在屈服阶段和破坏阶段产生突跳，改进了FLAC3D渗流模拟中渗透性不变的缺陷。

(2)H-M耦合并考虑渗透性变化的监测点沉降值要远远大于不考虑渗流的沉降值，并且与实测结果较为吻合，说明富水区盾构隧道开挖时不能忽略H-M耦合的影响，单元的渗透性与单元状态是紧密相关的，不能按常数处理。

(3)计算结果显示，开挖面的剪切破坏区域主要集中在开挖面的底部，而拉伸破坏区域出现在开挖面中心附近。改变开挖面的支护压力比。随着支护压力比的增加，开挖面的破坏率逐渐减小，说明根据实际情况适当调整支护压力比是控制开挖面破坏率有效的途径。

扩充了岩土材料的安全评价体系，为流固耦合过程中渗透系数的计算提供了新的方法，改进了渗流数值模拟中渗透性不变的缺陷。使得计算结果更加准确、合理。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (4)

1.一种隧道流固耦合稳定性分析方法，其特征在于，包括：

根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，所述弹性阶段的状态指标包括剪切状态和拉伸状态；

根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标；

根据所述弹性阶段的状态指标、所述屈服阶段以及破坏阶段的状态指标确定所述材料单元耦合过程中渗透系数的变化值；

根据所述渗透系数进行隧道的流固耦合稳定性分析。

2.根据权利要求1所述的方法，其特征在于，所述根据材料单元实际的应力状态点计算所述材料单元弹性阶段的状态指标，包括：

将所述应力状态点进行受力分析；

根据空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离与该点在相同罗德角方向上到屈服面的距离之比确定所述状态指标为

其中，所述ZSI为状态指标，所述O′B为空间应力状态点相应的最稳定参考点沿最不利路径到屈服面的距离，所述O′P为该点在相同罗德角方向上到屈服面的距离，所述τ'_π为子午面GBO′面上B点纵坐标，所述τ_π为子午面GBO′面上P点纵坐标，所述I₁为应力张量第一不变量，所述为内摩擦角，所述θ_σ为应力罗德角，所述J₂为偏应力张量第二不变量，所述c为粘聚力。

3.根据权利要求2所述的方法，其特征在于，所述根据所述材料单元的破坏临界点的等效塑性应变值确定所述材料单元的屈服阶段和破坏阶段的状态指标，包括：

根据等效塑性剪切应变和极限等效塑性剪切应变确定所述材料单元在剪切屈服段内距离破坏段的接近程度；

根据所述接近程度确定屈服阶段和破坏阶段的状态指标为

$ZSI=min(1-\frac{{\mathrm{\&epsiv;}}^{ps}}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^{ps}},1-\frac{{\mathrm{\&epsiv;}}^{pt}}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^{pt}})---\left(2\right)$

其中，所述ε^ps为等效塑性剪切应变，所述为极限等效塑性剪切应变，所述ε^pt为等效塑性拉伸应变，所述为极限等效塑性拉伸应变。

4.根据权利要求2或3所述的方法，其特征在于，所述渗透系数的变化值为

$k=\left\{\begin{array}{cc}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI\&GreaterEqual;1\\ {\mathrm{\&xi;k}}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& 0\&le;ZSI<1\\ {\mathrm{\&xi;}}^{\&prime;}{k}_0\frac{{(1+{\mathrm{\&epsiv;}}_v/n_0)}^3}{1+{\mathrm{\&epsiv;}}_v}& ZSI<0\end{array}\right.---\left(3\right)$

其中，所述k为当前的渗透系数，所述k₀为初始渗透系数，所述ε_v为体应变，所述n₀为初始孔隙度，所述ξ为屈服阶段突跳系数，所述ξ'为破坏阶段突跳系数。