CN107169236B - 一种基于有限元与离散元耦合的虚拟三轴试验仿真方法 - Google Patents

Abstract

本发明公开了一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，首先是建立常规离散元虚拟三轴试样颗粒模型，然后借助离散元应力伺服原理，将试样伺服到初始固结应力状态，其次将侧向刚性墙体模型替换为有限元单元模型，建立有限元与离散元耦合虚拟三轴试验模型，最后施加试验围压与轴向加载速度，开展耦合虚拟三轴试验。与现有技术相比，本发明提供的虚拟三轴试验仿真方法能够在试验中实现侧向围压的柔性施加，保证压缩过程中试样侧向可以自由变形，可以真实地模拟室内三轴试验。通过虚拟数值试验仿真实例，验证了本发明的基于有限元与离散元耦合的虚拟三轴试验仿真方法的优势性。

Description

一种基于有限元与离散元耦合的虚拟三轴试验仿真方法

技术领域

本发明涉及一种虚拟三轴试验的仿真方法，具体地说，涉及一种基于有限元与离散元耦合的虚拟三轴试验的仿真方法。

本发明主要应用于岩体工程领域颗粒材料虚拟试验仿真方面，可以实现对颗粒材料室内三轴试验的计算机真实虚拟仿真模拟。

背景技术

近些年来随着计算机和数值计算方法的发展，虚拟试验仿真技术在岩土工程领域得到了广泛的应用。虚拟数值试验作为室内试验的一种计算机虚拟仿真，其在散体材料介质复杂力学行为研究方面具有得天独厚的优势。

颗粒离散元虚拟试验作为一种常用的虚拟仿真试验，其核心原理是基于非连续介质理论的数值计算方法，特别适用于颗粒材料这类非连续介质的模拟，尤其是对材料介质细观力学机理的研究。

目前，常规颗粒离散元虚拟三轴试验主要是采用无摩擦刚性墙体对试样施加侧向围压，这种围压加载方式属于刚性加载，试验压缩过程中限制了试样侧向的自由变形，这与室内三轴试验的柔性加载方式不同(即，通过试样外层橡胶皮套来对试样施加侧向围压，试验中允许试样侧向自由变形)，导致虚拟试验仿真结果与实际结果存在一定的偏差，故这种加载方式的虚拟三轴试验无法实现对室内三轴试验的真实仿真模拟。

发明内容

为了克服常规颗粒离散元虚拟三轴试验中无法施加侧向柔性应力，无法实现对室内三轴试验的真实仿真模拟，本发明的目的是提供一种基于有限元与离散元耦合的虚拟三轴试验仿真方法。该方法在试验中能对试样侧向围压进行柔性加载，可实现对室内三轴试验的真实仿真模拟。

为实现上述目的，本发明采用以下技术方案：一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，包括如下步骤：

(1)根据虚拟三轴试样尺寸，建立常规离散元三轴试验模型；

(2)利用应力伺服控制原理，将虚拟试样伺服到初始固结应力状态；

(3)利用有限元模型替换侧向墙体模型，建立有限元与离散元耦合的虚拟三轴试验模型；

(4)施加侧向围压和轴向加载速度，开展有限元与离散元耦合虚拟三轴试验数值仿真模拟。

其中，步骤(1)的具体步骤包括：

(1.1)根据虚拟三轴试样尺寸，建立无摩擦刚性墙体模型；

(1.2)根据试样颗粒级配要求，建立虚拟试样球体颗粒模型。

步骤(2)的具体步骤包括：

(2.1)计算当前时步侧向墙体和上下墙体的平均应力值，计算公式如下：

式中，S_r和S_h分别为当前时步球体颗粒作用于侧向和上下墙体的平均应力，

分别为当前时步球体颗粒作用于侧向、上部、下部墙体上的作用力，D_t为当前时步侧向墙体的直径，H_t为当前时步上下墙体间的高度；

(2.2)计算下一时步墙体的移动速度，计算公式如下：

式中，v_r为侧向圆柱形墙体径向移动速度，

和

为侧向和上下墙体的初始固结应力，G为增益系数，按如下公式计算：

式中，

为球体颗粒与侧向墙体间的法向接触刚度，N_r为与侧向墙体接触的球体颗粒数目，N_h为与上下墙体接触的球体颗粒总数，△t为当前计算时步，α为松弛系数，取值范围为(0,1]，通常可取0.5；

