FLAC3D Theory and Background • Constitutive Models

Pillar Stability with IMASS Model


To view this project in FLAC3D, use the menu command Help ► Examples…. Choose “Imass/ Piller” and select “pillar.prj” to load. The project’s main data file is shown at the end of this example.

In this example a single pillar is modeled with width-to-height ratio (w:h) of 1 in plane strain condition and uniaxially compressed in a quasi-static fashion at a constant strain rate until post-peak behavior is reached. The stress-strain response of the pillar and contour of emer_weak_sloss at approximately 4% axial strain are shown in Figure 1 and Figure 2, respectively. It can be seen that the shear bands are formed across the entire width of the pillar and the core is yielding. #modelimass3pillar-fig-minstress` and Figure 4 present contours of minimum principal (most compressible) stress and maximum shear strain, respectively. This is a good example to investigate the effect of rock mass brittleness on pillar strength and the overall behavior of the pillar (change in_weak_multecrit and see how it impacts the pillar response).

The pillar peak strength is in good agreement with the empirical limits. This however will change if the pillar w:h ratio is increased to values above 1.5-2 whereby pillar core will remain intact in the model and the estimated peak strength will be greater than those predicted by the empirical relationships. Experience suggests that squat pillars with w:h ratios greater than 1.5-2 can carry stresses much higher than suggested by the empirical pillar strength criteria. This is likely due to the fact that “failure” in these criteria is normally defined through visual observation, which is generally limited to the exterior of the pillar.


Figure 1: Stress-strain response of the pillar.


Figure 2: Contour of sloss.


Figure 3: Contour of minimum principal (most compressible) stress.


Figure 4: Contour of maximum shear strain.

Data Files


model new
model large-strain off
model configure imass
; geometry
fish define geometry
  ; input 
  global hp = 10. ; pillar height
  global wr = 10. ; room width
  global wh_ratio = 1.0 ;pillar width-height ratio
  global z_size = 0.5 ; zone size
  ; derived 
  local half_hp = hp/2.
  local wp = hp * wh_ratio
  local half_wp = wp/2.
  local half_wr = wr/2.
  local half_w = half_wr + half_wp
  local roof_height = half_hp*2.0
  global roof_top = 3.0*half_hp + roof_height
  local xel_side = int( half_w / z_size)
  local zel_side = int( roof_height / z_size)
  local top_zel_side = int(zel_side * 0.7)
  global area_pillar = wp * z_size
    zone create brick point 0 [-half_w] 0 0 ...
                      point 1 0 0 0 ... 
                      point 2 [-half_w] [z_size] 0 ...
                      point 3 [-half_w] 0 [roof_height] ...
                      ratio 1 1 1 size [xel_side] 1 [zel_side] group 'rock'
    zone group 'room' range position-x [-half_w] [-half_wp] position-z 0 [half_hp]
    zone create brick point 0 [-half_w] 0 [roof_height] ...
                      point 1 0 0 [roof_height] ...
                      point 2 [-half_w] [z_size] [roof_height] ...
                      point 3 [-half_w] 0 [roof_top] ...
                      ratio 1 1 1.05 size [xel_side] 1 [top_zel_side] group 'rock'
zone reflect origin 0 0 0 normal 1 0 0 
zone reflect origin 0 0 0 normal 0 0 -1 
zone face skin
; boudary condition
zone face apply velocity-x 0 range group 'West' or 'East'
zone face apply velocity-y 0 range group 'North' or 'South'
zone face apply velocity-z 0 range group 'bouthom' or 'top'
; assign cmodel and properties
zone delete range group 'room'
zone cmodel assign imass
zone initialize density 2500.0
zone property in_mod_youngintact 30e9
zone property in_stren_gsi 60.0
zone property in_stren_mi 15.0
zone property in_stren_ucsi 50e6
zone property in_weak_multecrit 1.0
[global gps_top = list(gp.list)(gp.isgroup(::gp.list,'Top'))]
fish define stress_top
  return list.sum(gp.force.unbal(::gps_top)->z)/area_pillar
[global gp_top1 = gp.near(0,0,roof_top)]
fish define strain_zz
  return strain_zz = -100.0*2.0*gp.disp.z(gp_top1)/roof_top
[global appliedVel = 1.e-6] ; loading velocity
zone face apply velocity-z [ appliedVel] range group 'Bottom'
zone face apply velocity-z [-appliedVel] range group 'Top'
zone mechanical damp combined
history interval 500
fish history stress_top
fish history strain_zz
model cycle 500000
model save 'pillarImass'