PFC5.0 2D单轴压缩模拟全代码分享

PFC5.0 2D Uniaxial Compression Simulation (单轴压缩模拟)

本文记录了使用 PFC 5.0 (Particle Flow Code) 进行二维单轴压缩试验模拟的全过程。模拟涵盖了试样生成、伺服固结、接触模型参数赋予、单轴加载以及裂纹监测等关键步骤。

1. 模拟流程概述

整个模拟过程主要分为以下几个阶段，分别对应不同的脚本文件：

  1. 试样生成 (main.p2dat): 生成圆柱形（2D中为矩形）试样，生成墙体和颗粒，并进行初步平衡。
  2. 伺服固结 (sifu.p2dat): 通过伺服机制控制墙体速度，使试样达到预定的各向同性应力状态（虽然是单轴，但预压过程有助于试样密实）。
  3. 参数赋予 (contact.p2dat): 赋予颗粒间平行粘结（Parallel Bond）模型及相应的微观参数。
  4. 单轴加载 (jaizai.p2dat): 移除侧向墙体，对上下墙体施加恒定速度进行加载，直至试样破坏。

此外，还包含了一些用于监测和后处理的 FISH 函数文件：

  • ss_wall.fis: 利用墙体监测应力-应变。
  • ss_gage.fis: 利用测量颗粒（Gage particles）监测应变。
  • fracture.p2fis: 监测和记录裂纹（张拉/剪切）的产生与发展。

2. 关键步骤详解

2.1 试样生成 (main.p2dat)

该脚本负责初始化模型。

  • 定义尺寸: 设置了圆柱体半径 cylinder_rad 和高度 height。
  • 生成墙体: 创建了上下加载板（id 1, 2）和侧向约束墙（id 3, 4）。
  • 生成颗粒: 使用 ball distribute 生成具有一定孔隙率和粒径分布的颗粒集合。
  • 初步平衡: 赋予初始摩擦和阻尼，通过 solve 求解至平衡状态。

2.2 伺服固结 (sifu.p2dat)

在赋予粘结模型之前，通常需要让试样在一定的围压下固结，或者至少保证颗粒间接触良好。

  • 伺服函数 (servo): 定义了一个伺服函数，根据当前墙体应力与目标应力（sxxreq, syyreq）的差值，自动调节墙体速度。
  • 迭代平衡: 通过 iterate 函数循环调用伺服机制，直到应力达到目标值（如 -1e5 Pa）。

2.3 接触模型与参数 (contact.p2dat)

这是模拟岩石或混凝土材料的关键。

  • 模型选择: 选择了 linearpbond (线性平行粘结模型)。
  • 参数赋予: 设置了模量 (emod, pb_modules)、刚度比 (kratio)、粘结强度 (pb_ten, pb_coh) 和摩擦角 (pb_fa)。这些微观参数决定了宏观材料的强度和变形特性。

2.4 单轴加载 (jaizai.p2dat)

模拟的核心加载过程。

  • 移除侧限: wall delete walls range id 3 4 删除了侧向墙体，实现单轴应力状态。
  • 施加荷载: 给上下墙体赋予恒定的相向速度 (yvel -0.05, yvel 0.05)。
  • 监测与结束: 调用 ss_wall.fis 和 fracture.p2fis 进行实时监测。定义了 loadhalt_wall 函数，当应力下降到峰值的一定比例（如 80%）时自动停止模拟。

2.5 裂纹监测 (fracture.p2fis)

利用 DFN (Discrete Fracture Network) 模块来记录粘结破坏。

  • add_crack 函数: 这是一个回调函数 (callback)，每当发生粘结破坏（bond_break）时自动触发。
  • 分类记录: 根据破坏模式（张拉 mode=1 或 剪切 mode=2）将裂纹分类记录，并统计裂纹数量和长度。

3. 所有代码文件

以下是本次模拟使用的所有源代码文件。

3.1 main.p2dat (试样生成)

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
new

title 'Testing Bonded Particle Model' 

[Rmin=1.5]

[Rmax=Rmin*1.66]

[cylinder_rad=250]    

[height=500]        



domain extent [-cylinder_rad*2] [cylinder_rad*2]  [-height*2]  [height*2] condition destroy



cmat default model linear method deform emod 1.0e9 kratio 0.0 

