[发明专利]力学渗流耦合的非饱和土数值模拟方法及系统

权利要求书

1.一种力学渗流耦合的非饱和土数值模拟方法，其特征在于，包括如下步骤：

S100：根据非饱和土的地质模型建立有限元数值模型；

S200：根据非饱和土的物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；

S300：根据有限元数值模型、参数特征、动静态模拟要求确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；

S400：以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。

2.根据权利要求1所述的力学渗流耦合的非饱和土数值模拟方法，其特征在于，以均匀时间步长递增时间的方式开展有限元循环计算，获取分析的计算步骤包括

S410：施加当前时间步的力学、孔隙流体的边界荷载增量；

S420：对有限单元计算应变增量、更新单元内的孔隙流体压力、流速；

S430：在单元内实施BBM一致性算法的计算，更新单元应力，计算单元一致性切线模量；

S440：计算单元节点的广义内力以及节点的广义不平衡力；

S450：判断节点的广义不平衡力是否满足平衡要求，如满足平衡要求则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算；

S460：计算单元刚度阵，由单元刚度阵累加获得总体刚度阵，由总体刚度阵和节点的广义不平衡力计算更新节点广义位移；

S470：判断节点的广义位移是否满足收敛条件，如满足则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算；

S480：跳出迭代循环至下一时间步时相应获得当前时间下模型的位移、应力、应变、孔隙流体压力、流速；

S490：判断是否达到总的计算时间步，未达到则开展新的时间步计算，达到则跳出有限元计算并获取相应的分析结果。

3.根据权利要求2所述的力学渗流耦合的非饱和土数值模拟方法，其特征在于，BBM一致性算法的计算步骤为：

基于Cosserat连续体理论，对非饱和土的固体骨架的变形计算，其运动微分方程可表示为：

式1中u_i和ω_i分别为固体骨架的位移和微转角，σ_jk和μ_jk分别为传统应力和偶应力，b_i和g_i分别为体力和体力偶，ρ和I_c分别是密度和微惯性矩，l_c为内尺度参数，δ_ik和e_ijk分别为克罗内克尔张量和置换张量，表示梯度算子；

考虑非饱和土中孔隙水和孔隙气体在渗流过程中满足的质量守恒方程和动量守恒方程，孔隙水和孔隙气体满足的控制方程为

其中,

式2～6中，下标w和a分别表示孔隙水和气，S_w、S_a是饱和度，ρ_w、ρ_a表示密度，表示达西速度，n为孔隙率；

针对单元内的应力更新计算，考虑在一个时间增量步上的塑性变形，当前应力表示为：

式7中，Δσ^e为当前时间步应力增量的弹性预测值，为应力作用引起的塑性应变增量，Δε_s为吸力引起的应变增量，D_e为线弹性本构模量矩阵；

应力增量的弹性预测值表示为

Δσ^e＝D_eΔε -式8

式8中，Δε为当前时间步总的应变增量；

应用关联塑性流动准则，塑性应变增量则表示为

式9中Δλ为塑性乘子，s为偏应力向量；

基于牛顿-拉弗森方法迭代计算当前时间步的塑性乘子Δλ，考虑在从第v步迭代到第v+1步时，屈服面方程应得到满足，即

由此，当前时间步的塑性乘子表示为

Δλ_v+1＝Δλ_v+δ(Δλ)_v -式12

式10～12中，下标v和v+1分别代表第v步和第v+1步的变量，屈服面偏导数F_,λ表示为

式13中，F_,p为加载湿陷屈服面对等效塑性的偏导数；

当前时间步应力向量的计算式为

σ_v+1＝σ_v-(2G^*s+b₂m)|_vδ(Δλ)_v -式14

当前等效塑性应变的更新表示为

其中，修正的剪切模量和体积模量分别表示为

其中，相关参数表示为

推导得到的应力对应变和吸力偏导量的一致性弹塑性切线模量矩阵和表示如下：

式22～24中，