PFC岩石双轴

目录

岩石双轴

加胶结

围压与加载

裂纹显示

---

岩石双轴

---

加胶结

---

 成样

new
def par
    width=0.4
    height=width*2
    
    rdmin=0.006
    rdmax=0.009
    
    poro=0.12
end
@par

domain extent [-width*2.0] [width*2.0] ...
        [-height*2.0] [height*2.0]
set random 10001

wall generate box [-width*0.5] [width*0.5] ...
            [-height*0.5] [height*0.5] expand 1.5
ball distribute porosity @poro radius [rdmin] [rdmax] box ...
    [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] 
ball attribute density 2e3 damp 0.7
cmat default model linear method deform emod 100e6 kratio 1.5
cycle 2000 calm 50
ball property fric 0.5
solve
save sample

预压

restore sample

[txx=1e4]
[tyy=1e4]
def wall_ini
    wpup=wall.find(3)
    wpdown=wall.find(1)
    wpleft=wall.find(4)
    wpright=wall.find(2)
end
@wall_ini

def computer_chicun
    wlx=wall.pos.x(wpright)-wall.pos.x(wpleft)
    wly=wall.pos.y(wpup)-wall.pos.y(wpdown)
end
def computer_stress
    computer_chicun
    wyss=0.5*(wall.force.contact.y(wpdown)-wall.force.contact.y(wpup))/wlx
    wxss=0.5*(wall.force.contact.x(wpleft)-wall.force.contact.x(wpright))/wly
end

[servo_factor=0.8]
def get_g

    zongKNY=100e6*2.0
    loop foreach ct wall.contactmap(wpup)
        zongKNY+=contact.prop(ct,"kn")
    endloop
    loop foreach ct wall.contactmap(wpdown)
        zongKNY+=contact.prop(ct,"kn")
    endloop
    
    zongKNX=100e6*2.0
    loop foreach ct wall.contactmap(wpleft)
        zongKNX+=contact.prop(ct,"kn")
    endloop
    loop foreach ct wall.contactmap(wpright)
        zongKNX+=contact.prop(ct,"kn")
    endloop
    
    gx=1.0*servo_factor*wly/(zongKNX*global.timestep) 
    gy=1.0*servo_factor*wlx/(zongKNY*global.timestep)  
end

[sevro_freq=100]
[time_record=global.step-1]
def sevro_walls
    computer_stress
    if global.step>time_record then
        get_g
        time_record=global.step+sevro_freq
    endif    
    yvel=gy*math.abs(math.abs(wyss)-tyy)
    if math.abs(wyss)<tyy then
        wall.vel.y(wpup)=-yvel
        wall.vel.y(wpdown)=yvel
    else
        wall.vel.y(wpup)=yvel
        wall.vel.y(wpdown)=-yvel
    endif  
    xvel=gx*math.abs(math.abs(wxss)-txx)
    if math.abs(wxss)<txx then
        wall.vel.x(wpleft)=xvel
        wall.vel.x(wpright)=-xvel
    else
        wall.vel.x(wpleft)=-xvel
        wall.vel.x(wpright)=xvel
    endif  
end

set fish callback -1.0 @sevro_walls 
history id 1 @wxss
history id 2 @wyss
cycle 10
solve
save yuya

加胶结

contact method bond gap [rdmin*0.2]

添加胶结，胶结间距一般取rdmin*0.2。间距越大，获得的胶结越多；间距过小可能会出现悬浮颗粒。

restore yuya

cmat add 1 model linearpbond method deform emod 1e8 kratio 1.5 ...
                            pb_deform emod 1e8 kratio 1.5 ...
                            property pb_coh 1e6 pb_ten 1e6 pb_fa 40 fric 0.5 ...
                            range contact type ball-ball
；add接触模型必须指定范围（range）
cmat apply
；改变当前接触模型
cycle 1
solve 
save temp1
contact method bond gap [rdmin*0.2]
；先平衡，平衡之后再施加胶结
；当胶结距离比较小的时候会出现一些悬浮颗粒，需要适当调大胶结距离

