这里成样时用的空隙率分别是0.5和0.7，其它参数一致。根据需求自己调松密，摩擦也是。

1. sample

new

def par

sampleRadius=0.025

height=sampleRadius*4

rdmin=0.002

rdmax=0.002

poro=0.50

end

@par

domain extent [-sampleRadius*2.0] [sampleRadius*2.0] ...

[-sampleRadius*2.0] [sampleRadius*2.0] ...

[-height*1.0] [height*1.0]

set random 10001

wall generate cylinder base 0 0 [-height*0.5*1.4] axis 0 0 1 ...

radius [sampleRadius] height [height*1.4] cap false false

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]

ball attribute density 2600 damp 0.7

cmat default model linear method deform emod 1e8 kratio 1.0

cycle 2000 calm 50

ball delete range cylinder ...

end1 0 0 [-height*0.5] end2 0 0 [height*0.5] radius [sampleRadius] not

solve aratio 5e-3

save sample

2. yuya

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 get_sample_ball_volume

global sample_ball_volume = 0

loop foreach local bp ball.list

sample_ball_volume = (math.pi*4*ball.radius(bp)^3)/3 + sample_ball_volume

endloop

end

set fish callback -1.0 @get_sample_ball_volume

;------------------------------------------------------------------------------------------------------------------

def get_voidratio

calChiCun

sample_volume = math.pi*wlr^2*wlz

sample_void_volume = sample_volume - sample_ball_volume

porosity = sample_void_volume/sample_volume

;voidratio = porosity/(1-porosity)

end

set fish callback -1.0 @get_voidratio

;------------------------------------------------------------------------------------------------------------------

def getg

voidratio = porosity/(1-porosity)

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 aratio 5e-3

save yuya

3. jiechu

restore yuya

cmat default model linear method deform emod 1e8 kratio 1.0 prop fric 0.3 type ball-ball

cmat default model linear method deform emod 1e7 kratio 0 prop fric 0.0 type ball-facet

cmat apply

cycle 1

solve aratio 5e-3

save jiechu

4. weiya

restore jiechu

[trr=1e5]

[tzz=1e5]

cycle 1

solve

save weiya

5. jiazai

restore weiya

[do_zServo=false]

[do_rSevro=false]

[wallVel1= 1e-2]

set fish callback -1.0 remove @sevro_wall

ball attribute displacement multiply 0

wall attribute velocity multiply 0

set mech age 0

wall attribute zvel [wallVel1] range id 3

wall attribute zvel [-wallVel1] range id 2

;------------------------------------------------------------------------------------------------------------------

def changeWallVel

WallVel=wlr/wlz*wallVel1

loop foreach local ve wall.vertexlist(wpCe)

rr=math.unit(vector(wall.vertex.pos.x(ve),wall.vertex.pos.y(ve)))

wall.vertex.vel(ve)=WallVel*rr

endloop

end

set fish callback -1.0 @changeWallVel

;------------------------------------------------------------------------------------------------------------------

[Iz0=wlz]

[Ir0=wlr]

[wsrr0=wsrr]

def computer_strain

time = mech.age

calChiCun

calStress

wezz=(Iz0-wlz)/Iz0

werr=(Ir0-wlr)/Ir0

wevol=wezz+werr+werr

pwp=wsrr0-wsrr

q = wzss -wsrr

p_e = (wzss+2.0*wsrr)/3.0

p_t = p_e+pwp

end

@computer_strain

set fish callback -1.0 @computer_strain

;------------------------------------------------------------------------------------------------------------------

history delete

history id 1 @wsrr ;the stress of radius direction

history id 2 @wzss ;the stress of z direction

history id 3 @wezz ;the strain of z direction

history id 4 @werr ;the strain of radius direction

history id 5 @wevol ;the strain of volume

history id 7 @pwp ;the pressure of water

history id 8 @time

history id 9 @q

history id 10 @p_e

history id 11 @p_t

history id 12 @voidratio

def loadhalt_wall

loadhalt_wall = 0

local abs_ = math.abs(wezz)

if abs_ > 0.10

loadhalt_wall = 1

end_if

end

solve fishhalt @loadhalt_wall

save jieguo

Leo

简介

擅长领域