数值模拟方法

4．2．1　控制方程

4．2．1．1　有效应力原理

对于多孔介质材料，有效应力原理可写为：

式中：σ为总应力为有效应力；uw、ug分别为水压与气压；I为单位矩阵；χ因子与饱和度有关，完全饱和时取1。当不考虑大气边界条件以及忽略气压项对变形的影响时，上式变为：

分别考虑吸附水与自由水时，有效应力可分解为作用在土体骨架上的应力和作用在吸附水上的压力：

式中为土骨架上有效应力为作用有吸附水上的平均压力；nt为吸附水体积与总体积的比值。

以下给出下文中将涉及的几个参数定义。

孔隙率：

对于大变形问题，采用上标为0表示材料的初始构形，上式可改写为：

式中：

当用孔隙率时，有以下关系：

饱和度定义为：

自由水孔隙率为：

自由水的吸附水总的孔隙率为：

nf＝sn＋nt

4．2．1．2　多孔介质流固耦合控制方程

平衡方程的虚功形式为：

式中：δυ为虚速度场；δε＝sym 为虚变形率。体力项中包括水的重量fw和土骨架的重量f，fw＝（sn＋nt）ρwg。上式可改写为：

代入虚功方程，有：

4．2．1．3　多孔介质渗流连续方程

单元体内自由水和吸附水的质量为：

单位时间内单元体内水分质量变化率为：

单位时间内流入单元体的水分质量为

假定单位时间内流入单元体的水分等于其质量变化率，有：

由散度定理，上式改写为：

即有：

上式的虚功形式为：

上式在参考构形中可写为：

采用向后差分，上式可写为：

在当前构形中，上式写为：

4．2．2　本构关系

4．2．2．1　基本方程

根据弹塑性理论，总应变可分成弹性应变和塑性应变，其增量形式为：

弹性应变可应用广义虎克定律计算，其表达式为：

塑性应变可根据增量理论计算，需要确定材料的屈服函数、流动法则和硬化规律。

土可看作加工硬化材料，屈服面是应力和硬化参数的函数，其关系表示为：

塑性流动是由塑性势所引起的。塑性势面可用塑性函数g表示，它是应力状态的函数。塑性应变增量与应力｛σ｝之间的关系服从流动法则，其表达式为：

式（4－2－20）中dλ为比例常数。

若假定塑性势函数与屈服函数一致，即f（σij）＝g（σij），则为关联塑性流动法则。根据塑性应变增量理论，可得到弹塑性应力－应变的普遍关系式。

式（4－2－21）中［Dep］为弹塑性矩阵，可由下式表示：

式中：A＝F′。

4．2．2．2　Mohr－Coulomb模型

1）Mohr－Coulomb屈服准则

对于一般受力情况的岩土，其抗剪强度可用Mohr－Coulomb公式表示为：

式中：τn为极限抗剪强度；σn为受剪面上的法向应力，以拉为正。

Mohr－Coulomb抗剪强度曲线见图4－1。

图4－1　Mohr－Coulomb抗剪强度曲线

Mohr－Coulomb模型屈服面是一个不规则的六角形截面的角锥体表面，见图4－2，其在π平面上的投影如图4－3所示。

图4－2　Mohr－Coulomb模型屈服面

图4－3　π平面上的投影

式（4－2－22）还可采用应力不变量表示为以下形式：式（4－2－23）还可简化为：

式中：Rmc＝ sinθσsinφ ）。

Mohr－Coulomb模型屈服面在子午面内的投影可用Rmc·q－p应力平面表示，见图4－4。

图4－4　Rmc·q－p应力平面

2）流动法则

在子午面内，流动势函数采用双曲线函数表示（图4－5），在π平面内采用Menétrey和Willam（1995）提出的圆滑椭圆函数表示（图4－6），其表达式为：

图4－5　子午面内双曲线流动势函数

图4－6　π平面内Menétrey和Willam流动势函数

；φ为在Rmc·q－p应力平面内的膨胀角；c｜0为初始屈服应力；ε为系数，表示子午面内双曲线接近渐进线的程度；e为偏应力参数，一般取e＝

4．2．2．3　扩展Drucker－Prager模型

扩展Drucker－Prager模型（简称D－P模型）的屈服面在π平面上不是圆形的，屈服面在子午面上的形式有限性模型、双曲线模型和指数模型。扩展的D－P模型具有以下特点：

（1）适合用来模拟土体、岩石等摩擦材料，这些材料的屈服与围压有关，围压越大，材料的强度越高。

（2）适于用来模拟压缩屈服强度远大于拉伸屈服强度的材料。

（3）允许材料各向同性硬化或者软化。

（4）考虑了材料的剪胀性。

（5）可以模拟蠕变功能以描述材料的长期非弹性变形。

线性D－P模型屈服函数与流动势函数分别为：

式中：p，q，r分别为压力，第二与第三偏应力不变量；β，d分别为摩擦角与黏聚力；φ为剪胀角；K为材料三轴拉伸屈服应力与压缩应力比值；t表征屈服函数在子午π平面的形态，其表达式为：

4．2．2．4　混凝土损伤模型

该模型有以下特点：

（1）适用于混凝土和其他准脆性材料。

（2）各向同性损伤弹性定义，同时采用各向同性拉压塑性表征材料的非弹性性状。

（3）适用于单调荷载、循环荷载及低围压下的动力荷载。

（4）既可以用于素混凝土，也可用于钢筋混凝土材料。

模型应力应变关系、屈服函数及流动势函数分别为式（4－2－22）至式（4－2－24）。

