压缩试验模拟考虑的几个问题(本构模型和NMD算法)

1 引言

由于岩土材料问题的物理不稳定性，非线性材料的路径依赖性以及应力-应变响应的非线性特征，选择合适的本构模型会变得非常困难。《本构模型(Constitutive Models)选择》《IMASS---FLAC3D和3DEC新的本构模型(1)》《FLAC3D 7.0 新特性简介(P3)---新的本构模型》《FLAC2D---过去，现在和将来》。迄今为止，FLAC3D/3DEC已经内置了超过35种本构模型，当模拟一个问题时，我们不可能试验所有的本构模型。大多数情况下，总是从最简单的各向同性的弹性模型(Isotropic Elastic Model)入手。弹性模型(block zone cmodel assign elastic)只需要两个材料参数：体积模量和剪切模量(block zone prop  dens=0.0026 bulk=4000 shear=3000)，运行速度最快。在解决全尺寸的边值问题之前使用弹性模型测试能够对模拟的问题呈现出一个big picture，粗略判断出应力集中的位置，了解模型的预期响应，有助于重新定义网格的密度以及改变材料的本构模型。

对于一般的弹塑性问题，可以直接使用Mohr-Coulomb模型进行测试。Mohr-Coulomb模型是塑性模型组(Plastic Model Group)里最简单的模型。输入参数除了体积模量和剪切模量外，只需提供材料的密度，粘结力，内摩擦角和抗拉强度即可。

block zone property density 2.5E3 bulk 1.19E10 shear 1.1E10 friction 44 cohesion 2.72E5 tension 2E5

下面的例子演示了一个由Mohr-Coulomb材料组成的压缩试验。这个问题的物理意义是模拟一个矿柱的屈服行为。

2 模拟步骤

(1) 几何形状。几何形状必须在一定程度上代表着真实的物理问题，主要考虑以下几个方面：[1] 模型应加入多少细节(节理，断层，层面等)来表示地质结构，或者是否使用离散断裂网络DFN表示随机的岩体结构；[2] 模型边界的位置将如何影响模型结果？边界设置的太小，计算结果不准确，而边界设置的太大，计算时间将大大增加；[3] 如果使用可变形块体，在感兴趣的区域内，需要多大的网格密度才能得到精确的解答？这三个方面决定了合适模型的尺寸。block create 可以直接创建6种基元，它们分别是：polyhedron, brick, drum, prism, tunnel, tetrahedron, 在本例中使用drum产生一个圆柱体。

block create drum center-1 (0,0,0) center-2 (0,0,2) ...                     radius-1 1 radius-2 1 edges 16

(2) 划分网格。网格划分是一门艺术，没有万能的准则，在许多情况下可能需要不断调节尺寸以获得最优的结果。一些有限元软件，例如Rocscience RS2和MIDAS GTS NX(岩土工程有限元分析软件MIDAS GTS NX 2021 V1.1)中，确实内置了一些网格划分准则帮助用户进行网格划分，然而在FLAC3D/3DEC中没有这样显式的准则，不过提供了非常多的选项供用户自行选择，在本问题中，使用最大边长max-edge作为单元划分准则。

block zone generate-new max-edge 0.25

(3) 本构模型。如引言中所述，本构模型使用Mohr-Coulomb模型，同时设置该模型相应的材料参数。

block zone cmodel assign mohr-coulombblock zone property dens 2500 bulk 1.19e10 shear 1.1e10block zone property cohesion 2.72e5 friction 44 tension 2e5

(4) 边界条件。为了模拟一个矿柱的受力状态，在模型的顶部和底部施加数值相等的速度边界，类似于单轴抗压强度试验。

block gridpoint apply velocity (0, 0, 1e-3) range position-z 0 block gridpoint apply velocity (0, 0, -1e-3) range position-z 2

(5) 监测点布置。布置3个监测点。第一个点监测模型底部中心点的垂直位移，第二个点检测模型内部中点的垂直应力，第三个点检测模型侧表面中点的垂直应力。

block history displacement-z position (0,0,0)block history stress-zz position (0,0,1)block history stress-zz position (1,0,1)

3 模拟结果

下图所示的是在运行一定时步后(model step 20000)的位移矢量图和位移等值线图。

下图所示的是两个垂直应力监测点随矿柱顶板位移增加的变化。可以看出，矿柱内部的应力随着位移的增加而增加，而矿柱表面的应力增加不是非常显著。

4 塑性结果

四面体单元计算出的塑性结果不如六面体单元或高阶单元计算的塑性结果准确。如果需要非常精确的塑性计算，应该使用六面体单元。由于3DEC是四面体单元，因此为了得到精确的塑性计算结果，需要启用节点混合离散化(Nodal Mixed Discretization, NMD)算法。NMD算法在大多数情况下能够产生更精确的塑性计算结果，因此为四面体网格提供了一个更精确的解决方案。默认状态下，这个命令是关闭的(off), 需要使用下述命令打开:

block zone nodal-mixed-discretization on

比较结果显示，NMD on对最终位移的影响不大，但对应力有显著影响。NMD on使得在相同位移情况下塑性应力的值增大，因而这个NMD on导致计算结果和设计更保守一些。