目录

常规三轴

成样

预压

围压加载

二维直剪

成样

预压

围压加载

---

常规三轴

---

成样

---

cylinder keyword ...(3D ONLY)

Generate a cylinder in 3D. If the name keyword has not been specified, then s = cylinderWall. By default, the cylinder is closed and each side of the cylinder is a separate wall. In this case, the wall name is concatenated with the wall side name (e.g., top, bottom, or side) and the wall ID. For instance, if s = fred, the name of the top wall is fredTopXX, where XX is the wall ID. A unique ID, starting with id, is assigned to each wall. Use the one-wall keyword to specify that all cylinder sides belong to a single wall.

        cap bbottom <btop >

        Indicate the inclusion status of the end caps. If one boolean is specified, then this state         is given to both end caps. If two booleans are specified, then the first one applies to         the bottom cap and the second to the top cap. By default, both caps are created.

new
def par
    sampleRadius=0.08  
    height=sampleRadius*4
    
    rdmin=0.006
    rdmax=0.009
    
    poro=0.12
end
@par

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

wall generate cylinder base 0 0 [-height*0.5*1.2] axis 0 0 1 ...
        radius [sampleRadius] height [height*1.2] cap false false  
；生成空心圆筒，base指底部位置
wall generate plane position 0 0 [height*0.5] ddir 0 dip 0 ；倾角
wall generate plane position 0 0 [-height*0.5] ddir 0 dip 0
ball distribute porosity @poro radius [rdmin] [rdmax] range cylinder ...
        end1 0 0 [-height*0.5+rdmin] end2 0 0 [height*0.5-rdmin] radius [sampleRadius-rdmin]
；把半径缩小一点，防止颗粒生成在wall上
        
        
ball attribute density 2e3 damp 0.7
cmat default model linear method deform emod 100e6 kratio 1.5
cycle 2000 calm 50
ball delete range cylinder ...
        end1 0 0 [-height*0.5] end2 0 0 [height*0.5] radius [sampleRadius] not
；删除wall之外的颗粒
wall property fric 0.2   
ball property fric 0.5
solve
save sample

---

预压

---

        二维改三维

restore sample

[trr=1e4]
[tzz=1e4]
[do_rSevro=true]
[do_zServo=true]  ；给出径向和法向伺服的开关
[servo_factor=0.5]
[vpOnWall = wall.vertex.find(1)]
[wp=wall.find(1)]
def wallinit
    wpCe=wall.find(1)
    wpUp=wall.find(2)
    wpDown=wall.find(3)
end
@wallinit
def calChiCun
    pos_xy=vector(wall.vertex.pos.x(vpOnWall),wall.vertex.pos.y(vpOnWall),0)
    wlr=math.mag(pos_xy)  ；计算圆环的半径
    wlz=wall.pos.z(wpUp)-wall.pos.z(wpDown)
end

def calStress
    calChiCun
    sumForce=0
    loop foreach ct wall.contactmap(wpCe)    
        sumForce+=contact.force.normal(ct) ；得到侧向应力
    endloop
    wsrr=sumForce/(2*math.pi*wlr*wlz)
    wzss=(-1)*0.5*(wall.force.contact.z(wpdown)-wall.force.contact.z(wpup))/(math.pi*wlr*wlr)
end

def getg
    zongKNR=100e6*2.0
    zongKNZ=100e6*2.0
    loop foreach ct wall.contactmap(wpCe)
        zongKNR+=contact.prop(ct,"kn")
    endloop
    loop foreach ct wall.contactmap(wpup)
        zongKNZ+=contact.prop(ct,"kn")
    endloop
    loop foreach ct wall.contactmap(wpdown)
        zongKNZ+=contact.prop(ct,"kn")
    endloop
    gz=1.0*servo_factor*math.pi*wlr*wlr/(zongKNZ*global.timestep)  
    gr=1.0*servo_factor*2*math.pi*wlr*wlz/(zongKNR*global.timestep)  
end

