土工格栅拉拔试验2

define make_clump_bars(bar_id_start)

   global maxn_row_0=5

   global maxn_col_0=7

   global clump_idc=bar_id_start+1

   global bar_idc_start=clump_idc

   global ball_overlap_ratio_in_bar=1e-15

   global rc_joint_real=rc_joint*(1.0+ball_overlap_ratio_in_bar)

   global rc_rib_real=rc_rib*(1.0+ball_overlap_ratio_in_bar)

   maxn_row=maxn_row_0

   maxn_col=maxn_col_0

   global xc_joint=-90/unit

   global yc_joint=-80/unit

   global zc_joint=0/unit

   loop n_row(1,maxn_row)

       loop n_col(1,maxn_col)

           pebble_pos_vec1=vector(xc_joint,yc_joint,zc_joint)

           pebble_pos_vec2=vector(xc_joint-rc_joint,yc_joint,zc_joint)

           pebble_pos_vec3=vector(xc_joint+rc_joint,yc_joint,zc_joint)

           pebble_pos_vec4=vector(xc_joint,yc_joint-rc_joint,zc_joint)

           pebble_pos_vec5=vector(xc_joint,yc_joint+rc_joint,zc_joint)

           pebble_pos_vec6=vector(xc_joint-rc_joint,yc_joint-rc_joint,zc_joint)

           pebble_pos_vec7=vector(xc_joint+rc_joint,yc_joint-rc_joint,zc_joint)

           pebble_pos_vec8=vector(xc_joint-rc_joint,yc_joint+rc_joint,zc_joint)

           pebble_pos_vec9=vector(xc_joint+rc_joint,yc_joint+rc_joint,zc_joint)

           pebble_pos_vec10=vector(xc_joint-rc_joint-rc_rib,yc_joint,zc_joint)

           pebble_pos_vec11=vector(xc_joint-rc_joint-rc_rib,yc_joint-rc_rib,zc_joint)

           pebble_pos_vec12=vector(xc_joint-rc_joint-rc_rib,yc_joint+rc_rib,zc_joint)

           pebble_pos_vec13=vector(xc_joint+rc_joint+rc_rib,yc_joint,zc_joint)

           pebble_pos_vec14=vector(xc_joint+rc_joint+rc_rib,yc_joint-rc_rib,zc_joint)

           pebble_pos_vec15=vector(xc_joint+rc_joint+rc_rib,yc_joint+rc_rib,zc_joint)

           pebble_pos_vec16=vector(xc_joint,yc_joint-rc_joint-rc_rib,zc_joint)

           pebble_pos_vec17=vector(xc_joint-rc_rib,yc_joint-rc_joint-rc_rib,zc_joint)

           pebble_pos_vec18=vector(xc_joint+rc_rib,yc_joint-rc_joint-rc_rib,zc_joint)

           pebble_pos_vec19=vector(xc_joint,yc_joint+rc_joint+rc_rib,zc_joint)

           pebble_pos_vec20=vector(xc_joint-rc_rib,yc_joint+rc_joint+rc_rib,zc_joint)

           pebble_pos_vec21=vector(xc_joint+rc_rib,yc_joint+rc_joint+rc_rib,zc_joint)

           command

               ;;pause key

               ;ball delete range sphere center [pebble_pos_vec1] radius [rc_joint_real]

               ;ball delete range sphere center [pebble_pos_vec2] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec3] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec4] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec5] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec6] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec7] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec8] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec9] radius [rc_rib_real]

               ;

               ;ball delete range sphere center [pebble_pos_vec10] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec11] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec12] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec13] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec14] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec15] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec16] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec17] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec18] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec19] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec20] radius [rc_rib_real]

               ;ball delete range sphere center [pebble_pos_vec21] radius [rc_rib_real]

               clump create id=[clump_idc] pebbles 21   &

                           @rc_joint_real [pebble_pos_vec1] &

                           @rc_rib_real   [pebble_pos_vec2] &

                           @rc_rib_real   [pebble_pos_vec3] &

                           @rc_rib_real   [pebble_pos_vec4] &

                           @rc_rib_real   [pebble_pos_vec5] &

                           @rc_rib_real   [pebble_pos_vec6] &

                           @rc_rib_real   [pebble_pos_vec7] &

                           @rc_rib_real   [pebble_pos_vec8] &

                           @rc_rib_real   [pebble_pos_vec9] &

                           @rc_rib_real   [pebble_pos_vec10] &

                           @rc_rib_real   [pebble_pos_vec11] &

                           @rc_rib_real   [pebble_pos_vec12] &

                           @rc_rib_real   [pebble_pos_vec13] &

                           @rc_rib_real   [pebble_pos_vec14] &

                           @rc_rib_real   [pebble_pos_vec15] &

                           @rc_rib_real   [pebble_pos_vec16] &

                           @rc_rib_real   [pebble_pos_vec17] &

                           @rc_rib_real   [pebble_pos_vec18] &

                           @rc_rib_real   [pebble_pos_vec19] &

                           @rc_rib_real   [pebble_pos_vec20] &

                           @rc_rib_real   [pebble_pos_vec21] &

                           calculate 0.1

               ;ball delete range sphere center [xc_joint] [yc_joint] [zc_joint] radius [rc_joint]

               ;ball delete range sphere center [xc_joint-rc_joint] [yc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint+rc_joint] [yc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint] [yc_joint-rc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint] [yc_joint+rc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint-rc_joint] [yc_joint-rc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint+rc_joint] [yc_joint-rc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint-rc_joint] [yc_joint+rc_joint] [zc_joint] radius [rc_rib]

               ;ball delete range sphere center [xc_joint+rc_joint] [yc_joint+rc_joint] [zc_joint] radius [rc_rib]

               ;clump create id=[clump_idc] pebbles 9 &

               ;            @rc_joint [xc_joint] [yc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint-rc_joint] [yc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint+rc_joint] [yc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint] [yc_joint-rc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint] [yc_joint+rc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint-rc_joint] [yc_joint-rc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint+rc_joint] [yc_joint-rc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint-rc_joint] [yc_joint+rc_joint] [zc_joint] &

               ;            @rc_rib [xc_joint+rc_joint] [yc_joint+rc_joint] [zc_joint] &

               ;            calculate 0.1

               ;pause ;key

           end_command

           clump_idc=clump_idc+1

           xc_joint=xc_joint+40/unit

       end_loop

       yc_joint=yc_joint+40/unit

       xc_joint=-90/unit

   end_loop

   ;command

   ;    pause

   ;end_command

   maxn_row=maxn_row_0

   maxn_col=maxn_col_0-1

   global xc_joint=-90/unit

   global yc_joint=-80/unit

   global zc_joint=0/unit

   global maxnum=int((40/unit-d_joint-d_rib*2)/rc_rib)+1

   loop n_row(1,maxn_row)

       xc_joint=-90/unit

       loop n_col(1,maxn_col)

           xc_rib=xc_joint+rc_joint+d_rib  ;+rc_rib

           yc_rib=yc_joint

           loop num(1,maxnum)

               pebble_pos_vec1=vector(xc_rib,yc_rib,zc_rib)

               pebble_pos_vec2=vector(xc_rib,yc_rib-rc_rib,zc_rib)

               pebble_pos_vec3=vector(xc_rib,yc_rib+rc_rib,zc_rib)

               bp=ball.create(rc_rib_real,pebble_pos_vec1,clump_idc)

               clump_idc=clump_idc+1

               bp=ball.create(rc_rib_real,pebble_pos_vec2,clump_idc)

               clump_idc=clump_idc+1

               bp=ball.create(rc_rib_real,pebble_pos_vec3,clump_idc)

               ;command

               ;    ;pause key

               ;    ball delete range sphere center [pebble_pos_vec1] radius [rc_rib_real]

               ;    ball delete range sphere center [pebble_pos_vec2] radius [rc_rib_real]

               ;    ball delete range sphere center [pebble_pos_vec3] radius [rc_rib_real]

               ;    clump create id=[clump_idc] pebbles 3 &

               ;                @rc_rib_real [pebble_pos_vec1] &

               ;                @rc_rib_real [pebble_pos_vec2] &

               ;                @rc_rib_real [pebble_pos_vec3] &

               ;                calculate 0.1

               ;end_command

               clump_idc=clump_idc+1

               xc_rib=xc_rib+rc_rib

           end_loop

           ;command

           ;    pause key

           ;end_command

           xc_joint=xc_joint+40/unit

       end_loop

       yc_joint=yc_joint+40/unit

   end_loop

   maxn_row=maxn_row_0-1

   maxn_col=maxn_col_0

   global xc_joint=-90/unit

   global yc_joint=-80/unit

   global zc_joint=0/unit

   global maxnum=int((40/unit-d_joint-d_rib*2)/rc_rib)+1

   loop n_col(1,maxn_col)

       yc_joint=-80/unit

       loop n_row(1,maxn_row)

           yc_rib=yc_joint+rc_joint+d_rib    ;+rc_rib

           xc_rib=xc_joint

           loop num(1,maxnum)

               pebble_pos_vec1=vector(xc_rib,yc_rib,zc_rib)

               pebble_pos_vec2=vector(xc_rib-rc_rib,yc_rib,zc_rib)

               pebble_pos_vec3=vector(xc_rib+rc_rib,yc_rib,zc_rib)

               bp=ball.create(rc_rib_real,pebble_pos_vec1,clump_idc)

               clump_idc=clump_idc+1

               bp=ball.create(rc_rib_real,pebble_pos_vec2,clump_idc)

               clump_idc=clump_idc+1

               bp=ball.create(rc_rib_real,pebble_pos_vec3,clump_idc)

               ;command

               ;    ;pause key

               ;    ball delete range sphere center [pebble_pos_vec1] radius [rc_rib_real]

               ;    ball delete range sphere center [pebble_pos_vec2] radius [rc_rib_real]

               ;    ball delete range sphere center [pebble_pos_vec3] radius [rc_rib_real]

               ;    clump create id=[clump_idc] pebbles 3 &

               ;                @rc_rib_real [pebble_pos_vec1] &

               ;                @rc_rib_real [pebble_pos_vec2] &

               ;                @rc_rib_real [pebble_pos_vec3] &

               ;                calculate 0.1

               ;end_command

               clump_idc=clump_idc+1

               yc_rib=yc_rib+rc_rib

           end_loop

           ;command

           ;    pause key

           ;end_command

           yc_joint=yc_joint+40/unit

       end_loop

       xc_joint=xc_joint+40/unit

   end_loop

   bar_idc_end=clump_idc-1

   ;loop foreach clp clump.list

   ;    clump.density(clp)=bar_dens

   ;    clump.damp(clp)=loc_damp

   ;    clump.group(clp) = 'bar_joint'

   ;end_loop

