目录

双轴压缩

伺服原理

参数化成样

应力计算并预压

加围压

重定义函数与加载

双轴压缩

伺服原理

        通过位移给墙体施加一个力。（墙体没有质量，只能根据重叠量来受力）

        使用墙给一个颗粒施加力

new domain extent -10 10 ball create position 0 0 radius 1.0 ball attribute density 2e3 damp 0.7 ball fix vel spin cmat default model linear property kn 1e7 wall create vertices -1 1 1 1 def cal_chongdie chongdieliang=1e4/1e7 end @cal_chongdie ；计算重叠量 wall attribute yvel -1.0 range id 1 ；给一个y方向的速度 set timestep fix 1e-5 solve time [chongdieliang]

        使用墙给多个颗粒施加力

讯享网new domain extent -10 10 ball create position -2 0 radius 0.5 ball create position -1 0 radius 0.5 ball create position -0 0 radius 0.5 ball create position 1 0 radius 0.5 ball create position 2 0 radius 0.5 ball attribute density 2e3 damp 0.7 ball fix vel spin cmat default model linear property kn 1e7 wall create vertices -3 0.5 3 0.5 def cal_chongdie chongdieliang=1e4/(1e7*5) end @cal_chongdie ；重叠量为力除以刚度，墙体的力为1e4N ；一般不给试样施加力，施加应力 wall attribute yvel -1.0 range id 1 set timestep fix 1e-5 solve time [chongdieliang] 

        对上面的加载进行优化，使得cycle 1 步就能得到目标值。

讯享网

new domain extent -10 10 ball create position -2 0 radius 0.5 ball create position -1 0 radius 0.5 ball create position -0 0 radius 0.5 ball create position 1 0 radius 0.5 ball create position 2 0 radius 0.5 ball attribute density 2e3 damp 0.7 ball fix vel spin cmat default model linear property kn 1e7 wall create vertices -3 0.5 3 0.5 def cal_chongdie chongdieliang=(1e4-0)/(1e7*5) sudu=chongdieliang/global.timestep command wall attribute yvel [-sudu] range id 1 endcommand end ；对墙底施加速度 set fish callback 1.0 @cal_chongdie ；在时步计算之后进行回调 cycle 1 

        球体串联情况

讯享网new domain extent -10 10 ball create position 0 0 radius 1.0 ball create position 0 2 radius 1.0 ball attribute density 2e3 damp 0.7 ball fix vel spin range id 1 cmat default model linear property kn 1e7 wall create vertices -1 3 1 3 ;k=k1*k2/(k1+k2) ;k=k1/2 ；串联时两个球体刚度计算 def cal_chongdie chongdieliang=1e4/(1e7*0.5) end @cal_chongdie wall attribute yvel -1.0 range id 1 set timestep fix 1e-5 solve time [chongdieliang] wall attribute yvel -0.0 range id 1 solve

        三角形分布球体伺服

        whilestepping，stop_me

new domain extent -10 10 ball create position -1 0 radius 1.0 ball create position 1 0 radius 1.0 ball create position 0 [math.sqrt(3)] radius 1.0 ；生成三个颗粒，根号3 ball attribute density 2e3 damp 0.7 ball fix vel spin range id 1 ball fix vel spin range id 2 cmat default model linear property kn 1e7 wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1] ;1、力小的时候，下压 ;2、力大了，上浮 [wp=wall.find(1)] ；定义墙体指针 def get_force whilestepping wall_force=wall.force.contact.y(wp) end ；whilestepping的效果和回调函数一样，即函数结果将会随着计算的更新而更新 def apply_vel whilestepping if wall_force<1e4 then wall.vel.y(wp)=-1.0 else wall.vel.y(wp)=1.0 endif end ；wall_force向上，为正值 ；力小的话则墙体下压，否则墙体上浮 [stop_me=0] def stop_me if math.abs(wall_force-1e4)<1 then stop_me=1 endif end ；stop_me函数判断截止条件，被solve调用，solve fishhalt @stop_me调用了停止函数 ；当值为0时继续运行，当为1时停止运行 set timestep fix 1e-5 ；时步不能太大，否则容易使得力来回变动 ;solve fishhalt @stop_me cycle 1 ；使用cycle 1 来打破平衡 solve ；两个要求，一施加力达到要求，二达到平衡 ;cycle 2.2w 

        伺服优化（速度+伺服系数）

