为什么出PFC6.0？看完它与FlAC耦合的柔性三轴实验就明白

model new def chicun_par    sample_rad=1    sample_hight=sample_rad*4
    keli_rdmin=0.04    keli_rdmax=0.06end@chicun_par

model domain extent [-sample_rad*2] [sample_rad*2] [-sample_rad*2] [sample_rad*2]   [-sample_hight*1.5] [sample_hight*1.5]

[n=1.4]model random 10001wall generate cylinder base 0 0 [-sample_hight*0.5*n] axis 0 0 1 ...                    height [sample_hight*n] radius [sample_rad] resolution 0.3 cap false false

wall generate plane position 0 0 [sample_hight*0.5] dip 0 dip-direction  0wall generate plane position 0 0 [-sample_hight*0.5] dip 0 dip-direction  0

ball distribute group "shiyang" radius [keli_rdmin] [keli_rdmax] porosity 0.28 ...    range cylinder end-1 0 0 [sample_hight*0.5-keli_rdmin] ...    end-2 0 0 [-sample_hight*0.5+keli_rdmin] radius [sample_rad-keli_rdmin]

cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 property  fric 0.2cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5  ball attribute density 2.7e3 damp 0.7model cycle 2000 calm 50

model solve

model save "sample"

model restore "sample"

ball property "fric" 0.5def wp_wall    wp_up=wall.find(2)    wp_down=wall.find(3)    wp_rr=wall.find(1)        loop foreach vt wall.vertexlist(wp_rr)        vert_in_ce=vt             endloopend@wp_wall

[tzz=-1e4][trr=-1e4][sevro_fac=0.5]

[do_zservo=true][do_rservo=true]



[sevro_freq=100][timestepNow=global.step-1]def servo_walls    computer_wallStress    if timestepNow<global.step then        get_gain(sevro_fac)        timestepNow+=sevro_freq    endif    if do_zservo=true then        z_vel=gz*(wszz-tzz)        wall.vel.z(wp_up)=-z_vel        wall.vel.z(wp_down)=z_vel    endif
    if do_rservo=true then        r_vel_mag=(-1)*gr*(wsrr-tzz)        loop foreach vt wall.vertexlist(wp_rr)            mag=math.sqrt(wall.vertex.pos.x(vt)^2+wall.vertex.pos.y(vt)^2)            fang_normal_x=wall.vertex.pos.x(vt)/mag            fang_normal_y=wall.vertex.pos.y(vt)/mag                        r_vel=vector(fang_normal_x,fang_normal_y,0)*r_vel_mag                        wall.vertex.vel(vt)=r_vel         endloop    endif    
end



def computer_chicun    x_pos=wall.vertex.pos.x(vert_in_ce)    y_pos=wall.vertex.pos.y(vert_in_ce)    wlr=math.sqrt(x_pos^2+y_pos^2)    wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down)end

def computer_wallStress    computer_chicun    ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianji    ce_mianji=2*math.pi*wlr*wlz    wsrr=0    loop foreach ft wall.facetlist(wp_rr)        ft_fangxiang=wall.facet.normal(ft)        loop foreach ct wall.facet.contactmap(ft)            force_in_facet=contact.force.global(ct)            wsrr+=-(math.dot(force_in_facet,ft_fangxiang))/ce_mianji               endloop    endloop
end

def get_gain(fac)    gz=0    gr=0    zonggangZ=0    zonggangR=0    loop foreach ct wall.contactmap(wp_up)        zonggangZ+=contact.prop(ct,"kn")    endloop
    loop foreach ct wall.contactmap(wp_down)        zonggangZ+=contact.prop(ct,"kn")    endloop
    loop foreach ct wall.contactmap(wp_rr)        zonggangR+=contact.prop(ct,"kn")    endloop
    gz=fac*ding_yuanmianji/(zonggangZ*global.timestep)    gr=fac*ce_mianji/(zonggangR*global.timestep)    end





fish callback add @servo_walls -1.0

