yade-users team mailing list archive

[Shiva Prashanth Kumar <question661683@xxxxxxxxxxxxxxxxxxxxx>]

New question #661683 on Yade:
https://answers.launchpad.net/yade/+question/661683

Hello, 

I am trying to simulate the direct shear test using PFC (Python scripting). When I plot the graph of vertical strain and shear strain, it is not starting from zero. Please go through the code and do the needful in this regard. Please find the python script for DST below. 

new
set random 10001
domain extent -0.2 0.2 
[ballFriction=0.0]
[wallFriction=0.0]
cmat default model linear method deform emod 1.0e8 kratio 2.5  

def setup
 height=0.053               
 height_=-0.053
 width=0.053                
 width_=-0.053
 new_height_1=height+0.02  
 new_height_2=-(height+0.02)
 new_width_1=width+0.02     
 new_width_2=-(width+0.02)  
 new_width_3=-(width+0.05) 
 new_width_4=(width+0.05)
 fc=0.5                  
 a_damp=0.7              
 p_density=1200                
end
@setup

def make_walls
 command
  wall create vertices @new_width_2, @height  @new_width_1, @height  id 1, name top_wall 
  wall create vertices @new_width_2, @height_  @new_width_1, @height_ id 2, name bot_wall
  wall create vertices @width_,  0,  @width_, @new_height_1 id 3 
  wall create vertices @width,  0,  @width, @new_height_1 id 4   
  wall create vertices @width_,  @new_height_2  @width_, 0 id 5  
  wall create vertices @width,   @new_height_2  @width,  0 id 6   
  wall create vertices @new_width_3, 0 @width_, 0 id 7  
  wall create vertices @width, 0 @new_width_4, 0  id 8
 end_command
end
@make_walls

wall property kn=1e8

def wll_wp
      global wpy1 = wall.find(2)  
      global wpy2 = wall.find(1)  
      global wpx1 = wall.find(3)  
      global wpx2 = wall.find(4)  
      global wpx11 = wall.find(5) 
      global wpx22 = wall.find(6)            
end
@wll_wp 


ball distribute porosity 0.175 radius 0.001 0.004 box -0.053 0.053 -0.053 0.053 
ball attribute density 1700 damp 0.7
ball property fric @ballFriction
wall property fric @wallFriction

cycle 10 calm 10
set timestep scale
solve arat 1e-5
set timestep auto
cycle 10
solve arat 1e-5

define compute_vessel_dimensions
  global et3_xlen = wall.pos.x(wall.find(4)) - wall.pos.x(wall.find(3))
  global et3_ylen = wall.pos.y(wall.find(1)) - wall.pos.y(wall.find(2))  
  x_wpx1= wall.pos.x(wpx1)
  x_wpx2= wall.pos.x(wpx2)
  y_wpy1= wall.pos.y(wpy1)
  y_wpy2= wall.pos.y(wpy2)    
end

def _et3_wstresses                  
         et3_wss             
          _xfob = 0.5 * (wall.force.contact(wpx1,1)+wall.force.contact(wpx11,1) - wall.force.contact(wpx2,1)-wall.force.contact(wpx22,1))
          _yfob = 0.5 * (wall.force.contact(wpy1,2) - wall.force.contact(wpy2,2))
        if _et3_wax # 0.0
         et3_wsxx = _xfob/_et3_wax
        else
         et3_wsxx = 0.0
        end_if
        if  _et3_way # 0.0
         et3_wsyy = _yfob/_et3_way
        else
         et3_wsyy = 0.0
        end_if  
end

define et3_wss
    compute_vessel_dimensions
    _et3_way = et3_xlen*1.0
    _et3_wax = et3_ylen*1.0    
end

define et3_servo_gain_poly          
  if et3_servo_xon = 1 then
     sum_knx = 0.0
    clist = wall.contactmap(wpx1)
    loop foreach local cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx11)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx2)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx22)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop    
    sum_knx = 0.5 * sum_knx 
    if sum_knx # 0.0 then
      et3_servo_gx = et3_servo_alpha * _et3_wax / (sum_knx * global.timestep)
    else
      et3_servo_gx = 0.0
    end_if
  end_if
  ;
  if et3_servo_yon = 1 then
    sum_kny = 0.0
    clist = wall.contactmap(wpy1)
    loop foreach cp clist
      sum_kny = sum_kny + contact.prop(cp,'kn')
    end_loop
    clist = wall.contactmap(wpy2)
    loop foreach cp clist
      sum_kny = sum_kny + contact.prop(cp,'kn')
    end_loop
    sum_kny = 0.5 * sum_kny  
    if sum_kny # 0.0 then
      et3_servo_gy = et3_servo_alpha * _et3_way / (sum_kny * global.timestep)
    else
      et3_servo_gy = 0.0
    end_if
  end_if
end

