示例 2: 超级椭球三轴试验模拟

该示例给出了超级椭球的三轴试验模拟

该例子给出了超级椭球的三轴固结及压缩试验模拟。该模拟使用 CompressionEngine 引擎来完成, 同时用户可通过 PyRunner 引擎来自定义 CompressionEngine 引擎。固结和压缩引擎的基本设置如下:

  • 固结

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

        triax=CompressionEngine(
        goalx = 1.0e5, # 围压 100 kPa
        goaly = 1.0e5,
        goalz = 1.0e5,
        savedata_interval = 10000,
        echo_interval = 2000,
        continueFlag = False,
        max_vel = 100.0,
        gain_alpha = 0.5,
        f_threshold = 0.01, 
        fmin_threshold = 0.01, #偏应力与平均应力的比值
        unbf_tol = 0.01, #不平衡力比值
        file='record-con'
    )

    • goalx = goaly = goalz = 围压

    • savedata_interval: 保存模拟信息 (轴应变,体应变,平均应力和偏应力) 到文件的DEM时步周期，在固结引擎中该选项已被弃用。

    • echo_interval: 显示模拟信息 (时步, 不平衡力比值, 各个方向的应力) 到命令行的DEM时步周期。

    • continueFlag: 用来表示继续固结或压缩的保留关键字。

    • max_vel: 墙体最大移动速度，可在固结时用于限制墙体的移动。

    • gain_alpha $\alpha$: 用来控制墙体下一步的移动速度。 $$v_w = \alpha \frac{f_w - f_g}{\Delta t \sum k_n}$$ 其中， $v_w$ 表示墙体在下一时步的移动速度; $f_w$ 和 $f_g$ 分别表示当前和目标接触力; $\Delta t$ 表示时步; $\sum k_n$ 表示墙体和所有接触颗粒的接触刚度之和。

    • f_threshold $\alpha_f$: 归一化的接触力偏离比值, 即, $$\alpha_f = \frac{|f_w-f_g|}{f_g}$$

    • fmin_threshold $\beta_f$: 偏应力 $q$ 和平均主应力 $p$ 的比值, 即, $$\beta_f = \frac{q}{p},\quad p= \frac{\sigma_{ii}}{3},\quad q=\sqrt{\frac{3}{2}\sigma_{ij}^{'}\sigma_{ij}^{'}}$$ $$\sigma_{ij} = \frac{1}{V}\sum_{c\in V}f_i^cl_j^c,\quad \sigma_{ij}^{'}=\sigma_{ij}-p\delta_{ij}$$ 其中 $\sigma_{ij}$ 表示颗粒集合体的应力张量, $\sigma_{ij}^{'}$ 表示应力张量的偏应力部分; $V$ 表示试样体积; $f^c$ 和 $l^c$ 分别表示接触 $c$ 的接触力和接触支向量; $\delta_{ij}$ 表示克罗纳克尔。

    • unbf_tol $\gamma_f$: 表示不平衡力比值，由下式给出 　　$$\gamma_f = \frac{\sum\limits_{B_i\in V} |\sum\limits_{j\in B_i}f^c_j+f^b_j|}{2\sum\limits_{c\in V}|f^c|}$$ 需要说明的是只有当 $\alpha_f$, $\beta_f$ 和 $\gamma_f$ 都达到设定值，固结过程才会结束。

    • file: 模拟信息输出的文件相关信息

  • 压缩

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

        triax=CompressionEngine(
        z_servo = False,
        goalx = 1.0e5,
        goaly = 1.0e5,
        goalz = 0.01,
        ramp_interval = 10,
        ramp_chunks = 2,
        savedata_interval = 1000,
        echo_interval = 2000,
        max_vel = 10.0,
        gain_alpha = 0.5,
        file='sheardata',
        target_strain = 0.4,
        gain_alpha = 0.5
    )

    • z_servo: 布尔值，用来控制z方向是应变加载还是应力加载，在应变加载状态下 goalz 是对应的应变加载速率。

    • ramp_interval: 逐级加载到稳定应变速率的时步。该变量用来模拟实际试验中逐级稳定加载的过程，用户可以通过设定一个很小的值来忽略这个变量的影响。

    • ramp_chunks: 逐级加载到稳定应变速率的区间数。与 ramp_interval 类似，用户可以通过设定较小值来忽略它的影响。

    • target_strain: 加载目标应变 (对数应变)，在应变达到该值的时候剪切模拟停止。

      注: 如果用户需要进行准静态的三轴剪切模拟，剪切速率要低于一定的阈值。特别地，我们使用加载惯性系数来表征加载的准静态性， 即当惯性系数 $I$ 低于 $10^{-3}$时可视为准静态剪切，其中惯性系数由下式给出 $$I = \dot{\epsilon_1}\frac{\hat{d}}{\sqrt{\sigma_0/\rho}}$$ 其中 $\dot{\epsilon_1}=\mathrm{d}\epsilon_1/\mathrm{d}t$ 表示剪切应变加载速率; $\hat{d}$ 和 $\rho$ 分别表示颗粒的平均直径和平均密度， $\sigma_0$ 表示围压。