其中：d为损伤变量，为等效塑性应变、温度及其他内变量的函数；¯p，¯q分别为有效压应力和有效第二偏应力不变量为最大有效主应力值分别为有效拉应力及有效压应力；σb0／σc0为初始压缩双轴屈服应力与单轴屈服应力的比值，默认值为1．16；Kc为初始屈服时，拉伸子午面上第二应力不变量与压缩子午面上第二应力不变量比值，0．5＜Kc≤1．0；φ为剪胀角；σt0为材料单轴拉伸强度；ξ为偏心率。

4．2．2．5　多孔介质渗流本构关系

多孔介质材料在低流速下的渗流本构关系采用达西定律，该定律认为通过单位时间内单位面积的流体与孔压的负梯度成比例时，即：

式中：s为饱和度；n为孔隙率；＾K为渗透系数；φ为水头高度，定义为φ＝z＋，在非饱和情况下，其渗透系数与饱和度及介质特性有关，表示为，＾K＝ksk；ks（s）是与饱和度有关的参数，ks（1）＝1．0，k是多孔介质完全饱和时的渗透系数。

4．2．3　材料开裂模拟方法

可以采用考虑孔压的黏性单元（Cohesive Element）模拟在孔压作用下桩土界面材料的开裂性状。黏性单元将材料断裂能与损伤结合在一起，能模拟在孔压作用下材料的损伤起裂及裂纹扩展过程，在研究材料黏结性能、岩石破裂等方面得到了较为广泛的应用。

4．2．3．1　界面本构关系

桩土界面特性对基桩承载性能有重要影响，与计算模型尺寸相比，界面物理厚度可以忽略不计，采用基于力－位移关系本构模型（Seperation－Traction），该模型有以下特点：

（1）可以描述复合材料开裂过程。

（2）断裂能可以假定为界面法向与切向变形率的函数。

（3）材料损伤前是线性弹性关系。

（4）材料破坏是由于材料刚度逐渐弱化引起。

（5）可以定义多个损伤机制。

线弹性关系为：

式中：t为应力矢量；δn，δs，δt分别为法向与切向开度；T0为黏性单元的初始厚度。

4．2．3．2　起裂准则

起裂准则包括最大应力、最大应变及二次名义应力等，本书采用最大名义应力准则，当满足式（4－2－34）时材料起裂。

4．2．3．3　损伤演化准则

1）有效应力

采用标量损伤因子D代表材料总有效损伤，该有效损伤是多种损伤机制的综合反应。起裂后，在加载过程中D值由0增加到1，在此过程中应力矢量计算方法为：

式中，，分别为未损伤状态下弹性本构关系得到的应力矢量。

2）有效位移

为了描述法向与切向应力耦合作用下的损伤演化过程，采用有效位移，定义为：

3）混合模式（Mixed－Mode）

在法向与切向应力共同作用下存在不同的损伤机理，有必要对此进行区分。通过定义混合模式可以定量法向与切向变形各自的比例关系。可以采用两种方式定义混合模式：基于能量的方法与基于力（Traction）的方法。基于能量的方法表达式为：

式中：GT＝Gn＋Gs＋Gt，Gn、Gs、Gt分别是法向力、第一切向力、第二切向力与其共轭的位移所做的功。显然，剪切力所做功为m2＋m3（或1－m1）。数值计算时，需据此定义材料损伤演化参数。基于力的混合模式定义见图4－7、图4－8。

图4－7　混合模式

图4－8　混合受力模式对材料损伤破坏的影响

式中：τ＝为有效剪切应力。要注意的是，上述两种定义方法在某些情况下是不同的。比如，材料承受切向与法向力，但只有法向变形。据能量方法切向功为零，不存在混合模式，但据力的方法，则存在混合模式。

损伤演化模型中两个参数需定义，一个是材料完全破坏时的有效位移δfm（或者是完全破坏所需断裂能Gc），参见图4－9。另外一个参数是损伤标量D，该值可以有多种定义方法：线性、指数的函数形式或者是直接定义为有效位移的插值形式。一般情况下，两个参数均是混合模式、温度或其他场变量的函数。

图4－9　软化模型

本书采用线性软化模型，损伤因子采用下式计算：

式中：Gc、分别为总断裂能与起裂时的有效应力；为加载过程中的最大有效位移。当满足下式时，材料完全破坏：

4．2．3．4　渗流关系

采用考虑孔压的黏性单元可以保持渗流的连续性，可以模拟孔压致裂现象。流体本构关系包括沿黏性单元内间隙的切向与法向渗流定义，如图4－10所示。

图4－10　切向与法向渗流模型

1）切向渗流

对于牛顿流体，有以下关系式：

式中：Q为切割流速；p为沿黏性单元孔压梯度；d为黏性单元开度，有d＝dcurr－dorig＋dinit；Kt为渗透系数，Kt＝；μ为流体的黏性。

2）法向渗流

图4－11中，pt、pb分别为单元顶部与底部界面孔压值，pi为单元中部孔压值。

3）黏性单元组合渗流关系（表4－1）

图4－11　法向渗流计算模型

表4－1　孔压单元渗流定义

4．2．4　接触面计算模型

界面接触问题采用约束法，其接触算法需要解决以下问题：①确定接触区域；②判断界面接触状态；③界面不同接触状态下的力学行为。

接触界面有许多定义方法，本处使用点面接触单元。典型的三结点接触单元如图4－12所示。其中，I－J为目标面，N为从节点。在外载作用下，从结点相对于目标面发生滑动，在移动过程中从节点与不同的目标面发生接触或脱离，从而模拟界面受力性状。

图4－12　接触计算示意图
（a）三结点接触单元图；（b）MC本构关系图

接触界面本构关系采用弹性库仑摩擦模型，界面剪应力和法向应力关系如图4－12所示。剪切应力与法向应力的函数关系为：式中：τ为剪应力；p为法向应力；ω为界面相对位移；μ为界面摩擦系数；ωs为弹性极限相对位移。