history id 1 @wszzhistory id 2 @wsrr

model cycle 1model solvemodel save "yuya"

model restore "yuya"

[tzz=-3e5][trr=-3e5]



model cycle 1 model solve

model save "weiya"

model restore "weiya"

fish callback remove @servo_walls -1.0[segments = 12][rad=wlr*1.01]

geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start ([rad*(-1)],0,[-wlz*0.6]) end (0,[rad*(-1)],[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start (0,[rad*(-1)],[-wlz*0.6]) end ([rad],0,[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start ([rad],0,[-wlz*0.6]) end (0,[rad],[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start (0,[rad],[-wlz*0.6]) end ([rad*(-1)],0,[-wlz*0.6]) ...   segments [segments];extrude the edges to make a cylindergeometry generate from-edges extrude (0,0,[wlz*1.2]) segments [segments*2]

structure shell import from-geometry 'Default' element-type dkt-cst structure shell property isotropic (1e6, 0.0) thick 0.25 density 930.0

structure node group 'top'  range position-z [wlz*0.48] [wlz*0.6]structure node group 'btm'  range position-z [-wlz*0.6] [-wlz*0.48]structure node group 'mid'  range position-z [-wlz*0.48] [wlz*0.48]structure shell group 'mid' range position-z [-wlz*0.48] [wlz*0.48]wall deletewall generate id 1 cylinder base 0 0 [wlz*0.5] axis 0 0 1 radius [wlr] height [keli_rdmax*10] one-wallwall generate id 2 cylinder base 0 0 [-wlz*0.5] axis 0 0 -1 radius [wlr] height [keli_rdmax*10] one-wall

structure damp localstructure shell apply [trr] range group 'mid'structure node fix velocity rotation range group "top"structure node fix velocity rotation range group "btm"wall servo force (0,0,[-tzz*math.pi*wlr*wlr]) activate true ...                 range id 2wall servo force (0,0,[tzz*math.pi*wlr*wlr]) activate true ...                 range id 1wall-structure create

measure create id 1 position 0 0 0  radius [wlr*0.4][mp=measure.find(1)]

def wp_wall    wp_up=wall.find(1)    wp_down=wall.find(2)end@wp_walldef jiance_measure       stressXX=measure.stress.xx(mp)    stressYY=measure.stress.yy(mp)    stressZZ=measure.stress.zz(mp)
    ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianjiendfish callback add @jiance_measure -1history deletehistory id 10 @stressXXhistory id 11 @stressYYhistory id 12 @stressZZhistory id 1 @wszz

model cycle 200model solve ratio-average 6e-5

model save "rouxing"

model restore "rouxing"

wall servo activate false range id 2wall servo activate false range id 1

history delete ball attribute displacement multiply 0[strainRate=5e-2]

wall attribute velocity-z [strainRate*wlz] range id 2wall attribute velocity-z [-strainRate*wlz] range id 1



structure node initialize velocity-z [strainRate*wlz] range group "btm"structure node initialize velocity-z [-strainRate*wlz] range group "top"

[Iz0=wall.pos.z(wp_up)-wall.pos.z(wp_down)-keli_rdmax*10]def jiance       stressXX=measure.stress.xx(mp)    stressYY=measure.stress.yy(mp)    stressZZ=measure.stress.zz(mp)    wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down)-keli_rdmax*10    wezz=(wlz-Iz0)/Iz0    ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianjiendfish callback add @jiance -1

history id 10 @stressXXhistory id 11 @stressYYhistory id 12 @stressZZhistory id 1 @wszzhistory id 2 @wezz

[stop_me=0]def stop_me    if wezz<-20e-2 then        stop_me=1    endifend

[baocunpinlv=1e-2][time_record=wezz+1][count=0]def savefile        if time_record-wezz >= baocunpinlv then        filename=string.build("jieguo_%1",count)        command           model save @filename        endcommand        time_record=wezz        count +=1    endif    endfish callback add @savefile -1.0

model solve fishhalt @stop_memodel save "result"