我们这里给出了一个三轴剪切加载的例子作为参考:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

##################################################
from sudodem import plot, _superquadrics_utils
from sudodem._superquadrics_utils import *

import math
import random as rand
rand.seed(11245642)
####################################################
epsilon = 1.4
isSphere = False
boxsize=[16.2e-3,16.2e-3,16.2e-3]
############################
#一些辅助函数
#生成单一粒径的超椭球颗粒，它们的初始位置和朝向是随机的
def gen_sample(width, depth, height, r_max = 1.0e-3):
   b_num_tot = 5000
   for j in range(b_num_tot):
      rb = r_max*0.5
      x = rand.uniform(rb, width-rb)
      y = rand.uniform(rb, depth - rb)
      z = rand.uniform(rb, height - rb)
      r = r_max*0.5
      b = NewSuperquadrics2(r*1.25,r,r,epsilon,epsilon,mat,True,isSphere)
      b.state.pos=(x, y, z)
      O.bodies.append(b)
#材料参数定义
mat = SuperquadricsMat(label="mat1",Kn=3e4,Ks=3e4,frictionAngle=math.atan(0.1),density=2650e4) #define Material with default values
wallmat1 = SuperquadricsMat(label="wallmat1",Kn=1e6,Ks=0.,frictionAngle=0.) #define Material with default values
wallmat2 = SuperquadricsMat(label="wallmat2",Kn=1e6,Ks=0.,frictionAngle=0.) #define Material with default values
O.materials.append(mat)
O.materials.append(wallmat1)
O.materials.append(wallmat2)
#############################################
#####生成真三轴的墙体
############################################
O.bodies.append(utils.wall(0,axis=0,sense=1, material = 'wallmat1'))#left x
O.bodies.append(utils.wall(boxsize[0],axis=0,sense=-1, material = 'wallmat1'))#right x
O.bodies.append(utils.wall(0,axis=1,sense=1, material = 'wallmat1'))#front y
O.bodies.append(utils.wall(boxsize[1],axis=1,sense=-1, material = 'wallmat1'))#back y
	
O.bodies.append(utils.wall(0,axis=2,sense=1, material = 'wallmat2'))#bottom z
O.bodies.append(utils.wall(boxsize[2],axis=2,sense=-1, material = 'wallmat2'))#top z
#####生成颗粒

gen_sample(boxsize[0], boxsize[1], boxsize[2])

####用来保存模拟数据的函数
def savedata():
        O.save(str(O.time)+'.xml.bz2')  
        
#########################################  
####三轴引擎
#注意: 压缩引擎默认前六个形状对象是墙体，所以在生成颗粒和墙体时，墙体总是要先于颗粒生成
triax=CompressionEngine(
   goalx = 1.0e5, # confining stress 100 kPa
   goaly = 1.0e5,
   goalz = 1.0e5,
	savedata_interval = 10000,
	echo_interval = 2000,
	continueFlag = False,
	max_vel = 10.0,
   gain_alpha = 0.5,
   f_threshold = 0.01, #1-f/goal
   fmin_threshold = 0.01, #the ratio of deviatoric stress to mean stress
   unbf_tol = 0.01, #unblanced force ratio
	file='record-con',
)

