重要通知：本公众号已开通赞赏功能，各位读者老爷有钱的可以给呱太捧个钱场！推荐使用微信PC端进行阅读，以避免命令流显示畸变。文章所附命令流均可直接复制到FLAC3D 6.0中运行，若有报错，请手动重输命令流。如果我的推文对你有帮助，请关注我，并点赞，在看，分享。======正文======在FLAC3D中，施加与坐标轴成线性关系的项目时，需要使用gradient关键字进行施加。对于线性项，其位于(x,y,z)处的项值计算公式如下：S = S(0) + grad(x)*x + grad(y)*y+ grad(z)*z        (1)式中S(0)为该项在全局坐标系下原点处的值；grad(i)分别为x/y/z方向的梯度。

一般遇到的线性项仅与z轴成线性，因此本文仅讨论线性项与z轴成线性关系的情况，此时，计算式位于(x,y,z)处的项值计算公式变为：

S = S(0) + grad(z)*z                       (2)以初始化浅埋工况下的自重应力为例，示意见图1，施加的命令流为：

zone initialize stress-zz S(0) grad 0 0 grad(z)zone initialize stress-xx μ*S(0) grad 0 0 μ*grad(z);μ为侧压力系数zone initialize stress-yy μ*S(0) grad 0 0 μ*grad(z)

图1 浅埋自重应力示意图

由此可见，在一般的梯度施加中，只需要确定S(0)与grad(z)即可，梯度计算示意图见图2：

图2 梯度计算示意图由图可知，在计算梯度时，首先需要确定线性荷载两端的两个坐标值；以图1的竖向应力为例，易知位于模型顶端坐标为(0,0)，底端坐标为(z,ρgz)。确定了坐标之后即可求得该线性项关于z的方程：由几何关系易知，梯度值grad(z)为方程的斜率，令z = 0 即可求得S(0)。======案例演示======

以初始化浅埋工况的自重应力为例(侧压力系数取0.5)，模型尺寸为10m*10m*10m，模型密度取为1000kg/m^3，重力加速度取为10m/s^2。分两个工况模拟，工况1的底面z坐标为0，工况2的底面z坐标为-3。

工况1：即z2=0，z1=10；模型顶部的自重应力为S1=0，底部的自重应力为S2=-1e5Pa。由式(3)计算得到梯度grad = (-1e5-0)/(0-10)=1e4Pa/m；S(0) = -1e5。因此工况1的命令流如下：

;模型底面z坐标为0model newzone create brick size 10 10 10zone cmodel assign mohr-coulombzone property bulk 3e8 shear 1e8 cohesion 1e6 friction 15 zone initialize density 1000zone initialize stress-zz -1e5 grad 0 0 1e4zone initialize stress-xx -5e4 grad 0 0 5e3zone initialize stress-yy -5e4 grad 0 0 5e3

运行上述命令流后，模型的三个正应力云图如下：

图3 stress-zz云图(工况1)

图4 stress-xx云图(工况1)

图5 stress-yy云图(工况1)工况2：即z2=-3，z1=7；模型顶部的自重应力为S1=0，底部的自重应力为S2=-1e5Pa。由式(3)计算得到梯度grad = (-1e5-0)/(-3-7)=1e4Pa/m；S(0) = -1e4*(-3)+(-1e5)=-7e4Pa。因此工况2的命令流如下：

;模型底面z坐标为-3model newzone create brick point 0(0,0,-3) point 1(10,0,-3) point 2(0,10,-3) point 3(0,0,7) size 10 10 10zone cmodel assign mohr-coulombzone property bulk 3e8 shear 1e8 cohesion 1e6 friction 15 zone initialize density 1000zone initialize stress-zz -7e4 grad 0 0 1e4zone initialize stress-xx -3.5e4 grad 0 0 5e3zone initialize stress-yy -3.5e4 grad 0 0 5e3

运行上述命令流后，模型的三个正应力云图如下：

图6 stress-zz云图(工况2)

图7 stress-xx云图(工况2)

图8 stress-yy(工况2)

======总结======

计算梯度时，首先要知道线性项的首位坐标以及对应的项值，通过直线方程求得斜率(grad)；令z=0即可求得S(0)。施加项可包括应力、速度、孔压等项目，当所施加项与x/y/xyz轴成线性时，计算方法相同。

-