cmat default property dp_nratio 0.5  type ball-ball



wall create vertices [-cylinder_rad*0.7],[height/2]   [cylinder_rad*0.7],[height/2]  id 1    

wall create vertices [-cylinder_rad*0.7],[-height/2]  [cylinder_rad*0.7],[-height/2] id 2    

wall create vertices [-cylinder_rad/2],[-height*0.6]  [-cylinder_rad/2],[height*0.6] id 3     

wall create vertices [cylinder_rad/2],[-height*0.6]   [cylinder_rad/2],[height*0.6] id 4     



set random 10001

ball distribute porosity 0.1 radius [Rmin] [Rmax] box  [-cylinder_rad/2] [cylinder_rad/2]  [-height/2]  [height/2]

ball attribute density 2700 damp 0.7



cycle 1000 calm 10

set timestep scale

solve aratio 1e-5

set timestep auto

calm

save unbonded
call sifu.p2dat
call contact.p2dat

ret

3.2 sifu.p2dat (伺服固结)

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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

new

rest unbonded



def get_ss 

    wsyy  = 0.5*(wall.force.contact.y(wadd2) - wall.force.contact.y(wadd1)) / (cylinder_rad)       

    wsxx  = 0.5*(wall.force.contact.x(wadd3) - wall.force.contact.x(wadd4)) / (height)         

end

 

def get_gain 

    global alpha = 0.5

    global avg_stiff = 0.0

    

    loop foreach contact wall.contactmap(wadd1)

        avg_stiff = avg_stiff + contact.prop(contact,"kn")

    endloop

    

    loop foreach contact wall.contactmap(wadd2)

        avg_stiff = avg_stiff + contact.prop(contact,"kn")

    endloop

    

    if avg_stiff =0

        avg_stiff = 1e9

    endif

    gy = alpha *cylinder_rad/ (avg_stiff * mech.timestep)

    

    avg_stiff = 0.0

    loop foreach contact wall.contactmap(wadd3)

        avg_stiff = avg_stiff + contact.prop(contact,"kn")

    endloop

    

    loop foreach contact wall.contactmap(wadd4)

        avg_stiff = avg_stiff + contact.prop(contact,"kn")

    endloop

    

    if avg_stiff =0

        avg_stiff =1e9

    endif

    gx = alpha *height/ (avg_stiff * mech.timestep)

end



[max_v = 0.5 ] 

def servo

    get_ss                 ; compute stresses & strains

    if x_servo = 1         ; switch stress servo on or off

        udx = math.sgn(gx*(wsxx - sxxreq))*math.min(math.abs(gx*(wsxx - sxxreq)),max_v)

        wall.vel.x(wadd3) =  udx

        wall.vel.x(wadd4) = -udx

    end_if

    

    if y_servo = 1         ; switch stress servo on or off

        udy =  math.sgn(gy*(wsyy - syyreq))*math.min(math.abs(gy*(wsyy - syyreq)),max_v)

        wall.vel.y(wadd2) =  udy

        wall.vel.y(wadd1) = -udy

    end_if

end

set fish callback -1.0 @servo



def iterate

    loop while 1 # 0

        get_gain

        if math.abs((wsxx - sxxreq)/sxxreq) < sig_tol then

            if math.abs((wsyy - syyreq)/syyreq) < sig_tol then

                exit

            end_if

        end_if

        command

            cycle 100

        end_command

    end_loop

end



def wall_addr

    wadd1 = wall.find(1)

    wadd2 = wall.find(2)

    wadd3 = wall.find(3)

    wadd4 = wall.find(4)

end

@wall_addr



history id 11 @wsxx

history id 12 @wsyy



SET @sxxreq = -1e5 @syyreq = -1e5 @sig_tol = 0.05 @x_servo = 1 @y_servo = 1

@iterate  

solve aratio 1e-4

save sifu_comp

ret

3.3 contact.p2dat (接触参数)

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

new

rest sifu_comp



calm

set fish callback -1.0 remove  @servo

wall attribute velocity multi 0



cmat default model linear method deform emod 1.0e9 kratio 1.0  type ball-facet

cmat default model linear method deform emod 1.0e9 kratio 1.0  type ball-ball