define et3_servo_poly          
  _et3_wstresses
  if et3_servo_xon = 1
    local _sgn = math.sgn( et3_wsxx - et3_wsxx_req )
    if et3_wsxx # 0.0 then
      wall.vel(wpx1,1) = et3_servo_gx * (et3_wsxx - et3_wsxx_req)
      wall.vel(wpx11,1) = et3_servo_gx * (et3_wsxx - et3_wsxx_req)      
    else
      wall.vel(wpx1,1) = _sgn * _vmax
      wall.vel(wpx11,1) = _sgn * _vmax
    end_if
    if math.abs(wall.vel(wpx1,1)) > _vmax then
      wall.vel(wpx1,1) = _sgn * _vmax
      wall.vel(wpx11,1) = _sgn * _vmax
    end_if
    wall.vel(wpx2,1) = -wall.vel(wpx1,1)
    wall.vel(wpx22,1) = -wall.vel(wpx1,1)
  end_if
  ;
  if et3_servo_yon = 1
    _sgn = math.sgn( et3_wsyy - et3_wsyy_req )
    if et3_wsyy # 0.0 then
      wall.vel(wpy1,2) = et3_servo_gy * (et3_wsyy - et3_wsyy_req)
    else
      wall.vel(wpy1,2) = _sgn * _vmax
    end_if
    if math.abs(wall.vel(wpy1,2)) > _vmax then
      wall.vel(wpy1,2) = _sgn * _vmax
    end_if
    wall.vel(wpy2,2) = -wall.vel(wpy1,2)
  end_if
  ;
end

define et3_servo_gain

  _et3_gain_cnt = _et3_gain_cnt + 1
  if _et3_gain_cnt >= et3_servo_gain_cyc then
    _et3_gain_cnt = 0
    et3_servo_gain_poly         
  end_if
end

def m_tol
 y_tol=(et3_wsyy - et3_wsyy_req)/et3_wsyy_req
 x_tol=(et3_wsxx - et3_wsxx_req)/et3_wsxx_req 
 mech_solve=mech.solve("aratio")
end

history id 104 @y_tol
history id 105 @x_tol
history id 107 @mech_solve

set fish callback  1.0 @m_tol

[et3_wsxx_req = -75.0e3]
[et3_wsyy_req = -75.0e3]      
[et3_servo_xon=1.0]
[et3_servo_yon=1.0]           

set fish callback  1.0 @et3_servo_poly      

history id 101 @et3_wsxx 
history id 102 @et3_wsyy  
history id 103 @mech_solve

set orientation on
calm

[tol =  5e-3]
[_et3_gain_cnt = 0]
[et3_servo_gain_cyc = 50]
[et3_servo_alpha = 0.5]
[et3_servo_vmax = 10.0]     
[_vmax=et3_servo_vmax]

define stop_me
  et3_servo_gain
  if math.abs((et3_wsyy - et3_wsyy_req)/et3_wsyy_req) > tol
    exit
  endif
  if math.abs((et3_wsxx - et3_wsxx_req)/et3_wsxx_req) > tol
    exit
  endif  
  if mech.solve("aratio") > 1e-5
    exit
  endif
  command
   sav sample_ini.p2sav
  end_command   
  stop_me = 1
end 

set grav 0 0

solve fishhalt @stop_me

res sample_ini.p2sav

def ini_shear       
 et3_xlen_0=et3_xlen
 et3_ylen_0=et3_xlen
end
@ini_shear 

def get_ss 
wall_x_dis=(wall.disp.x(wpy2)+wall.disp.x(wpx1)+wall.disp.x(wpx2))/3.0   
wall_x_str=wall_x_dis/et3_ylen_0  
wall_y_dis=(wall.disp.y(wpy2)+wall.disp.y(wpx1)+wall.disp.y(wpx2))/3.0   
wall_y_str=wall_y_dis/et3_ylen_0                                        
Fx1=wall.force.contact(wpx1,1)
Fx2=wall.force.contact(wpx2,1)
Fx=Fx1+Fx2
stress_x = -Fx/et3_xlen_0
vol_str = (et3_xlen*et3_ylen-et3_xlen_0*et3_ylen_0)/(et3_xlen_0*et3_ylen_0)         
end

define stop_shear
  gain_cnt = gain_cnt + 1
  if gain_cnt >= et3_servo_gain_cyc then
    et3_servo_gain_poly      
    gain_cnt = 0
  endif
 if math.abs(wall_x_str) >= target then
    stop_shear = 1
  endif
end

set fish callback  1.0 @get_ss

history id  1   @wall_x_dis    
history id  2   @wall_x_str    
history id  3   @stress_x      
history id  4   @vol_str       
history id  5   @wall_y_str


[et3_servo_xon=0]     

ball prop fric 0.5
[rate = 1.0e-3]     
[target = 2.0e-1]   
[stop_shear = 0]
[gain_cnt=0]

wall attribute xvelocity [rate] range id 1 
wall attribute xvelocity [rate] range id 3 
wall attribute xvelocity [rate] range id 4 
wall attribute xvelocity [rate] range id 7 

set echo off
solve fishhalt @stop_shear
set echo on


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.