############################


      
newton=NewtonIntegrator(damping = 0.3,gravity=(0.,0.,0.),label="newton",isSuperquadrics=1)
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Superquadrics_Aabb(),Bo1_Wall_Aabb()],verletDist=0.1e-3),
   InteractionLoop(
      [Ig2_Wall_Superquadrics_SuperquadricsGeom(), Ig2_Superquadrics_Superquadrics_SuperquadricsGeom()], 
      [Ip2_SuperquadricsMat_SuperquadricsMat_SuperquadricsPhys()],
      [SuperquadricsLaw()]
   ),
   triax,
   newton,
   PyRunner(command='quiet_ball()',iterPeriod=2,iterLast = 50000,label='calm',dead = False),
   PyRunner(command='savedata()',realPeriod=7200.0,label='savedata',dead = False)#saving data every two hours
]


#删除在真三轴剪切盒之外的颗粒
def del_balls():
        num_del = 0
        #墙体位置
        left = O.bodies[0].state.pos[0]#x
        right = O.bodies[1].state.pos[0]
        front = O.bodies[2].state.pos[1]#y
        back = O.bodies[3].state.pos[1]
        bottom = O.bodies[4].state.pos[2]#z
        top = O.bodies[5].state.pos[2]
        for b in O.bodies:
                if isinstance(b.shape,Superquadrics):
                        (x,y,z) = b.state.pos
                        if x > right or x < left or y < front or y > back or z < bottom or z > top:
                                O.bodies.erase(b.id)
                                num_del += 1
        print str(num_del)+" particles are deleted!"           
                               


def quiet_ball():
    global calm_num
    
    newton.quiet_system_flag = True
    if calm_num > 2000:
        O.engines[-2].iterPeriod= 5
        calm_num += 5
    else:
        calm_num += 2
    if calm_num > 40000:
        if calm_num < 40010:
            print "calm procedure is over"
            del_balls()
        if calm_num > 50000:
            O.engines[-2].dead = True
            O.save("init_assembly.xml.bz2")
            #固结引擎开始
            triax.dead=False
            print "consolidation begins"
    
    
O.dt=1e-4


triax.dead=True
calm_num = 0

命令行将会有如下信息显示

calm procedure is over
0 particles are deleted!
consolidation begins
Iter 50000 Ubf 0.483161, Ss_x:9.67986e-05, Ss_y:4.15364e-05, Ss_z:6.1184e-05, e:0.992618
Iter 52000 Ubf 0.461741, Ss_x:1.2412, Ss_y:1.06812, Ss_z:1.10116, e:0.854445
Iter 54000 Ubf 0.263699, Ss_x:1.75334, Ss_y:1.98567, Ss_z:1.62181, e:0.757469
Iter 56000 Ubf 0.118329, Ss_x:3.06302, Ss_y:3.54818, Ss_z:3.22986, e:0.679106
Iter 58000 Ubf 0.0110747, Ss_x:22.0177, Ss_y:21.3303, Ss_z:21.3485, e:0.614824
Iter 60000 Ubf 0.000577066, Ss_x:91.4199, Ss_y:89.6039, Ss_z:89.1057, e:0.58984
Iter 62000 Ubf 0.000128293, Ss_x:100.012, Ss_y:99.7251, Ss_z:99.6123, e:0.587807
consolidation completed!

Packing of superellipsoids after consolidation.