[emod000=1e9]
[kratio=1.0]
[pb_modules=1e9]
[pb_kratio=1.0]
[ten_=13e6]
[coh_=50e6]
[pb_fa=30]
[fric=0.3]
contact model linearpbond range contact type ball-ball

contact method deform emod [emod000]  krat [kratio] range contact type ball-ball;emod 0.9e9 krat 2.5

contact method pb_deform emod [pb_modules] krat [kratio]  range contact type ball-ball;pb_deform emod 4.5e9 krat 2.5

contact property pb_ten [ten_] pb_coh [coh_] pb_fa [pb_fa] range contact type ball-ball

contact property dp_nratio 0.5  dp_sratio 0.0  range contact type ball-ball

contact property fric [fric] range contact type ball-ball

contact method bond gap 0.5e-4 range contact type ball-ball



ball attribute displacement multiply 0.0

contact property lin_force 0.0 0.0 lin_mode 1

ball attribute contactforce multiply 0.0 contactmoment multiply 0.0 

cycle 1

solve aratio 1e-5

save parallel_bonded


call jaizai.p2dat
;save ucs1


ret

3.4 jaizai.p2dat (加载)

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

new

res parallel_bonded



set echo off

call ss_wall.fis

call fracture.p2fis

@track_init

set echo on



set energy mech on 



wall delete walls range id 3 4 

wall attribute vel multiply 0.00 

wall attribute displacement multiply 0.00 

ball attribute displacement multiply 0.0



@setup_wall

wall attribute yvel -0.05 range id 1

wall attribute yvel  0.05 range id 2



[ball_left  =  ball.near(-cylinder_rad/2,0) ]

[ball_right =  ball.near( cylinder_rad/2,0) ]



def jiancei

   Elasticity_mod  = axial_stress_wall/(axial_strain_wall + 1e-6)

   Poisson_ratio   = math.abs((ball.disp.x(ball_right)-ball.disp.x(ball_left))/cylinder_rad/(axial_strain_wall + 1e-6))

end

set fish callback -1.0 @jiancei



history id 1 @axial_stress_wall

history id 2 @axial_strain_wall

history id 3 @Elasticity_mod

history id 4 @Poisson_ratio

history id 5 mech energy eboundary

history id 6 mech energy epbstrain 

history id 7 mech energy eslip

history id 8 mech energy ekinetic

history id 9 mech energy estrain


