new def par    rdmin=6e-3    rdmax=9e-3        width=4    height=width/3.0end@par

domain extent [-width] [width]

wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5

ball distribute porosity 0.5 radius [rdmin] [rdmax] box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] 

cmat default model linear

cmat add 1  model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 range contact type ball-ball

cmat add 2 model linear method deformability emod 100e6 kratio 1.5  range contact type ball-facet

ball attribute density 2.7e3 damp 0.7cycle 1000 calm 50cycle 2000set gravity 9.8solvesave sample

restore sample

[diding=height*0.5*0.22]ball delete range y [diding] [height*0.5]

cycle 1solvesave diji

restore diji[nianlun_rad=height*0.1][force=3.8e3*9.8*math.pi*nianlun_rad^2]wall generate id 10 name nianlun circle position [-width*0.5+nianlun_rad] [diding+nianlun_rad*1.01] radius [nianlun_rad] wall delete range id 3wall servo activate on gainupdate 100 vmax [height] yforce [-force] range id 10ball property fric 0.5wall property fric 0.5 range id 5 68cmat modify 2 model linear method deformability emod 100e6 kratio 1.5  range contact type ball-facet

cycle 1000solve aratio 1e-5save create_nianlun

restore create_nianlunset mech age 0ball attribute displacement multiply 0

[yisu=0.2][spin=yisu/nianlun_rad][cishu=5][mea_rad=width*0.05]

wall attribute xvel [yisu] spin [-spin] range id 10

[wp_nianlun=wall.find(10)]def updata_center    whilestepping    wall.rotation.center(wp_nianlun)=wall.pos(wp_nianlun)end[onetime=(width-2*nianlun_rad)/yisu][timerecord=0][file_count=1]

measure deletemeasure create id 1 position [-width*0.5*0.2] [-height*0.5*0.3] radius [mea_rad]measure create id 2 position [-width*0.5*0.1] [-height*0.5*0.3] radius [mea_rad]measure create id 3 position [width*0.5*0.1] [-height*0.5*0.3]  radius [mea_rad]measure create id 4 position [width*0.5*0.2] [-height*0.5*0.3]  radius [mea_rad] def mp_get    mp1 = measure.find(1)    mp2 = measure.find(2)    mp3 = measure.find(3)     mp4 = measure.find(4) end@mp_getdef jiance    whilestepping        measure_poros1=measure.porosity(mp1)    measure_poros2=measure.porosity(mp2)    measure_poros3=measure.porosity(mp3)    measure_poros4=measure.porosity(mp4)        measure_av=(measure_poros1+measure_poros2+measure_poros3+measure_poros4)*0.25    forceLunY = wall.force.contact.y(wp_nianlun)    forceLunX = wall.force.contact.x(wp_nianlun)    time_now = mech.age    weiyiX = yisu * mech.age    nianlun_pos = wall.pos.x(wp_nianlun)end@jiancedef circle_nianya    whilestepping    if mech.age-timerecord>onetime then        wall.vel.x(wp_nianlun)=-wall.vel.x(wp_nianlun)        wall.spin(wp_nianlun)=-wall.spin(wp_nianlun)        timerecord=mech.age        filename=string.build("jieguo%1",file_count)        command            save @filename        endcommand        file_count+=1    endifendhistory delete history id 1 @time_nowhistory id 2 @weiyiXhistory id 3 @forceLunYhistory id 4 @forceLunXhistory id 5 @measure_poros1history id 6 @measure_poros2history id 7 @measure_poros3history id 8 @measure_poros4history id 9 @measure_avsolve time [cishu*onetime]