三维砂土边坡的模拟

new def par    width=1.0    height=width*0.5      length=width*0.5    rdmax=0.009    rdmin=0.006end@pardomain extent -1 1wall generate box [-width*0.5] [width*0.5] [-length*0.5] [length*0.5] [-height*0.5] [height*0.5] expand 1.5ball distribute porosity 0.4 radius @rdmin @rdmax box [-width*0.5+rdmin] [width*0.5-rdmin] ...                                                        [-length*0.5+rdmin] [length*0.5-rdmin] ...                                                        [-height*0.5+rdmin] [height*0.5-rdmin]cmat default model linear method deformability emod 100e6 kratio 1.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solve save sample

restore sample ball property fric 0.5[txx=-12.5e3][tyy=-12.5e3][tzz=-12.5e3][sevro_factor=0.1][do_xSevro=true][do_ySevro=true][do_zSevro=true][sevro_freq=100][timestepNow=global.step-1]def sevro_walls    compute_stress     if timestepNow<global.step then        get_g(sevro_factor)        timestepNow+=sevro_freq    endif    if do_xSevro=true then        Xvel=gx*(wxss-txx)        wall.vel.x(wpRight)=-Xvel        wall.vel.x(wpLeft)=Xvel    endif        if do_zSevro=true then        Zvel=gz*(wzss-tzz)        wall.vel.z(wpUp)=-Zvel        wall.vel.z(wpDown)=Zvel    endif     if do_ySevro=true then        Yvel=gy*(wyss-tyy)        wall.vel.y(wpFront)=Yvel        wall.vel.y(wpBack)=-Yvel    endifenddef wp_ini    wpDown=wall.find(1)    wpUp=wall.find(2)    wpLeft=wall.find(3)    wpRight=wall.find(4)    wpFront=wall.find(5)    wpBack=wall.find(6)end@wp_inidef computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wlz=wall.pos.z(wpUp)-wall.pos.z(wpDown)    wly=-wall.pos.y(wpFront)+wall.pos.y(wpBack)enddef compute_stress    computer_chiCun    wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/(wly*wlz)    wzss=-(wall.force.contact.z(wpUp)-wall.force.contact.z(wpDown))*0.5/(wlx*wly)    wyss=(wall.force.contact.y(wpFront)-wall.force.contact.y(wpBack))*0.5/(wlx*wlz)enddef get_g(fac)    gx=0    gy=0    gz=0    zongKNX=1e6    zongKNY=1e6    zongKNZ=1e6    loop foreach ct wall.contactmap(wpLeft)        zongKNX+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpRight)        zongKNX+=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    loop foreach ct wall.contactmap(wpFront)        zongKNY+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpBack)        zongKNY+=contact.prop(ct,"kn")    endloop        gx=fac*wly*wlz/(zongKNX*global.timestep)    gy=fac*wlx*wlz/(zongKNY*global.timestep)    gz=fac*wlx*wly/(zongKNZ*global.timestep)end@get_g(@sevro_factor)@compute_stressset fish callback -1.0 @sevro_wallsdef addbond    command        contact method bond gap [MH_smalaradius*0.5]    endcommandendhistory id 1 @wxsshistory id 2 @wysshistory id 3 @wzsscycle 1solvesave yuya

restore yuyaset gravity [9.8*80.0/wlx]wall delete walls range id 2wall attribute vel 0set fish callback -1.0 remove @sevro_wallscycle 1 solvesave zizhong

restore zizhongdef sousuoFanwei    local xyzMinMax=array.create(3,2)        local x_min=1e100    local y_min=1e100    local z_min=1e100    local x_max=-1e100    local y_max=-1e100    local z_max=-1e100        loop foreach bp ball.list        if x_min>ball.pos.x(bp)-ball.radius(bp) then            x_min=ball.pos.x(bp)-ball.radius(bp)         endif        if x_max<ball.pos.x(bp)+ball.radius(bp) then            x_max=ball.pos.x(bp)+ball.radius(bp)         endif                if y_min>ball.pos.y(bp)-ball.radius(bp) then            y_min=ball.pos.y(bp)-ball.radius(bp)         endif        if y_max<ball.pos.y(bp)+ball.radius(bp) then            y_max=ball.pos.y(bp)+ball.radius(bp)         endif                if z_min>ball.pos.z(bp)-ball.radius(bp) then            z_min=ball.pos.z(bp)-ball.radius(bp)         endif        if z_max<ball.pos.z(bp)+ball.radius(bp) then            z_max=ball.pos.z(bp)+ball.radius(bp)         endif    endloop    xyzMinMax(1,1)=x_min    xyzMinMax(1,2)=x_max    xyzMinMax(2,1)=y_min    xyzMinMax(2,2)=y_max    xyzMinMax(3,1)=z_min    xyzMinMax(3,2)=z_max    sousuoFanwei=xyzMinMaxenddef ranse    local xyz_min_max=sousuoFanwei    local NWidth=10.0 ;横向方块数目    local NLength=5.0    local NHeight=5.0;纵向数目        local stepwidth=(xyz_min_max(1,2)-xyz_min_max(1,1))/NWidth    local steplength=(xyz_min_max(2,2)-xyz_min_max(2,1))/NLength    local stepheight=(xyz_min_max(3,2)-xyz_min_max(3,1))/NHeight    loop local n (1,NWidth)        local minwidth=xyz_min_max(1,1)+stepwidth*(n-1)                loop local m (1,NHeight)            local minlength=xyz_min_max(2,1)+steplength*(m-1)            loop k(1,NHeight)                 local minheight=xyz_min_max(3,1)+stepheight*(k-1)                if (n+m+k)/2-(n+m+k)/2.0=0 then                    command                        ball group gg1 slot 5 range x [minwidth] [minwidth+stepwidth] ...                                        y [minlength] [minlength+steplength] ...                                        z [minheight] [minheight+stepheight]                    endcommand                else                    command                        ball group gg2 slot 5 range x [minwidth] [minwidth+stepwidth] ...                                        y [minlength] [minlength+steplength] ...                                        z [minheight] [minheight+stepheight]                    endcommand                endif            endloop        endloop    endloopend@ransesave ranse

restore ranseball attribute displacement multiply 0set mech age 0[kengjiaoX=wlx*0.5*0.1][kengjiaoZ=-wlz*0.5*0.1][pojiao=70]ball delete range plane origin @kengjiaoX 0 @kengjiaoZ dip @pojiao dd 90 above ...        plane origin @kengjiaoX 0 @kengjiaoZ dip 0 dd 90 above        cycle 1solve time 10save qiepo