岩土颗粒材料的连续离散耦合数值模拟

第1章连续离散耦合分析方法 连续离散耦合分析方法结合了有限单元法和离散单元法，将基于连续介质力学的有限单元法与基于接触检索、接触力计算和显式动力学求解的离散单元法结合在一起。一个典型的连续离散耦合分析可能包含数以千计的离散颗粒，每个颗粒被离散为单独的有限元网格。颗粒集合体的能量耗散包括：弹性滞后能、塑性应变能、颗粒破碎耗散能和颗粒间摩擦耗散能。颗粒在接触力和边界约束下的变形包括两个部分：转动和有限应变。有限应变与颗粒材料的力学性质有关，此外考虑材料的非线性，如断裂、破碎和磨损。颗粒材料由连续状态向非连续状态转化，导致颗粒形状逐渐变化、颗粒数目逐渐增多。 1.1 接触力模型 1.1.1 接触问题的一般描述 在连续离散耦合分析方法中，颗粒间的接触满足互不侵入条件并传递法向和切向接触力，如何处理颗粒间的接触是关键问题之一。从算法的角度来看，可以将与颗粒接触有关的问题分为两个方面：接触检索和接触力模型。接触检索是在颗粒集合体中找出发生接触的颗粒对，避免在相距较远的颗粒间进行接触力分析。一旦接触的颗粒对被检索到，就采用接触力模型计算颗粒间传递的接触力。颗粒间的接触通常是面接触，而非点接触，接触面的形状往往是不规则的，颗粒通过接触面上的节点传递法向和切向接触力。当颗粒间的法向接触力较小时，接触面的面积较小，而随着法向接触力的增加，颗粒表面粗糙的起伏发生弹塑性变形，导致接触面的面积逐渐增大。基于对颗粒形状、分布和表面粗糙度的假定，学者提出了颗粒间接触的理论模型或者微观力学模型。在计算力学中，通过变分形式进一步简化颗粒间接触的理论假定，并且认为法向接触力是法向侵入量的函数，而切向接触力是法向接触力和接触状态的函数。 将接触的边界值问题转化为在接触边界域Γ构造泛函Π及其变分形式，在接触边界域Γ，接触颗粒的位移场满足：C()u=0（1.1.1）式中：C是位移函数；u是接触边界域内的坐标。 接触问题的变分形式需要在接触边界域Γ上构造一个泛函，通过寻找泛函的驻值来满足不可贯入条件。根据构造泛函的形式不同，可分为*小二乘法、拉格朗日乘子法和罚函数法。为了接触边界域Γ上满足接触约束条件，*小二乘法定义泛函Π为 （1.1.2） 可以看出，当接触约束条件不能完全满足时，泛函始终为正值。当泛函Π取极小值时，即可得 （1.1.3） 拉格朗日乘子法引入接触约束条件的附加泛函Π1，与式（1.1.2）中的泛函相加，可得 （1.1.4） 式中：λ 是接触边界域Γ 上的独立函数，也称为拉格朗日乘子。 在有限元逼近中，通常采用基函数构造拉格朗日乘子，但会产生额外的未知变量，其物理意义是接触力的大小。在拉格朗日乘子法中，通过增加未知变量的个数来执行接触约束条件。在静态和隐式动力学问题中，通过求解线性代数方程组来求解未知变量，使接触约束条件可以严格满足。而在瞬时动力学问题中，由于没有形成和求解代数方程组，拉格朗日乘子法只能近似地满足。 为了消除和弥补拉格朗日乘子法的缺陷和不足，罚函数法引入接触约束条件的附加泛函Π2 ：（1.1.5） 式中：p是惩罚项。将附加泛函Π与式（1.1.2）中的泛函Π相加，得 （1.1.6） 由于 （1.1.7） 如果泛函Π在接触边界域Γ上为*小值，则p必须为正值。通过求解式（1.1.6）中修正泛函的极小值，近似满足接触约束条件。p越大，接触约束条件的满足程度越好，当p 无穷大时，接触约束条件能够精确满足。在静态或隐式动力学问题中，通过迭代求解的方法来精确满足不可贯入条件。而在瞬时动力学问题中，放弃完全不可贯入条件，而采用足够大的罚函数，使接触的侵入量相对于颗粒尺寸来说可以忽略不计。 在连续离散耦合分析方法中，基于颗粒的有限元网格离散并结合接触势的概念进行接触力分析。由于每个颗粒都被离散为单独的有限元网格，在接触力分析中，可以方便地使用有限元节点的几何坐标来描述接触颗粒的几何形状，并且接触面上接触力的分布更加真实。更重要的是，接触边界附近的局部应变场的数值畸变性大为改善，当考虑颗粒材料的断裂和破碎时，这一点尤为重要。 1.1.2 二维情况下的接触力计算 罚函数法认为接触的两个颗粒可以相互侵入，并通过侵入产生接触力。罚函数法中，接触泛函的标准形式为 （1.1.8） 式中：rt 和rc 分别是目标颗粒和接触颗粒位于接触面上点的位置矢量；Uc 是接触泛函；Γc 是接触区域。当p趋近于无穷大时，颗粒间几乎不会产生侵入： （1.1.9） 但在显式时步积分时，过大的惩罚项会大大减小*大稳定时间步长Δtmax ，因此在实际应用中，罚函数法通常采用足够大的惩罚项使颗粒间能产生一定的重叠量，而重叠量相对于颗粒尺寸来说又可以忽略不计。当采用罚函数法求解颗粒间的接触力时，根据接触力的计算策略不同，可以分为集中式接触力和分布式接触力，如图1.1.1 所示。集中式接触力认为节点的接触力是节点侵入目标颗粒的侵入量的函数，而分布式接触力则与接触颗粒和目标颗粒间重叠区域的形状及大小有关。 图1.1.1 分布式和集中式接触力 在本书的研究中，采用分布式接触力的定义，接触的两个颗粒分别表示为接触（contactor）颗粒和目标（target ）颗粒。当发生接触时，目标颗粒和接触颗粒相互重叠，被边界Γ围成的重叠区域面积为A，如图1.1.2 所示。 图1.1.2 重叠的微元面积dA所产生的接触力df[1] 当接触颗粒侵入目标颗粒的微元面积为dA时，其所产生微小的接触力df为 （1.1.10）式中：P是重叠区域内的一点；.是相应的势函数，下标c 和t 分别表示接触颗粒和目标颗粒。 式（1.1.10）定义的接触力可以看作接触颗粒的微元先侵入目标颗粒，而目标颗粒的微元随后侵入接触颗粒。对于接触中的颗粒，通过接触颗粒的势函数梯度来计算其接触力，因此接触力场是一个守恒场。如果接触颗粒上的点Pc 以任意路径AG B 侵入目标颗粒，在此 过程中由势函数梯度产生的接触力所做的功与路径AB无关，而只与起点和终点的位置有关。如果A点和B点都在目标颗粒的边界上，接触颗粒的点从目标颗粒的A点开始接触，运动到B点时结束接触。为了保证系统在A点和B点时的总能量相同，接触颗粒的点从A点运动到B点时，接触力所做的功为零，此时目标颗粒边界上的A点和B点处的势函数满足： （1.1.11） 反之，对于接触颗粒边界上的A点和B点有 （1.1.12） 将式（1.1.10）中微元dA产生的接触力df在重叠区域S上积分，即可得到总接触力： （1.1.13） 或是在重叠区域S的边界Γ上积分也可得到总接触力： （1.1.14） 式中：nΓ是边界Γ上的外向单位法向矢量。 在连续离散耦合分析中，每个颗粒都被离散为独立的有限元网格，因此颗粒可以表示为其所离散的有限单元组成的集合： （1.1.15） 式中：βc 和βt 分别是接触颗粒和目标颗粒；m和n分别是接触颗粒和目标颗粒的有限单元个数。 此时，颗粒的势函数可以表示为其所离散的所有单元的势函数之和： （1.1.16） 此时，接触力在重叠区域S上的积分表达式可以表示为 （1.1.17） 将有限单元域上的积分等价为单元边界上的积分，可得接触力为 （1.1.18） 将在重叠区域S的边界Γ上的积分转化为在有限单元边界上的积分之和，也就是说通过对重合区域内有限单元边界上的积分求和来计算颗粒间的接触力。在二维情况下，采用*简单的三角形常应变单元离散颗粒集合体，单元的边界是直线并且其几何信息可由三个节点的坐标表示，如图1.1.3 所示。 如本节所述，颗粒边界上的势函数φ应为常数以保证能量平衡，这里定义三角形单元中一点P的势函数φ为P （1.1.19） 式中：A(i=1，2，3)是如图1.1.3 所示的子三角形的面积，当点P位于三角形的形心时，φ(P)= 1 当点P位于三角形的边上时，当点P位于三角形以外时。 可以看出，势函数无条件满足颗粒边界上势函数为常数的条件。根据式（1.1.18），将两个三角形单元的接触问题转化为接触三角形与目标三角形的边之间的接触，以及目标三角形与接触三角形的边之间的接触，如图1.1.4所示。 图1.1.3 三角形常应变单元中一点P的势函数φ 图1.1.4 接触三角形与目标三角形的接触 图1.1.5 为目标三角形与接触三角形的边AB之间的接触，为了便于分析，定义局部坐标系(u， v)，通过式（1.1.20）将整体坐标系下的坐标转化到局部坐标系中： （1.1.20）式中：ri和rA是点i和点A在整体坐标系中的位置矢量；pi是点i在局部坐标系中的位置矢量。 图1.1.5 目标三角形与接触三角形的边AB接触所产生的分布式接触力 在局部坐标系中，确定边AB与目标三角形的特性交点并计算相应的势函数值。目标三角形施加在边AB上的总接触力为 （1.1.21） 式中：v是沿接触边界L上的积分操作；由于局部坐标系中矢量u和v不是单位矢量，故在式（1.1.21）中有u2。 如图1.1.5所示，在特征点之间势函数.线性变化，因此在边AB上的积分简化为求解边AB上势函数值围成的阴影部分的面积。将计算得到的总接触力转化为节点A和B上的等效节点力，采用相同的方法计算接触颗粒的其余边上的接触力。考虑由于目标颗粒侵入接触颗粒所产生的接触力，按照上面的方法依次计算接触颗粒施加在目标颗粒边上的接触力。*后更新接触颗粒和目标颗粒节点上的接触力。 1.1.3 三维情况下的接触力计算 在三维问题中，颗粒间接触力的计算与二维情况类似，通过将每个颗粒离散为独立的有限元网格，进而将颗粒内的势函数按有限单元分片。在连续离散耦合分析方法中，颗粒的形状不局限于规则的圆球，而可以是更接近真实颗粒的任意形状，如不规则多面体。为了减小算法实现的难度和提高计算效率，采用线性四面体单元（图1.1.6）离散颗粒。 基于接触力势概念的接触力模型能够保证颗粒接触处在有限侵入量时能保持能量和动量平衡。如本节所述，通过有限元离散网格可以方便地定义颗粒上的势函数.。本节给出了基于四面体单元离散网格的接触力模型，为了简化几何和数值计算上的困难，在每个四面体有限单元上定义势函数：