end

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

def get_ss ; determine average stress and strain at walls

   new_height = wall.pos.z(wadd202) - wall.pos.z(wadd201)

   zdif  = new_height - height

   ;wszz  = 0.5*(wall.force.contact.z(wadd201) - wall.force.contact.z(wadd202)) / area

   wszz  = - wall.force.contact.z(wadd202) / area

   wezz  = 2.0 * zdif / (height + new_height)

   wevol = wezz

end

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

def get_gain ; determine servo gain parameters for axial and lateral motion

   global alpha = 0.5 ; relaxation factor

   global count = 0

   global avg_stiff = 0

   loop foreach contact wall.contactmap(wadd202)

       avg_stiff = avg_stiff + contact.prop(contact,"kn")

       count = count + 1

   endloop

   if get_bottom_avg_stiff = 1         ; switch stress servo on or off

       loop foreach contact wall.contactmap(wadd201)

           avg_stiff = avg_stiff + contact.prop(contact,"kn")

           count = count + 1

       endloop

   end_if

   ;loop foreach cp contact.list('ball-facet')

   ;    fp=contact.end2(cp)

   ;    if type.pointer(fp)='facet'

   ;        wp=wall.facet.wall(fp)  

   ;        if get_bottom_avg_stiff = 1         ; switch stress servo on or off

   ;            if wp = wadd201

   ;                count = count + 1

   ;                avg_stiff = avg_stiff + contact.prop(cp,'kn')

   ;            end_if

   ;        end_if

   ;        if wp = wadd202

   ;            count = count + 1

   ;            avg_stiff = avg_stiff + contact.prop(cp,'kn')

   ;        end_if

   ;    end_if

   ;end_loop

   ;

   ;

   ;loop foreach cp contact.list('pebble-facet')

   ;    fp=contact.end2(cp)

   ;    if type.pointer(fp)='facet'

   ;        wp=wall.facet.wall(fp)  

   ;        if get_bottom_avg_stiff = 1         ; switch stress servo on or off

   ;            if wp = wadd201

   ;                count = count + 1

   ;                avg_stiff = avg_stiff + contact.prop(cp,'kn')

   ;            end_if

   ;        end_if

   ;        if wp = wadd202

   ;            count = count + 1

   ;            avg_stiff = avg_stiff + contact.prop(cp,'kn')

   ;        end_if

   ;    end_if

   ;end_loop

   if count = 0 then

       exit

   end_if

   if get_bottom_avg_stiff = 1

       ncount = count / 2.0

   else

       ncount = count

   end_if

   avg_stiff = avg_stiff / count

   gz = alpha * length * width/ (avg_stiff * ncount * mech.timestep)