(2.3)离散元循环迭代一个时步，更新侧向圆柱形墙体直径D_t和上下墙体间的高度H_t以及球体颗粒位置；

(2.4)重复执行上述步骤(2.1)至(2.3)，直到当前时步侧向墙体平均应力值S_r和上下墙体平均应力值S_h与其初始固结应力值

的误差小于指定容差ε时，计算迭代终止，最终获得试样初始固结应力状态，其中，ε按如下公式进行计算：

式中，ε取值范围为[0.001,0.005]。

步骤(3)的具体步骤包括：

(3.1)删除侧向墙体模型，在原墙体模型位置处建立一层有限元单元模型替代墙体模型；

(3.2)根据建立的有限元单元模型，在单元内侧建立耦合控制墙体单元；

(3.3)检索与存储每一控制墙体顶点与有限元模型单元节点的对应关系。

步骤(4)的具体步骤包括：

(4.1)在有限元单元模型外侧边界施加三轴试验侧向围压；

(4.2)设置上下墙体加载板的轴向加载速度；

(4.3)在动力求解模式下，将有限元和离散元设置为相同的分析时步；

(4.4)开始耦合虚拟三轴试验仿真模拟，每一时步有限元与离散元彼此间进行数据交换，并在每隔一定时步记录试验的轴向偏应力与轴向应变；

(4.5)当轴向应变达到指定目标值，停止试验加载，保存试验计算结果。

更进一步，步骤(4.4)中有限元与离散元耦合虚拟三轴试验是基于有限元(FEM)和离散元(DEM)相互耦合计算实现的，其中FEM是作为服务器(server)，DEM作为客户端(client)，两者耦合分析同步进行，在每一步耦合计算过程中，FEM先计算循环一步，更新单元节点位置，并将耦合节点(有限元单元内侧节点)位置发送给DEM，随后，DEM更新控制墙体的位置并循环计算一步，更新颗粒的位置，接着将颗粒作用于控制墙体上的作用力发送给FEM，通过作用力映射原理获得控制墙体对应耦合节点上的等效节点力；接着FEM再循环计算一步，更新单元节点位置，将耦合节点速度再次发送给DEM，DEM再次更新控制墙体的位置并循环计算一步，更新离散域内颗粒的位置，接着DEM将控制墙体的作用力再次发送给FEM；依次类推，耦合逐步计算下去。

更进一步的，步骤(4.4)中作用力映射原理核心是将球体颗粒作用于控制墙体上的力转为对应的耦合节点上等效节点力；

假设某个球体颗粒作用于控制墙体上的作用力为(F_x，F_y，F_z)，作用点在控制墙体上的位置为(x_c，y_c，z_c)，相应地，作用力等效到控制墙体所对应的三个耦合节点上的等效节点力可由如下公式计算：

式中，(F_xi，F_yi，F_zi)为作用力等效到控制墙体第i个顶点的等效力，w_i为第i个顶点面积权重系数，可由下式计算：

式中，A为控制墙体面积，A_i为作用点与控制墙体三个顶点组成的子三角形，满足

有益效果：与现有技术相比，本发明基于有限元与离散元耦合的虚拟三轴试验仿真方法可以实现侧向围压的柔性施加，能够真实地模拟室内三轴试验。本发明方法原理简单、数值计算效率高且稳定，适用于散体材料(如，土体、土石混合料、堆石料)宏细观力学特性的分析研究。

附图说明

图1是本发明基于有限元与离散元耦合的虚拟三轴试验仿真方法流程图；

图2是本发明实施例虚拟试样模型示意图；

图3A是根据图2建立的离散元虚拟三轴试验墙体模型图；

图3B是根据图2建立的离散元虚拟三轴试样球体颗粒模型图；

图4是图3A、图3B应力伺服后的虚拟三轴试验模型图；

图5是根据图4建立的有限元单元模型图；

图6是沿图5中有限元单元边建立的控制墙体模型图；

图7是虚拟三轴试验加载条件示意图；

图8是虚拟三轴试验有限元与离散元耦合计算示意图；

图9是球体作用力等效节点力分配示意图；

图10是虚拟三轴试验获得的轴向应变与轴向偏应力曲线图；

图11A是三轴试验后试样颗粒位移云图；

图11B是三轴试验后有限元单元位移云图。

具体实施方式

下面结合附图和具体实施例进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，本发明公开了一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，包括以下步骤：

S1：根据虚拟三轴试样尺寸，建立常规离散元三轴试验模型。具体的步骤为：