讯享网new domain extent -10 10 ball create position -1 0 radius 1.0 ball create position 1 0 radius 1.0 ball create position 0 [math.sqrt(3)] radius 1.0 ball attribute density 2e3 damp 0.7 ball fix vel spin range id 1 ball fix vel spin range id 2 cmat default model linear property kn 1e7 wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1] ;1、力小的时候，下压 ;2、力大了，上浮 [wp=wall.find(1)] def get_force whilestepping wall_force=wall.force.contact.y(wp) end [g=5e-4] def apply_vel whilestepping vel=g*math.abs(wall_force-1e4) ；速度=伺服系数乘以力的差值 if wall_force<1e4 then wall.vel.y(wp)=-vel else wall.vel.y(wp)=vel endif end [stop_me=0] def stop_me if math.abs(wall_force-1e4)<1 then stop_me=1 endif end set timestep fix 1e-5 ;solve fishhalt @stop_me cycle 1 solve ;g=1e-3 cycle 5k8 ;g=1e-2 波动 不平衡 伺服系数g不能太大 ;g=5e-3 cycle 1w4 ;g=2e-3 cycle 2w2 ;g=5e-4 cycle 1w4 ；伺服系数需要试算 

        试算伺服系数

new domain extent -10 10 ball create position -1 0 radius 1.0 ball create position 1 0 radius 1.0 ball create position 0 [math.sqrt(3)] radius 1.0 ball attribute density 2e3 damp 0.7 ball fix vel spin range id 1 ball fix vel spin range id 2 cmat default model linear property kn 1e7 wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1] ;1、力小的时候，下压 ;2、力大了，上浮 [wp=wall.find(1)] def get_force whilestepping wall_force=wall.force.contact.y(wp) end ；重叠量=力除以刚度 ；速度=重叠量除以时步 [servo_factor=0.8] def get_g ct_wp=contact.find("ball-facet",1) gangdu=contact.prop(ct_wp,"kn") g=1.0*servo_factor/(gangdu*global.timestep) end set fish callback 1.0 @get_g def apply_vel whilestepping vel=g*math.abs(wall_force-1e4) if wall_force<1e4 then wall.vel.y(wp)=-vel else wall.vel.y(wp)=vel endif end [stop_me=0] def stop_me if math.abs(wall_force-1e4)<1 then stop_me=1 endif end set timestep fix 1e-5 ;solve fishhalt @stop_me cycle 10 solve 

参数化成样

        加载板和试样之间的摩擦力为0，而颗粒与颗粒之间为0.5。（两者摩擦力取最小值）

讯享网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 ；给定一个random的数值，使得每次运行得到的结果相同 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 ；method可以用来定义一些宏观的参数，但emod即弹性模量，但是弹性模量也会自动被计算为刚度 cycle 2000 calm 50 ；让cycle尽量大，使得颗粒散开 ball property fric 0.5 ；最好给出fric solve save sample 

应力计算并预压

        预压模拟的是土体的先前条件。伺服频率

        对于墙体和颗粒接触，模量乘以2对应kn的数值。

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 ；假设初始时有一定的约束，防止颗粒过于松散，kn为0 loop foreach ct wall.contactmap(wpup) ；遍历，ct为临时变量，它会把墙体接触依次定义 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) ；wly表示对应的面积，把应力转化为力 gy=1.0*servo_factor*wlx/(zongKNY*global.timestep) ；当试样特别松散的时候，刚度可能出现为0的情况 ；时步正比于(质量/刚度) 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 ；每100步计算一次get_g 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 ；记录x和y方向应力的变化 cycle 10 ；多运行几下，打破平衡 solve save yuya

加围压

        利用callback和fish函数。

讯享网restore yuya ；采用伺服原理 [txx=100e3] [tyy=100e3] ；设置为100kpa ；改变txx和tyy会直接代入yuya命令部分的函数进行计算 ；因为callback存在，所以运行过程中会不断更新数值 cycle 1 solve save weiya

重定义函数与加载

        对于预压部分的伺服函数进行重定义，删掉竖向伺服，添加竖向墙体速度。同时设定hist并添加加载截止条件。

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 ；获得g这个伺服系数 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 [stop_me=0] def stop_me if weyy>0.2 then stop_me=1 endif end solve fishhalt @stop_me ；人工定义截止条件，当应变大于0.2时就停止