end

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

def servo

   if z_servo_close=1 then

       wall.vel.z(wadd201) = 0.0

       wall.vel.z(wadd201) = 0.0

   else

       get_ss                 ; compute stresses & strains

       udz = gz * (wszz - szzreq)

       if count = 0 then

           udz=math.abs(udz0)

       ;    wall.vel.y(wadd105) = udz

       ;    wall.vel.y(wadd106) = -udz

       ;else

       ;    wall.vel.y(wadd205) = 0.0

       ;    wall.vel.y(wadd206) = 0.0

       end_if

       wall.vel.z(wadd202) = -udz

       if bottom_servo_close = 1         ; switch stress servo on or off

           wall.vel.z(wadd201) = 0.0

       else

           wall.vel.z(wadd201) = udz

       end_if

   end_if

end

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

def get_wszz_tol

   cycle_cnt=cycle_cnt+1

   cycle_mod=int((float(cycle_cnt)/cycle_cnt_dist-cycle_cnt/cycle_cnt_dist)*cycle_cnt_dist)

   if cycle_mod=0

       get_gain

   end_if

   wszz_tol=math.abs((wszz - szzreq)/szzreq)

end

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

def iterate_limit

   if wszz_tol < sig_tol then

       iterate_limit=1

   end_if

end

;==============================================================================================

domain extent [-length*0.6] [length*1.1] &

             [-width*0.6]  [width*0.6]  &

             [-height*0.7] [height*0.7]

domain condition destroy

set orientation on

[rad=1.0]

clump template create name rock                     ...

           pebbles 4                               ...

           [rad] [-rad*0.5] 0 0                     ...

           [0.7*rad] [rad*0.6] 0 0                      ...

           [0.9*rad] 0 0 [rad*0.9]                      ...

           [0.6*rad] 0 [rad*0.2] [rad*1.6]                      ...

           pebcalculate 0.1

;geometry import rock.stl

;clump template create name rock                    ...

;            geometry rock                         ...

;            bubblepack ratio 0.3 distance 90       ...

;            surfcalculate

;pause ;key