S1.1：根据虚拟的三轴试样尺寸，建立无摩擦刚性墙体模型；

图2所示为虚拟试样几何模型，其直径为D，高度为H。根据虚拟试样尺寸建立三轴试验墙体模型，如图3A所示，墙体模型体是由侧向圆柱形墙体和上下两个平板墙体组成，其中圆柱形墙体的直径即为试样的直径D，上下墙体间的相对距离即为试样的高度H。

S1.2：根据试样颗粒级配要求，建立虚拟试样球体颗粒模型；

虚拟试样中的颗粒是采用不同大小球体单元模拟，球体的大小分布满足规定的颗粒级配要求。如图3B所示为建立的虚拟三轴试样球体颗粒模型，模型中球体颗粒彼此相互接触，且与墙体模型也相互接触。

S2：利用应力伺服控制原理，将虚拟试样伺服到初始固结应力状态。具体的步骤为：

S2.1：根据公式(1)计算当前时步侧向墙体和上下墙体的平均应力值；

式中，S_r和S_h分别为当前时步球体颗粒作用于侧向和上下墙体的平均应力，

分别为当前时步球体颗粒作用于侧向、上部、下部墙体上的作用力，D_t为当前时步侧向墙体的直径，H_t为当前时步上下墙体间的高度；

S2.2：根据公式(2)计算下一时步墙体的移动速度；

式中，v_r为侧向圆柱形墙体径向移动速度，

和

为侧向和上下墙体的初始固结应力，G为增益系数，按如下公式计算：

式中，

为球体颗粒与侧向墙体间的法向接触刚度，N_r为与侧向墙体接触的球体颗粒数目，N_h为与上下墙体接触的球体颗粒总数，△t为当前计算时步，α为松弛系数，通常可取0.5；

S2.3：离散元循环迭代一个时步，更新侧向圆柱形墙体直径D_t和上下墙体间的高度H_t以及球体颗粒位置；

S2.4：重复执行上述步骤S2.1至步骤S2.3，直到当前时步侧向墙体平均应力值S_r和上下墙体平均应力值S_h与其初始固结应力值

的误差小于指定容差ε(通常指定容差取为0.001-0.005)时，计算迭代终止，最终获得试样初始固结应力状态。

在离散元中，应力伺服的原理是通过反复调整墙体的位置使试样内部应力逐步逼近目标应力值。如图3A、图3B所示，通过步骤S2将图2中初始虚拟试样伺服到指定的固结应力状态(试验围压)，此过程相当于室内试验的固结过程；在应力伺服后，虚拟试样发生了变化，其中直径变为D′，高度变为H′，如图4所示。

通常，当前时步侧向墙体平均应力值和上下墙体平均应力值与其初始固结应力值的误差小于指定容差值ε时，最终获得试样初始固结应力状态，其中，ε按如下公式进行计算：：

式中，ε取值范围为[0.001,0.005]。

S3：利用有限元模型替换侧向墙体模型，建立有限元与离散元耦合的虚拟三轴试验模型。具体的步骤为：

S3.1：删除侧向墙体模型，并在原墙体模型位置处建立一层有限元单元模型替代墙体模型；

如图5所示，删除图4中侧向圆柱形墙体，在原墙体位置建立一层有限元单元模型，模型单元内侧直径与圆柱形墙体模型的直径相同，且每个单元均为一个规则六面体单元。

S3.2：根据建立的有限元单元模型，在单元内侧建立耦合控制墙体单元；

如图6所示，在图5中单元模型内侧边界建立耦合控制墙体单元；每个单元内侧四边形均可拆分建立两个三角形的耦合控制墙体，且控制墙体顶点与单元顶点相互重合。

S3.3：检索与存储每一控制墙体顶点与有限元模型单元节点的对应关系。

如图6所示，根据控制墙体与有限元单元边的建立关系，确定和存储每个控制墙体顶点与有限元模型单元节点的对应关系，用于后续步骤S4中耦合数据的传递与交换。

S4：施加侧向围压和轴向加载速度，开展有限元与离散元耦合虚拟三轴试验数值仿真模拟。具体步骤为：

S4.1：在有限元单元模型外侧边界施加三轴试验侧向围压；

本发明中，虚拟三轴试验的围压是通过试样外层有限元模型间接施加到试样上，这与室内三轴试样外层的橡胶皮套功能类似。如图7所示，在有限元单元模型外侧边界施加均匀的法向应力，模拟室内试验的侧向围压。