[nnnflag111=0
[nnnflag222=0]
define compute_elastic_modulus
axial_strain_wall1=math.abs(axial_strain_wall)
   if axial_strain_wall1>1.5e-3 then
     if nnnflag111 = 0 then
       strain1=axial_strain_wall1
       stress1=math.abs(peak_stress)
       nnnflag111=1
     endif
   endif
     if axial_strain_wall1>2e-3 then
       if nnnflag222 =0 then
         strain2=axial_strain_wall1
         stress2=math.abs(peak_stress)
         nnnflag222=1
         compute_elastic_modulus=(stress2-stress1)/(strain2-strain1)
           command
           set fish callback 9 remove @compute_elastic_modulus
           endcommand
        endif
    endif
end
set fish callback 9 @compute_elastic_modulus

cyc 1000

SET @peak_fraction = 0.8

solve fishhalt @loadhalt_wall
list @peak_stress
list @compute_elastic_modulus


save ucs1

ret

3.5 fracture.p2fis (裂纹监测)

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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
define add_crack(entries)
    local contact    = entries(1)
    local mode       = entries(2)
    local frac_pos   = contact.pos(contact)
    local norm       = contact.normal(contact)
    local dfn_label  = 'crack'
    local frac_size
    local bp1 = contact.end1(contact)
    local bp2 = contact.end2(contact)
    local ret = math.min(ball.radius(bp1),ball.radius(bp2));contact.method(contact,'pb_radius')
    frac_size = ret
    local inDir = vector(-comp.y(norm),comp.x(norm))
    local vert1 = frac_pos + inDir * frac_size
    local vert2 = frac_pos - inDir * frac_size
    local arg = array.create(4)
    arg(1) = 'vertices'
    arg(2) = 2
    arg(3) = vert1
    arg(4) = vert2
    crack_num = crack_num + 1
    jiaojiepohuaineng+= entries(4)
    
    
    pb_cohesion=contact.prop(contact,"pb_coh")
   if mode = 1 then
        ; failed in tension

         crack_tension_num+=1
        dfn_label = dfn_label + '_tension'

    else if mode = 2 then
        ; failed in shear
        crack_shear_num+=1
        failforce= entries(3)
        if failforce<pb_cohesion then
        
            dfn_label = dfn_label + '_shear_t'
            crack_shearT_num+=1
        else 
            dfn_label = dfn_label + '_shear_c'
            crack_shearC_num+=1
        endif
    endif 
    
    global dfn = dfn.find(dfn_label)
    if dfn = null then
        dfn = dfn.add(0,dfn_label)
    endif
    local fnew = dfn.addfracture(dfn,arg)
     if mode =1 then
        crack_tension_length+= dfn.fracture.len(fnew)
    endif
     if mode =2 then
        crack_shear_length+= dfn.fracture.len(fnew)
    endif
    crack_length+=dfn.fracture.len(fnew)
    dfn.fracture.prop(fnew,'age')  = mech.age
    dfn.fracture.extra(fnew,1) = bp1
    dfn.fracture.extra(fnew,2) = bp2
    crack_accum += 1
    if crack_accum > 25
        if frag_time < mech.age
            frag_time = mech.age
            crack_accum = 0
            command 
                fragment compute
            endcommand
            getfragInfo
            ; go through and update the fracture positions
            loop for (local i = 0, i < 2, i = i + 1)
                local name = 'crack_tension'
                if i = 1
                    name = 'crack_shear'
                endif
               
                dfn = dfn.find(name)
                if dfn # null
                    loop foreach local frac dfn.fracturelist(dfn)
                        local ball1 = dfn.fracture.extra(frac,1)
                        local ball2 = dfn.fracture.extra(frac,2)
                        if ball1 # null
                            if ball2 # null
                                local len = dfn.fracture.len(frac)/2.0
                                local pos = (ball.pos(ball1)+ball.pos(ball2))/2.0
                                if comp.x(pos)-len > xmin
                                    if comp.x(pos)+len < xmax
                                        if comp.y(pos)-len > ymin
                                            if comp.y(pos)+len < ymax
                                                dfn.fracture.pos(frac) = pos
                                            endif
                                        endif
                                    endif
                                endif
                            endif
                        endif
                    endloop
                endif
            endloop
        endif
    endif
end

define track_init
    command
        dfn delete
        ball result clear
        fragment clear
        fragment register ball-ball
        fragment compute
    endcommand
  ; activate fishcalls
    command
        set fish callback bond_break remove @add_crack
        set fish callback bond_break @add_crack
    endcommand
    ; reset global variables
    global crack_accum = 0
    global crack_num = 0
    global crack_tension_num=0
    global crack_shear_num=0
    global crack_shearT_num=0
    global crack_shearC_num=0
    global crack_tension_length=0
    global crack_shear_length=0
    global crack_length=0
    global fragNum=0
    global max_frage_vol=1
    global jiaojiepohuaineng=0
    global track_time0 = mech.age
    global frag_time = mech.age
    global xmin = domain.min.x()
    global ymin = domain.min.y()
    global xmax = domain.max.x()
    global ymax = domain.max.y()
    
   ; command
      ;  history id 50 @crack_num
      ;  history id 51 @crack_tension_num
       ; history id 52 @crack_shear_num
       ; history id 53 @crack_length
       ; history id 54 @crack_tension_length
       ; history id 55 @crack_shear_length
      ;  history id 56 @fragNum
      ;  history id 57 @max_frage_vol       
     ;   history id 58 @jiaojiepohuaineng 
      ;  history id 59 @crack_shearT_num
     ;   history id 60 @crack_shearC_num
  ;  endcommand

end

def getfragInfo
    local frag1 = fragment.find(1)   
    local catalognow=fragment.catalog.num(mech.age)
    fragNum=fragment.num(catalognow)
    max_frage_vol = fragment.vol(frag1,catalognow) / fragment.vol(frag1,1)      
    end

3.6 ss_wall.fis (墙体监测)

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
; fname: ss_wall.fis
;
;  Fish functions to record stress and strain in wall-loaded core samples
;========================================================================
def setup_wall
;
; Find the walls and get initial sample dimensions
;
  global wp_top = wall.find(1) ; assume wall 1 is the top wall
  global wp_bottom = wall.find(2) ; assume wall 2 is the bottom wall
  global vertical_direction = global.dim
  global sample_height = wall.pos(wp_top,vertical_direction) - wall.pos(wp_bottom,vertical_direction)
  
  ; assume x and y are approximately the same in 3D
  local xmin = 1.0e12
  local xmax = -1.0e12
  loop foreach bp ball.list
    local ball_xmin = ball.pos.x(bp) - ball.radius(bp)
    xmin = math.min(xmin,ball_xmin)
    local ball_xmax = ball.pos.x(bp) + ball.radius(bp)
    xmax = math.max(xmax,ball_xmax)
  end_loop
  
  local diameter_ = xmax - xmin
  if global.dim = 2
    global cross_sectional_area = diameter_
  else
    ; assume cylindrical sample in 3D
    cross_sectional_area = math.pi*0.25*diameter_*diameter_
  end_if
  
end
========================================================
def axial_stress_wall
;
; Compute axial stress (positive tension) using walls
;
; Assumes global variables sample_height and sample_width have been set
;
  
  local force1 = -wall.force.contact(wp_top,vertical_direction)
  local force2 = wall.force.contact(wp_bottom,vertical_direction)
  axial_stress_wall = 0.5*(force1+force2)/cross_sectional_area
end
========================================================
def axial_strain_wall
;
; Compute axial strain (positive tension) using walls
;
; Assumes global variable sample_width has been set
;
  axial_strain_wall = 2.0*wall.disp(wp_top,vertical_direction)/sample_height
end
==========================================================
def loadhalt_wall
;
; Function used to stop a test when axial stress decreases to some fraction of the peak
; Assumes axial_stress_wall is a function that returns axial stress
; 
; INPUT: peak_fraction - fraction of peak stress that dictates the stopping of the test
;                        (global float)
;
  loadhalt_wall = 0
  local abs_stress = math.abs(axial_stress_wall)
  global peak_stress = math.max(abs_stress,peak_stress)
  if abs_stress < peak_stress*peak_fraction
    loadhalt_wall = 1
  end_if
end
;==============================================================================
; eof: ss_wall.fis

3.7 ss_gage.fis (Gage颗粒监测)

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
; fname: ss_gage.fis
;
;  Fish functions to record strain using gage balls
========================================================
def setup_gage
; 
; Find gage particles to record axial strain when there are no walls
; and set initial sample height
;

  global vertical_direction = global.dim
  
  ; assume x and y are approximately the same in 3D
  local bottom_ = 1.0e12
  local top_ = -1.0e12
  loop foreach bp ball.list
    top_ = math.max(ball.pos(bp,vertical_direction),top_)
    bottom_ = math.min(ball.pos(bp,vertical_direction),bottom_)
  end_loop
  if global.dim = 2
    local top_loc = vector(0.0,top_)
    local bottom_loc = vector(0.0,bottom_)
  else
    top_loc = vector(0.0,0.0,top_)
    bottom_loc = vector(0.0,0.0,bottom_)
  end_if
  
  global gage_top = ball.near(top_loc)
  global gage_bottom = ball.near(bottom_loc)
  global sample_height = ball.pos(gage_top,vertical_direction) - ball.pos(gage_bottom,vertical_direction)
end
========================================================
def axial_strain_gage
;
; Compute axial strain using gage balls (positive tension)
;
; Assumes global variables gage_top, gage_bottom and sample_height have been set
;
  axial_strain_gage = (ball.disp(gage_top,vertical_direction) - ball.disp(gage_bottom,vertical_direction))/sample_height
end
==========================================================
def loadhalt_meas
;
; Function used to stop a test when axial stress decreases to some fraction of the peak
; Assumes measurement circle 1 is installed
; 
; INPUT: peak_fraction - fraction of peak stress that dictates the stopping of the test
;                        (global float)
;
  loadhalt_meas = 0
  local mp = measure.find(1)
  local abs_stress = math.abs(meas.stress(mp,global.dim,global.dim))
  global peak_stress = math.max(abs_stress,peak_stress)
  if abs_stress < peak_stress*peak_fraction
    loadhalt_meas = 1
  end_if
end
;==============================================================================
; eof: ss_gage.fis