[sevro_freq=100]
[time_record=global.step-1]
def sevro_wall
    calStress
    if global.step>time_record then
        getg
        time_record=global.step+sevro_freq
    endif
    if do_rSevro =true then 
        rvel=gr*(wsrr-trr)    
        loop foreach vt wall.vertexlist(wpCe)
            pos_xy=vector(wall.vertex.pos.x(vt),wall.vertex.pos.y(vt),0)    
            direct=math.unit(pos_xy)
            vel_vector=direct*rvel
            wall.vertex.vel(vt)=vel_vector
        endloop 
    endif
    
    if do_zServo=true then
        zvel=gz*math.abs(math.abs(wzss)-tzz)
        if math.abs(wzss)<tzz then
            wall.vel.z(wpup)=-zvel
            wall.vel.z(wpdown)=zvel
        else
            wall.vel.z(wpup)=zvel
            wall.vel.z(wpdown)=-zvel
        endif 
    endif
end

set fish callback -1.0 @sevro_wall

history id 1 @wsrr
history id 2 @wzss
cycle 1
solve

save yuya

---

围压加载

---

加载围压

restore yuya

[trr=1e5]
[tzz=1e5]

cycle 1
solve
save weiya

试样加载

restore weiya
[do_zServo=false]  ；关上轴向伺服
ball attribute displacement multiply 0
[strainRate=1e-2]
wall attribute zvel [strainRate*wlz] range id 3
wall attribute zvel [-strainRate*wlz] range id 2

[Iz0=wlz]
[Ir0=wlr]
def computer_strain
    wlz=wall.pos.z(wpup)-wall.pos.z(wpdown)
    wezz=(Iz0-wlz)/Iz0
    werr=(Ir0-wlr)/Ir0
    wevol=wezz+werr+werr ；监测应变
end
set fish callback -1.0 @computer_strain
history delete 
history id 1 @wsrr
history id 2 @wzss
history id 3 @wezz
history id 4 @werr
history id 5 @wevol

[stop_me=0]
def stop_me
    if wezz>0.2 then
        stop_me=1
    endif
end
solve fishhalt @stop_me

---

二维直剪

---

成样

---

加载板发生改变，对应代码也需要重新编写。

​​​​​​

 

new
def par
    width=0.9
    height=width*0.5
    
    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

[scale=1.5]  ；用来增加盒子的高度
wall generate id 1 box [-width*0.5*2.0] [-width*0.5]  ...
            0 [height*0.5*scale] onewall
            
wall generate id 2 box  [width*0.5] [width*0.5*2.0] ...
            0 [height*0.5*scale] onewall
            
wall generate id 3 box [-width*0.5*2.0] [-width*0.5]  ...
             [-height*0.5*scale] 0 onewall
            
wall generate id 4 box  [width*0.5] [width*0.5*2.0] ...
             [-height*0.5*scale] 0 onewall
             
wall generate id 5 box  [-width*0.5] [width*0.5] ...
             [height*0.5] [height*0.5*scale]  onewall
；生成上加载板
wall generate id 6 box  [-width*0.5] [width*0.5] ...
             [-height*0.5*scale] [-height*0.5]  onewall
；生成下加载板

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

[tyy=1e4]
def wall_ini
    wpup=wall.find(5)
    wpdown=wall.find(6)
    wpleftUp=wall.find(1)
    wpRightUp=wall.find(2)
    wpleftDown=wall.find(3)
    wpRightDown=wall.find(4)
end
；定义不同的wall
@wall_ini

def computer_stress
    wyss=0.5*(wall.force.contact.y(wpdown)-wall.force.contact.y(wpup))/width
end
；只用保留y向应力

[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
    
    
    gy=1.0*servo_factor*width/(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  
   
end

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

---

围压加载

---

restore yuya

[tyy=1e5]

cycle 1
solve
save weiya

restore weiya

set mech age 0
ball attribute displacement multiply 0
[strainRate=1e-3]
wall attribute xvel [strainRate*width] range id 1
wall attribute xvel [strainRate*width] range id 2
wall attribute xvel [strainRate*width] range id 5
；对上方的墙体施加x方向的速度

def jiance
    whilestepping
    weiyi=strainRate*width*mech.age
    wexx=weiyi/width
    wsxx=0
    wsxx+=wall.force.contact.x(wpup)
    wsxx+=wall.force.contact.x(wpleftUp)
    wsxx+=wall.force.contact.x(wpRightUp)
    wsxx=wsxx/(width-weiyi) ；求解应力和位移
end
history delete 
history id 1 @wsxx
history id 2 @weiyi
history id 3 @wexx

[stop_me=0]
def stop_me
    if weyy>0.2 then
        stop_me=1
    endif
end
solve fishhalt @stop_me