S4.2：设置上下墙体加载板的轴向加载速度；

如图8所示，虚拟三轴试验是通过速度加载方式进行应力加载，即对上部墙体施加一个恒定的向下的移动速度，对下部墙体施加一个相同大小恒定的向上的移动速度；同时将有限元单元底部和顶部节点的速度设置为与相应墙体相同的移动速度，即上部顶端节点的速度与上部墙体的移动速度一致，下部底端节点的速度与下部墙体移动速度一致。

S4.3：在动力求解模式下，将有限元和离散元设置为相同的分析时步；

本发明中，为保证有限元与离散元耦合分析同步进行，需将有限元和离散元分析时步设置为同一时步大小。

S4.4：开始耦合虚拟三轴试验仿真模拟，每一时步有限元与离散元彼此间进行数据交换，并在每隔一定时步记录试验的轴向偏应力与轴向应变；

本发明中，有限元与离散元耦合分析是同步进行的，即有限元计算一个时步后，紧接着离散元计算一个时步，每一个计算时步中进行数据间的相互交换，整个耦合计算过程见图8所示。

如图8所示，有限元与离散元耦合虚拟三轴试验是基于有限元(FEM)和离散元(DEM)相互耦合计算实现的，其中FEM是作为服务器(server)，DEM作为客户端(client)，两者耦合分析同步进行，在每一步耦合计算过程中，FEM先计算循环一步，更新单元节点位置，并将耦合节点(有限元单元内侧节点)位置发送给DEM，随后，DEM更新控制墙体的位置并循环计算一步，更新颗粒的位置，接着将颗粒作用于控制墙体上的作用力发送给FEM，通过作用力映射原理获得控制墙体对应耦合节点上的等效节点力；接着FEM再循环计算一步，更新单元节点位置，将耦合节点速度再次发送给DEM，DEM再次更新控制墙体的位置并循环计算一步，更新离散域内颗粒的位置，接着DEM将控制墙体的作用力再次发送给FEM；依次类推，耦合逐步计算下去。

在计算过程中，作用力映射原理核心是将球体颗粒作用于控制墙体上的力转为对应的耦合节点上等效节点力；

假设某个球体颗粒作用于控制墙体上的作用力为(F_x，F_y，F_z)，作用点在控制墙体上的位置为(x_c，y_c，z_c)，相应地，作用力等效到控制墙体所对应的三个耦合节点上的等效节点力可由如下公式计算：

式中，(F_xi，F_yi，F_zi)为作用力等效到控制墙体第i个顶点的等效力，w_i为第i个顶点面积权重系数，可由下式计算：

式中，A为控制墙体面积，A_i为作用点与控制墙体三个顶点组成的子三角形，满足

图9所示为耦合过程节点力分配权重系数计算原理，对于控制墙体每个顶点权重系数w_i＝A_i/A。

S4.5：当轴向应变达到指定目标值，停止试验加载，保存试验结果。

下面以土体虚拟三轴试验为例，说明本发明的基于有限元与离散元耦合的虚拟三轴试验仿真方法的具体实施过程：

1、根据虚拟三轴试样尺寸，建立常规离散元三轴试验模型。

实施例中选取的虚拟试样尺寸为5cm×10cm(D×H)，见图2所示。图3A、图3B分别所示为根据虚拟试样尺寸建立的三轴试验墙体模型和试样颗粒模型，其中试样颗粒模型由12849个球体组成，球体的最小半径为1.0mm，最大半径为1.6mm，平均半径为1.4mm。

2、将虚拟试样伺服到指定的初始固结应力状态。

实施例中三轴试验的围压选取为500.0kpa，利用应力伺服原理将虚拟试样伺服到初始围压固结状态，同时取应力伺服终止的指定容差ε＝0.005，最终伺服结束后，球体颗粒作用于侧向墙体和上下墙体的平均应力值分别为499.5kpa和497.6kpa。如图4所示，伺服后的虚拟试样尺寸变为4.68cm×9.63cm(D′×H′)。

3、建立有限元与离散元耦合虚拟三轴试验模型。

如图5所示，在应力伺服后，删除侧向圆柱形墙体模型，并在相应位置建立一层有限元单元模型，单元总数为800，所有单元大小一致，且单元内侧边界与原墙体模型位置重合；如图6所示，根据建立的有限元单元模型，沿单元内侧边界建立一组耦合控制墙体，且墙体节点与单元节点在位置上是重叠，存储单元节点与控制墙体顶点的对应关系。

