PFC_常规三轴压缩试验代码（刚性墙）

前言

难点在于圆柱形墙的伺服控制，网上公开代码较少。

制样

; fname: make_sample.p3dat
;
;  Create unbonded sample
;=============================================================================
model new
plot creat 'triaxial_test'
; Set the domain extent
model domain extent -0.05 0.05 -0.05 0.05 -0.1 0.1 condition destroy

contact cmat default model linear method deformability emod 1.0e9 kratio 0.0
contact cmat default property dp_nratio 0.5

; create walls extending past the edges of the sample
wall generate id 10 plane dip 0 dip-direction 0   position 0 0 0.05
wall generate id 20 plane dip 0 dip-direction 0   position 0 0 -0.05
wall generate id 30 cylinder axis 0 0 1 base 0 0 -0.06 height 0.12 radius 0.025 cap false
;pause key
model random 10002
ball distribute porosity 0.2 radius 1.5e-3 2.5e-3 ...
                box -0.025 0.025 -0.025 0.025 -0.05 0.05
ball attribute density 2850 damp 0.7
ball delete range cylinder end-1 0 0 -0.0475 end-2 0 0 0.0475 rad 0.0225 not
; Calm the system
model cycle 1000 calm 10
; Solve the system to a target limit (here the average force ratio)
; Use density scaling to quickly reach equilibrium
model mechanical timestep scale
model solve ratio-average 1e-3
model mechanical timestep auto
model calm

model save 'triaxial_unbonded'
;=============================================================================
; eof: make_sample.p3dat

平行粘结模型

; fname: parallel_bonded.p3dat
;
;  Create parallel bonded sample
;==============================================================================
model restore 'triaxial_unbonded'

contact model linearpbond range contact type 'ball-ball'
contact method bond gap 2.0e-4

; set linear stiffness using methods
contact method deform emod 1.0e9 krat 1.0

; set stiffness of bonds using methods
contact method pb_deform emod 1.0e9 krat 1.0

; set bond strengths
contact property pb_ten 10.0e6 pb_coh 10.0e6 pb_fa 0.0

; set some damping at the contacts
contact property dp_nratio 0.5

; set ball-ball friction to non-zero value
contact property fric 0.577 range contact type 'ball-ball'

; Reset ball displacement 
ball attribute displacement multiply 0.0

; Set the linear force to 0.0 and force a reset of the linear contact forces. 
contact property lin_force 0.0 0.0 0.0 lin_mode 1
ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 
model cycle 1
model solve ratio-average 1e-5
model save 'tri_parallel_bonded'

;==============================================================================
; eof: parallel_bonded.p3dat

预压（伺服）

建议看PFC帮助文档，搜wall servo command，你会清楚get_gain()究竟在干嘛。其他的应该都很好理解。控制好servo_fac的大小，太大了墙的速度会很大（一般都是小于的1才对），报错结果包括1、结果除以0。2、from no pointer to … class啥的。

; fname: ucs.p3dat
;
;  Perform an unconfined compressive strength test on a bonded sample
;==============================================================================
model restore 'tri_parallel_bonded'
; to get the wall point and obtain the vertexs of cylindrical wall

define wallpoint
    wp_up = wall.find(10)
    wp_down = wall.find(20)
    wp_rr = wall.find(30)
    loop foreach vt wall.vertexlist(wp_rr)
        vert = vt
    endloop
end
@wallpoint
;pause key
define size_compute
    x_pos=wall.vertex.pos.x(vert)
    y_pos=wall.vertex.pos.y(vert)
    wlr=math.sqrt(x_pos^2+y_pos^2)
    wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down)
end

define stress_compute
    size_compute
    lateral_area = 2*math.pi*wlr*wlz
    bottom_area = 0.25*math.pi*wlr^2