固结完成后 (如图 Fig. [3.7] 所示，(#figex2consol){reference-type=“ref” reference=“figex2consol”}), 我们可以保存固结完成后的颗粒状态到文件中

1

O.save("final_con.xml.bz2")

然后我们开始对固结后的试样进行剪切，参见下面的脚本：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

##################################################
from sudodem import plot, _superquadrics_utils

from sudodem._superquadrics_utils import *

import math

saving_num = 0
####################################################
def savedata():
        global saving_num
        saving_num += 1
        O.save('shear'+str(saving_num)+'.xml.bz2')  
        
def quiet_system():
    for i in O.bodies:
        i.state.vel=(0.,0.,0.)
        i.state.angVel=(0.,0.,0.)
suffix=""
O.load("final_con"+suffix+".xml.bz2")
triax = O.engines[3]
O.dt = 5e-5
friction = 0.5
O.materials[0].frictionAngle=math.atan(friction)
if 1:
   for itr in O.interactions:
      id = min(itr.id1,itr.id2)
      if id > 5:
         itr.phys.tangensOfFrictionAngle = friction

#########################
#剪切开始
#########################
def shear():
    quiet_system()
    triax.z_servo = False
    triax.ramp_interval = 10
    triax.ramp_chunks = 2
    triax.file="sheardata"+suffix
    triax.goalz = 0.01
    triax.goalx = 1e5
    triax.goaly = 1e5
    triax.savedata_interval = 2500
    triax.echo_interval = 2000
    triax.target_strain = 0.4
    
shear()

O.engines[-1].iterPeriod = 2500
O.engines[-1].realPeriod = 0.0

试样的剪切加载总应变为40%，剪切终止时的试样状态如图所示 Fig. 3.8{reference-type=“ref” reference=“figex2shear”}.

Packing of superellipsoids after shear.

试样在剪切阶段的宏观力学响应可以用应力张量来表征，该张量由下式子得到: $$\label{eqDEMstress} \sigma_{ij} = \frac{1}{V} \sum_{c \in V} f_i^c l_j^c$$ 其中 $V$ 表示试样总体积； $f^c$ 表示接触 $c$ 的接触力矢量； $l^c$ 表示接触对应的支向量，也就是连接接触两颗粒质心的矢量。 平均主应力 $p$ 和偏应力 $q$ 被定义为: $$p = \frac{1}{3}\sigma_{ii},\quad q = \sqrt{\frac{3}{2} \sigma_{ij}^{'}\sigma_{ij}^{'}}$$ 其中 $\sigma_{ij}^{'}$ 表示应力张量的偏应力部分 $\sigma_{ij}$.

考虑到立方体试样的加载是由六个刚性墙来伺服控制的，轴向应变 $\epsilon_z$ 和 体应变 $\epsilon_v$ 可以由墙体的位置近似的给出, 即, $$\epsilon_z = \int_{H_0}^{H}\frac{\mathrm{d}h}{h}=-\ln\frac{H_0}{H},\quad \epsilon_v = \epsilon_x+\epsilon_y+\epsilon_z=\int_{V_0}^{V}\frac{\mathrm{d}v}{v}=-\ln\frac{V_0}{V}$$ 其中 $H$ 和 $V$ 分别表示剪切状态下的墙体高度和体积, $H_0$ 和 $V_0$ 表示刚性强初始状态的高度和体积。 体应变正值表示试样的剪胀状态。

剪切过程的数据保存格式如下：

AxialStrain VolStrain meanStress deviatorStress
0.000511000266 0.000490480823 106090.98 15931.56
0.003011001516 0.002450117535 131061.41 85718.91
0.005511002766 0.003867756447 148076.65 131980.95
0.008011004016 0.004825115794 159924.28 163706.90 
0.010511005266 0.005409561745 168606.27 186520.14
0.013011006516 0.005673777079 175131.53 203243.28
0.015511007766 0.005650846281 180293.19 216248.17
0.018011009016 0.005393122254 184581.23 227010.23

进而我们可以绘制出三轴剪切过程的应力-应变曲线以及应力路径曲线，如图所示 Figs. 3.9{reference-type=“ref” reference=“figex2shearstress”} and 3.10{reference-type=“ref” reference=“figex2shearpath”}.

Variation of deviatoric stress ratio $q/p$ and volumetric strain $\epsilon_v$ with axial strain during shearing.

剪切过程中的应力应变曲线以及体应变曲线。

Loading path during shearing.

剪切过程中的应力路径变化。

在剪切过程中，试样状态会不断地被保存在文件中，(例如， "shearxxx.xml.bz2") 用于随后的后处理，比如我们可以利用试样加载的元数据来获得试样的配位数和组构等信息。