4、开展有限元与离散元耦合虚拟三轴试验仿真模拟。

首先，如图7中所示的围压和速度的施加方式，在有限元单元模型外侧边界施加500kpa围压，上下墙体竖向移动速度v设为0.0005m/s；其次将有限元和离散元设置为相同大小的分析时步，时步取为1.0e-5s；最后根据图8所的耦合计算原理及图9所示节点力分配原理进行虚拟三轴试验耦合模拟，在模拟过程中记录轴向应变与轴向偏应力，当轴向应变达到12％时，停止加载，保存计算结果。图10为虚拟三轴试验过程记录的轴向应变与轴向偏应力曲线图。图11A、图11B分别为三轴试验后试样颗粒和有限元单元变形云图，由图11B中看出，虚拟试验压缩过程中试样有限元单元的侧向变形是不一致，实现了室内三轴试验侧向围压柔性加载的目的，达到了本发明的最终目标，且由图11A中看出，试样压缩后发生了鼓胀破坏，这与试样室内三轴试验破坏形式相一致，也进一步验证了本发明的优势。

Claims (3)

1.一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，其特征在于：它包括如下步骤：

(1)根据虚拟三轴试样尺寸，建立常规离散元三轴试验模型；

(2)利用应力伺服控制原理，将虚拟试样伺服到初始固结应力状态；

(3)利用有限元模型替换侧向墙体模型，建立有限元与离散元耦合虚拟三轴试验模型，即；

(3.1)删除侧向墙体模型，在原墙体模型位置处建立一层有限元单元模型；

(3.2)根据建立的有限元单元模型，在单元内侧建立耦合控制墙体单元；

(3.3)检索与存储每一控制墙体顶点与有限元模型单元节点的对应关系；

(4)施加侧向围压和轴向加载速度，开展有限元与离散元耦合虚拟三轴试验数值仿真模拟，即：

(4.1)在有限元单元模型外侧边界施加三轴试验侧向围压；

(4.2)设置上下墙体加载板的轴向加载速度；

(4.3)在动力求解模式下，将有限元和离散元设置为相同的分析时步；

(4.4)开始耦合虚拟三轴试验仿真模拟，每一时步有限元与离散元彼此间进行数据交换，并在每隔一定时步记录试验的轴向偏应力与轴向应变；

(4.5)当轴向应变达到指定目标值，停止试验加载，保存试验计算结果。

2.根据权利要求1所述的一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，其特征在于，所述步骤(1)根据虚拟三轴试样尺寸，建立常规离散元三轴试验模型的具体步骤包括：

(1.1)根据虚拟三轴试样尺寸，建立无摩擦刚性墙体模型；

(1.2)根据试样颗粒级配要求，建立虚拟试样球体颗粒模型。

3.根据权利要求1所述的一种基于有限元与离散元耦合的虚拟三轴试验仿真方法，其特征在于，所述步骤(2)利用应力伺服控制原理，将虚拟试样伺服到初始固结应力状态的具体步骤包括：

(2.1)根据公式(1)计算当前时步侧向墙体和上下墙体的平均应力值；

式中，S_r和S_h分别为当前时步球体颗粒作用于侧向和上下墙体的平均应力，

分别为当前时步球体颗粒作用于侧向、上部、下部墙体上的作用力，D_t为当前时步侧向墙体的直径，H_t为当前时步上下墙体间的高度；

(2.2)根据公式(2)计算下一时步墙体的移动速度；

式中，v_r为侧向圆柱形墙体径向移动速度，

和

为侧向和上下墙体的初始固结应力，G为增益系数，按如下公式计算：

式中，

为球体颗粒与侧向墙体间的法向接触刚度，N_r为与侧向墙体接触的球体颗粒数目，N_h为与上下墙体接触的球体颗粒总数，△t为当前计算时步，α为松弛系数，α＝0.5；

(2.3)离散元循环迭代一个时步，更新侧向圆柱形墙体直径D_t和上下墙体间的高度H_t以及球体颗粒位置；

(2.4)重复执行上述步骤(2.1)至步骤(2.3)，直到当前时步侧向墙体平均应力值S_r和上下墙体平均应力值S_h与其初始固结应力值

的误差小于指定容差ε时，计算迭代终止，最终获得试样初始固结应力状态，其中，ε按如下公式进行计算：

式中，ε取值范围为[0.001,0.005]。

Images