;    ws1 =  wall.force.contact.z(wp_down)
;    ws2 = wall.force.contact.z(wp_up)
    wszz = (wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/bottom_area
    wsrr1 = 0
    loop foreach ft wall.facetlist(wp_rr)
        ft_direction = wall.facet.normal(ft)
        loop foreach ct wall.facet.contactmap(ft)
            ct_force = contact.force.global(ct)
            wsrr1 += -(math.dot(ct_force, ft_direction))/lateral_area
        end_loop
    endloop
end

def get_gain(factor)
    gz=0
    gr=0
    zonggangduZ=0
    zonggangduR=0
    loop foreach ct wall.contactmap(wp_up)
        zonggangduZ +=contact.prop(ct,"kn")
    endloop
    
    loop foreach ct wall.contactmap(wp_down)
        zonggangduZ +=contact.prop(ct,"kn")
    endloop
    
    loop foreach ct wall.contactmap(wp_rr)
        zonggangduR +=contact.prop(ct,"kn")
    endloop
    
    gz=factor/(zonggangduZ*global.timestep)
    gr=factor/(zonggangduR*global.timestep)    
end

[tzz=-5.0e6]
[trr=-5.0e6]
[sevro_fac=0.05]
[do_zservo=true]
[do_rservo=true]
[maxvel=1.0]
[sevro_freq=100]
[timestepNow=global.step-1]
def servo_walls
    stress_compute
    if timestepNow<global.step then
        get_gain(sevro_fac)
        timestepNow += sevro_freq
    endif
    if do_zservo=true then
        z_vel=gz*(wszz-tzz)*bottom_area
        z_vel=math.sgn(z_vel)*math.min(maxvel, math.abs(z_vel))
        wall.vel.z(wp_up)=-z_vel
        wall.vel.z(wp_down)=z_vel
    endif
    
    if do_rservo=true then
        r_vel_mag=-gr*(wsrr1-trr)*lateral_area
        r_vel_mag=math.sgn(r_vel_mag)*math.min(maxvel, math.abs(r_vel_mag))
        loop foreach vt wall.vertexlist(wp_rr)
            mag=math.sqrt(wall.vertex.pos.x(vt)^2+wall.vertex.pos.y(vt)^2)
            direction_normal_x=wall.vertex.pos.x(vt)/mag
            direction_normal_y=wall.vertex.pos.y(vt)/mag            
            r_vel=vector(direction_normal_x,direction_normal_y,0)*r_vel_mag            
            wall.vertex.vel(vt)=r_vel 
        endloop
    endif    
    
end

fish callback add @servo_walls -1.0

[stop_me=0]
[tol=0.005]
define stop_me
    if math.abs((wszz-tzz)/tzz)>tol
        exit
    endif
    if math.abs((wsrr1-trr)/trr)>tol
        exit
    endif
    stop_me = 1
end

; record histories
fish history name 10 @wszz
fish history name 20 @wsrr1
fish history name 30 @z_vel


model solve fishhalt @stop_me

model save 'preloading'

加载

model restore 'preloading'

[wszz0=wszz]
[wlr0=wlr]
[wlz0=wlz]
; compute stress and strain
define stress_strain
    zstrain=(wlz - wlz0)/wlz0
    rstrain=(wlr - wlr0)/wlr0
    vstrain=zstrain + 2*rstrain
    zstress=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/bottom_area
    deviatoric_stress = zstress - wszz0
    confined_stress = wsrr1
end

[do_zservo=false]

model calm
ball attribute displacement multiply 0.0


[vel0=5e-2]

wall attribute velocity-z [-vel0] range id 10
wall attribute velocity-z [vel0] range id 20

[stop_load=0]
define stop_load
    stress_strain
    if zstrain>=0.04 then
        stop_load=1
    endif
end

fish history name 100 @zstrain
fish history name 200 @rstrain
fish history name 300 @vstrain
fish history name 400 @deviatoric_stress
fish history name 500 @confined_stress

model solve fishhalt @stop_load
model save 'load'

存在问题

常规三轴加载时是柔性膜，刚性不太行。PFC帮助文档提供了PFC和FLAC耦合实现的常规三轴加载试验，感兴趣可以去学习学习。