cycle 1
solve
save jiajiaojie

怎么看是否有胶结（监测）——重点是移除墙体伺服

不加胶结：删除墙体后接触消失

加胶结之后：删除墙体之后仍然有接触

restore jiajiaojie 
wall delete 

set fish callback -1.0 remove @sevro_walls 
；需要把墙体的伺服给移除

cycle 1
solve

---

围压与加载

---

加围压

        x方向和y方向分别施加0.5MPa围压。

restore jiajiaojie
[txx=0.5e6]
[tyy=0.5e6]

cycle 1
solve

save weiya

加载——调试显示应力应变曲线

restore weiya

ball attribute displacement multiply 0
def sevro_walls
    computer_stress
    if global.step>time_record then
        get_g
        time_record=global.step+sevro_freq
    endif  
    xvel=gx*math.abs(math.abs(wxss)-txx)
    if math.abs(wxss)<txx then
        wall.vel.x(wpleft)=xvel
        wall.vel.x(wpright)=-xvel
    else
        wall.vel.x(wpleft)=-xvel
        wall.vel.x(wpright)=xvel
    endif  
end
[strainRate=1e-2]
wall attribute yvel [strainRate*wly] range id 1
wall attribute yvel [-strainRate*wly] range id 3

[Iy0=wly]
[Ix0=wlx]
def computer_strain
    weyy=(Iy0-wly)/Iy0
    wexx=(Ix0-wlx)/Ix0
    wevol=weyy+wexx
end
set fish callback -1.0 @computer_strain
history delete 
history id 1 @wxss
history id 2 @wyss
history id 3 @weyy
history id 4 @wexx
history id 5 @wevol

call fracture.p2fis
@track_init
[stop_me=0]
def stop_me
    if weyy>0.2 then
        stop_me=1
    endif
end
solve fishhalt @stop_me

---

裂纹显示

---

手册_Example Applications_Rock Testing

添加 fracture.p2.fis，相当于二维的工具库，进行裂纹的显示

在每个裂纹地方，都有一个dfn。

查看：dfn_text val中选取name（张拉裂纹或者剪切裂纹）和age（裂纹龄期）。

fragment可以查看胶结破坏的颗粒，会存在明显裂纹。

; fname: fracture.p3fis
;
; Simple environment to track fragmentation in a BPM.
; Track LinearPBond model "bond_change" events and turn them into fractures.
; Use Fragment logic and Ball Result logic to record fragemnt ids 
;
;==============================================================================

define add_crack(entries)    ；entries只是一个定义的名称，作为一个array
    local contact    = entries(1)
    local mode       = entries(2)  ；破坏模式，需要再手册里找破坏模式对应的array slot
    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')
    ;ret为半个裂纹长度
    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
    ；vert1、vert2裂纹端点信息
    local arg = array.create(4)
    arg(1) = 'vertices'
    arg(2) = 2
    arg(3) = vert1
    arg(4) = vert2
    crack_num = crack_num + 1
    
    if mode = 1 then 
        ; failed in tension
        dfn_label = dfn_label + '_tension'
    else if mode = 2 then
        ; failed in shear
        dfn_label = dfn_label + '_shear'
    endif
        
    global dfn = dfn.find(dfn_label)
    if dfn = null then
        dfn = dfn.add(0,dfn_label)
    endif
    local fnew = dfn.addfracture(dfn,arg)
    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 > 50    ；每大于50个，更新一次裂纹位置
        if frag_time < mech.age
            frag_time = mech.age
            crack_accum = 0
            command 
                fragment compute
            endcommand
            ; 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
        fragment register ball-ball
    endcommand
  ; activate fishcalls
    command
        set fish callback bond_break @add_crack  
        ；回调第一个函数，只有在胶结破坏的时候才会回调@add_crack
    endcommand
    ; reset global variables
    global crack_accum = 0
    global crack_num = 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()

end

;==============================================================================
; eof: fracture.p3fis