Tailing Dam Static Liquefaction (FLAC3D)

Problem Statement

Note

The project file for this example is available to be viewed/run in FLAC3D.[1] The project’s main data file is shown at the end of this example.

This example demonstrates the analysis of static liquefaction potential due to the construction of an additional layer (simplified by applying a surcharge). Key analysis procedures illustrated in this example include:

• Automatic grid generation by controlling zone size through the import of a soil profile geometry file

• Staged construction simulation

• Phreatic surface (water table) calculation

• Static liquefaction analysis

The model includes three basic materials: the starter dam, dikes, and tailings. The profile is shown in Figure 1.

Figure 1: Model geometry of the tailing dam.

Modeling Procedure

Grid Generation

The model geometry is based on a DXF file named “tailingsdam-geometry.dxf”. The model is created interactively using the Sketch pane. The State Record pane was used to convert the result to a data file called “sketch.dat”. The contents of the Sketch pane can be be viewed after restoring any of the resulting data files.

The specific steps are as follows:

• After creating the project file, go to the “File” menu, select “New”, then “Sketch Set”, and input a set name “tailings”. A sketch pane called “tailings” will be created.

• In the “tailings” Sketch pane, click the “Import background image” button, then select the file “tailingsdam-geometry.dxf”. The geometry with nodes, edges, and polygons will be displayed in the Sketch pane.

Once the background image is loaded into the Sketch pane, automatically create edges from this image:

• Click the “Automatically create edges from background geometry” button.

• Click the “Mesh Tools” button, select “Auto-size”, and set the zoning using “zone length” with an input value of “1.0”.

• Click the “Mesh Tools” button again and select “Mesh All Polygons”.

At this point, the grid is automatically created with an approximate zone size of 1.0, which will be replaced by a FISH variable later. Next, assign group names to the sketch blocks and edges:

• Click the “Selection tool” button, click the sketch block of the starter dam, and input the group name “StarterDam” in the “Object Properties”.

• Similarly, input group names “Dike1” to “Dike7” for the sketch blocks of the dikes; input group names “Tailings01” to “Tailings10” for the sketch blocks of the tailings.

• Click the sketch edge of the uppermost tailing and input the group name “Load”.

After assigning group names to some sketch blocks and edges, proceed with the extrusion steps:

• Click “Bottom” view in the Sketch pane.

• Click the right point of the middle line. In the “Object Properties”, modify the Position’s W-coordinate to 1.

• Zoom in on the sketch, click or select the middle line. In the “Object Properties”, modify the “Zone number” to 1 and “Zone length” to 1.

• Click the “Create zones” button, the grid as shown in Figure 2 is created.

Figure 2: Soil Profile of the tailing dam.

All above operations in the Sketch pane are recorded in the State Record pane. It is good practice to save these records into a data file so that later modifications can be made directly to the data file:

• Click the “State Record” tab. Right-click anywhere in the “State Record” window.

• Select “Save to file as datafile” and input a data file name, e.g., “sketch.dat”.

To ensure the sketch operations can be duplicated, follow these steps:

• Open the newly created data file and add the line “model new” at the very beginning.

• Run the data file; the sketch operations should be replicated.

• Add a line at the end to save the model, e.g., model save ‘sketch’.

If we need to adjust the zone size, we can assign a global FISH variable in the main datafile and modify the corresponding values into this FISH variable:

• In the main data file, master.dat, input the first three lines as follows:

model new
[global zone_size = 1.0]
program call "sketch.dat"

• Comment out the first line of model new in sketch.dat in case it has not been commented out. Otherwise, the FISH variable in the main data file cannot be transferred to sketch.dat.

• Modify the zone size in the corresponding lines to use the assigned FISH variable, like

; model new
...
sketch set automatic-zone direction construction edge [zone_size]
...
sketch segment-node id 2 move-to (75,30.25,[zone_size]) tolerance-merge 1
...
model save 'sketch'

By following these steps, we can easily refine or coarsen the grid by simply adjusting the zone size variable in the main data file and re-running the script. This approach allows for quick and flexible modifications to the model’s resolution.

Staged Construction

This simulation is scripted in the datafile construction.dat.

After restoring the previous save file of grid generation, the command

zone face skin

is used to automatically assign face groups to the skins of the model. The face groups can be checked by plotting “zone->Face”. These face groups will be used later.

The model then sets some global configurations, including “fluid-flow”, “large-strain” and “gravity” by the commands

model configure fluid-flow
model large-strain off
model gravity 9.81

The boundary conditions are simply the standard static boundary conditions

zone face apply velocity-x 0 range group 'west' ;or 'east' 
zone face apply velocity-y 0 range group 'south' or 'north'
zone face apply velocity (0,0,0) range group 'bottom'

Since staged construction is to be simulated, all zones are first nulled using the command

zone null

The staged construction is simulated by reassigning the nulled zones to non-null materials. In this example, for simplicity, all materials are assumed to follow the Mohr-Coulomb model during the construction.

The starter dam is constructed first by assigning the “StarterDam” group to the Mohr-Coulomb model and setting the material properties. The command zone initialize-stress is used to initialize the vertical stress to be gravitational and the horizontal stress to be a factor of the vertical stress. Use the command model solve-static convergence 1 to cycle through each staged construction. The scripts for this purpose are

zone cmodel assign mohr-coulomb range group 'StarterDam'
zone property density 2.0 bulk 3e5 shear 2e5 range group 'StarterDam'
zone property cohesion 20 friction 36 tension 1e20 range group 'StarterDam'
zone initialize-stresses ratio 0.6 range group 'StarterDam'
;
model solve-static convergence 1

After start dam construction. The displacement contour is shown in Figure 3.

Figure 3: Displacement contour after construction of the starter dam.

For tailings and dikes, similar procedures can be adopted. Alternatively, a FISH function called “add_layer” is used for this purpose. The construction sequence is “StarterDam”, “Talings01”, “Talings02”, “Talings03”, “Dike1”, “Taings04”, “Dike2”, “Tailing05”, “Dike3”, “Tailing06”, “Dike4”, “Tailing07”, “Dike5”, “Tailing08”, “Dike6”, “Tailings09”, “Dike7”, and “Tailings10”. For example, Figure 4 shows the displacement contour after construction of Dike3.

Figure 4: Displacement contour after construction of Dike3.

After construction completion. The displacement contour is shown in Figure 5.

Figure 5: Displacement contour after construction completion.

Phreatic Surface Calculation

So far, no fluid has been considered. The stable phreatic surface needs to be calculated. This simulation is scripted in the data file fluid.dat.

After restoring the previous save file of the construction, the zone displacement, velocity, and plastic state are all re-initialized using the following commands:

zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0

Fluid properties now need to be input, such as “porosity,” “hydraulic-conductivity,” and fluid density. Note that “hydraulic-conductivity” here refers to Darcy’s coefficient, measured in units of [length/time]. It is different from another fluid property called “mobility-coefficient”, which is the value of “hydraulic-conductivity” divided by the fluid’s unit weight.

zone fluid property porosity 0.3 hydraulic-conductivity 1e-9
zone fluid-density 1

The saturation is initialized to be zero, and the saturation model is set to “cutoff” with the default cutoff value of zero:

zone gridpoint initialize saturation 0
; satuation model is set to cutoff 
zone fluid unsaturated cutoff

The fluid boundary conditions assume a fixed water table at the ground level on the “west” side and zero pore pressure at the top surface. Because the model base is at zero elevation (datum = 0), we can apply a hydraulic head of 57.75 m on the “west” side:

zone face apply head 57.75 range group 'west'
zone face apply pore-pressure 0 range group 'top'

To calculate a steady-state fluid solution, use the command:

zone fluid steady-state

Since an unsaturated flow model is present, the nonlinear effects are handled via Picard iteration. Note that since there is no transient response, it is not necessary to specify a fluid modulus (or Biot modulus) in the model. This command does not specify “effective on,” which means the default “effective off” is used, so the solution will update the pore pressures but not the total stresses in the zones, resulting in a change of effective stress. The saturation contour and pore-pressure contour at the steady-state fluid solution are presented in Figure 6 and Figure 7, respectively.

Figure 6: Saturation contour at the steady-state fluid solution.

Figure 7: Pore-pressure contour at the steady-state fluid solution.

The above command is for fluid calculation only. A mechanical-only solution is required to cycle the model to mechanical equilibrium:

model solve-static convergence 1

Static Liquefaction

The model has completed staged construction, steady-state fluid calculation, and mechanical equilibrium. To analyze the potential for static liquefaction, the Mohr-Coulomb model for tailings is not sufficient because it does not incorporate the state parameter and thus cannot distinguish between dilative and contractive behaviors using one set of material constants. In this model, NorSand is used for tailings. This is implemented in the data file “toNorSand.dat” (see Figure 8).

Figure 8: NorSand model for tailings and Mohr-Coulomb model for starter dam and dikes.

After the tailings have been changed to the NorSand model, the NorSand properties need to be input. The NorSand properties for tailings are summarized in Table 1.

Table 1: NorSand Properties for Tailings

Property	Symbol	Value
shear-reference	\(G_{ref}\)	2e4
exponent	\(m\)	0.72
pressure-reference	\(p_{ref}\)	100
poisson	\(\nu\)	0.15
critical-state-1	\(C_1\)	0.74
critical-state-2	\(C_2\)	0.05
critical-state-3	\(C_3\)	0.0
ratio-critical	\(M_{tc}\)	1.2
factor-dilatancy	\(\chi_{tc}\)	7.2
factor-coupling	\(N\)	0.25
hardening-0	\(H_0\)	45
hardening-y	\(H_y\)	30
flag-inner		on

In this example, the property “flag-inner” is turned “on” so that the inner cap yielding is considered. Note that the property “state-parameter-initial” (or “void-initial”) is significant to the material behavior. Refer to the documentation of the NorSand model for more details. In this model, for demonstration purposes, all tailings are set with an initial state parameter of 0.1, except for the first two tailing layers (“Tailings01” and “Tailing02”), which are set to -0.1 (Figure 8).

Figure 9: Initial state parameters of NorSand for tailings.

Note that the NorSand model is stress-dependent because the moduli are a function of the current stress state. The previously calculated stress state is transferred to NorSand when the model is first assigned by a FISH function called “ini_stress.” This function uses the keyword “operator” instead of “define” and includes an argument for a specific zone (“zp”) so that it can run multithreadedly. Another argument (“modelName”) is included to designate the specific constitutive model:

fish operator ini_stress(modelName,zp)
  if zone.model(zp) = modelName then
    zone.prop(zp,'stress-xx-initial') = zone.stress.eff(zp)->xx
    zone.prop(zp,'stress-yy-initial') = zone.stress.eff(zp)->yy
    zone.prop(zp,'stress-zz-initial') = zone.stress.eff(zp)->zz
    zone.prop(zp,'stress-xy-initial') = zone.stress.eff(zp)->xy
    zone.prop(zp,'stress-yz-initial') = zone.stress.eff(zp)->yz
    zone.prop(zp,'stress-xz-initial') = zone.stress.eff(zp)->xz
  endif
end
[ini_stress("norsand", ::zone.list)]

After the NorSand model and properties for tailings are assigned, run the model to reach mechanical equilibrium.

The next analysis steps (as in the datafile “static.dat”) is to evaluate if the tailing dam is potential to static liquidation due to a relatively quick surcharge on the current tailing surface. A surcharge can be from additional tailing construction upon the current tailing surface.

First, reinitialize the zone displacement, velocity, and plastic state. Second, change the model to a relatively “undrained” condition. This requires two steps: (1) turn the pore-pressure generation to “on” with the command

zone fluid property pore-pressure-generation on

and (2) set the water bulk modulus to a non-zero value

zone fluid property fluid-modulus 2e5

Note that since the water is not confined on all boundaries, a compromised (reduced) water bulk modulus of 2e5 kPa is set, which is 1/10 of the real water bulk modulus of 2e6 kPa under perfectly confined conditions.

In this model, a “ramp” FISH function is used to gradually add the surcharge. This ramp function increases linearly from 0 to 1 in 6000 cycles steps and remains 1 thereafter:

[global iniStep = mech.step]
fish define ramp
  global cycles = mech.step - iniStep
  if cycles < 6000
    ramp = float(cycles)/float(6000)
  else
    ramp = 1.0
  endif
end

The surcharge is applied using the “zone face apply stress-normal” command, followed by a value of -500 multiplied by the ramp function, and applied to the face group “Load,” which is assigned in the data file sketch.dat:

; Apply surcharge
zone face apply stress-normal -500 fish ramp range group 'Load'

The model is cycled to 12000 steps. Figure 10 presents the shear band by plotting the contour of maximum shear strain increment. Figure 11 plots the contour of displace. Both figures clearly demonstrate sliding failure due to static liquefaction. The displacement history at the highest dike top with coordinates (61,0,60.5) is plotted in Figure 12, which shows an accelerating increase, indicating an unstable slope.

The stress paths of selected zones (see golden dots in Figure 10), as well as the critical state line in the \(q\) vs \(p'\) plane are plotted in Figure 13. All these selected zones are located within the shear bands, and their final shear strengths are lower than their initial strengths, which are clear indicators of static liquefaction.

Figure 10: Maximum shear strain contour indicating localized shear bands induced by static liquefaction.

Figure 11: Displacement contour indicating sliding failure induced by static liquefaction.

Figure 12: Displacement history at the highest dike top with coordinates (61,0,60.5).

Figure 13: Stress paths of \(q\) vs \(p'\) of selected zones.

Discussions

For this model with strain-softening material behavior, the results are grid-dependent. This grid-dependence can be checked by increasing or decreasing the global FISH variable “zone_size.” The shear band width shown in Figure 10 will increase or decrease accordingly.

When a shear band occurs, the model’s global behavior is what truly matters. This is because the stress-strain behavior at specific locations can vary, including plastic strain softening, elastic unloading, and even elastic or plastic loading. Only zones within the shear bands may exhibit partial undrained failure paths, which may be mixed with elastic unloading or reloading (see Figure 13). Although the entire model assumes an “undrained” setting, maintaining a strictly undrained condition (constant volume for the soil matrix) is usually not possible for individual zones. This is due to two main reasons: (1) the model is not fluid-confined due to the presence of the phreatic surface, and (2) excess pore pressures can be transferred to neighboring zones.

Once the model is in sliding failure, the absolute displacement magnitude is less meaningful because it will continue to increase with further cycling (see Figure 12). However, the accelerating rate of displacement is a good indicator of model instability.

Data File

master.dat

model new
[global zone_size = 1.0]
program call "sketch.dat"
program call "construction.dat"
program call "fluid.dat"
program call "toNorSand.dat"
program call "static.dat"
program call "check.dat"

construction.dat

model restore 'sketch'
;
zone face skin
;
model configure fluid-flow
model large-strain off
model gravity 9.81
;
zone face apply velocity-x 0 range group 'west' ;or 'east' 
zone face apply velocity-y 0 range group 'south' or 'north'
zone face apply velocity (0,0,0) range group 'bottom'
; Null all first
zone null
; Start with just the starter dam
; units in kPa
zone cmodel assign mohr-coulomb range group 'StarterDam'
zone property density 2.0 bulk 3e5 shear 2e5 range group 'StarterDam'
zone property cohesion 20 friction 36 tension 1e20 range group 'StarterDam'
zone initialize-stresses ratio 0.6 range group 'StarterDam'
;
model solve-static convergence 1
model save 'stage-1'
;
; now build up dam layer by layer
[scount_ = 1]
fish define add_layer(mat_,type_)
  scount_ += 1
  if type_ = 'dike' then
    command
      zone cmodel assign mohr-coulomb range group [mat_] 
      zone property density 2.5 bulk 2e5 shear 1.5e5 range group [mat_]
      zone property cohesion 25 friction 32 range group [mat_] 
      zone initialize-stresses ratio 0.6 range group [mat_]
    endcommand
  endif
  if type_ = 'tails' then
    command
      zone cmodel assign mohr-coulomb range group [mat_] 
      zone property density 2.5 bulk 1e5 shear 0.6e5 range group [mat_]
      zone property cohesion 0  friction 28 range group [mat_] 
      zone initialize-stresses ratio 0.8 range group [mat_]
    endcommand
  endif
  command
    model solve-static convergence 1
    ;If saving construction stages
    model save ['stage-' + string(scount_)]
  endcommand
end
; Simulate the construction stages
[add_layer('Tailings01','tails')]
[add_layer('Tailings02','tails')]
[add_layer('Tailings03','tails')]
[add_layer('Dike1','dike')]
[add_layer('Tailings04','tails')]
[add_layer('Dike2','dike')]
[add_layer('Tailings05','tails')]
[add_layer('Dike3','dike')]
[add_layer('Tailings06','tails')]
[add_layer('Dike4','dike')]
[add_layer('Tailings07','tails')]
[add_layer('Dike5','dike')]
[add_layer('Tailings08','tails')]
[add_layer('Dike6','dike')]
[add_layer('Tailings09','tails')]
[add_layer('Dike7','dike')]
[add_layer('Tailings10','tails')]

model save 'construction'

fluid.dat

model restore 'construction'
;
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0
;
zone fluid property porosity 0.3 hydraulic-conductivity 1e-9
zone fluid-density 1
;
zone gridpoint initialize saturation 0
; satuation model is set to cutoff 
zone fluid unsaturated cutoff
;
zone face apply head 57.75 range group 'west'
zone face apply pore-pressure 0 range group 'top'
; 
zone fluid steady-state
model save 'fluid-only'
; now mechanical
model solve-static convergence 1
model save 'fluid'

toNorSand.dat

model restore "fluid"
;
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0
;
zone nodal-mixed-discretization on
; change to NorSand model
zone group 'Tailing' range group 'Tailings01' or 'Tailings02' or ...
    'Tailings03' or 'Tailings04' or 'Tailings05' or 'Tailings06' or ... 
    'Tailings07' or 'Tailings08' or 'Tailings09' or 'Tailings10'
;
zone cmodel assign norsand range group 'Tailing'
zone property shear-reference 2e4 exponent 0.72 range cmodel 'norsand'
zone property pressure-reference 100 poisson 0.15 range cmodel 'norsand'
zone property critical-state-1 .74 critical-state-2 .05 range cmodel 'norsand'
zone property ratio-critical 1.2 range cmodel 'norsand'
zone property factor-coupling 0.25 factor-dilatancy 7.2 range cmodel 'norsand'
zone property hardening-0 45.0 hardening-y 350.0 range cmodel 'norsand' 
zone property flag-inner on range cmodel 'norsand'
zone property state-parameter-initial  0.10 range cmodel 'norsand'
zone property state-parameter-initial -0.10 range group 'Tailings01'
zone property state-parameter-initial -0.10 range group 'Tailings02'
;
fish operator ini_stress(modelName,zp)
  if zone.model(zp) = modelName then
    zone.prop(zp,'stress-xx-initial') = zone.stress.eff(zp)->xx
    zone.prop(zp,'stress-yy-initial') = zone.stress.eff(zp)->yy
    zone.prop(zp,'stress-zz-initial') = zone.stress.eff(zp)->zz
    zone.prop(zp,'stress-xy-initial') = zone.stress.eff(zp)->xy
    zone.prop(zp,'stress-yz-initial') = zone.stress.eff(zp)->yz
    zone.prop(zp,'stress-xz-initial') = zone.stress.eff(zp)->xz
  endif
end
[ini_stress("norsand", ::zone.list)]

model solve-static convergence 1
model save "tonorsand"

static.dat

model restore "tonorsand"
;
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0
; 
zone fluid property pore-pressure-generation on
zone fluid property fluid-modulus 2e5
;
[global iniStep = mech.step]
fish define ramp
  global cycles = mech.step - iniStep
  if cycles < 6000
    ramp = float(cycles)/float(6000)
  else
    ramp = 1.0
  endif
end
; Apply surcharge
zone face apply stress-normal -500 fish ramp range group 'Load'
zone face apply-remove head range group "West"
;
fish define locations
  local v1 = ( -2.5, 0.5, 37.0)
  local v2 = ( 15.5, 0.5, 24.5)
  local v3 = ( 66.3, 0.5, 15.6)
  global z1 = zone.near(v1)
  global z2 = zone.near(v2)
  global z3 = zone.near(v3)
  data.scalar.value(data.scalar.create(v1)) = 1.0
  data.scalar.value(data.scalar.create(v2)) = 1.0
  data.scalar.value(data.scalar.create(v3)) = 1.0
end
[locations]
;
fish define p1
  local sten1 = zone.stress.eff(z1)
  global p1 = tensor.trace(sten1) / (-3.0)
  global q1 = math.sqrt(3.0*tensor.j2(sten1))
  local sten2 = zone.stress.eff(z2)
  global p2 = tensor.trace(sten2) / (-3.0)
  global q2 = math.sqrt(3.0*tensor.j2(sten2))
  local sten3 = zone.stress.eff(z3)
  global p3 = tensor.trace(sten3) / (-3.0)
  global q3 = math.sqrt(3.0*tensor.j2(sten3))
end
history interval 10
fish history name 'cycles' cycles
zone history name='disp-top' displacement position (61,0,60.5)
fish history name 'p1' p1
fish history name 'q1' q1
fish history name 'p2' p2
fish history name 'q2' q2
fish history name 'p3' p3
fish history name 'q3' q3
;
model solve-static cycles 10000
;
history export 'q1' vs 'p1' table 'z1'
history export 'q2' vs 'p2' table 'z2'
history export 'q3' vs 'p3' table 'z3'
table 'CSL' add (0,0) (800,960)
;
model save "static"

Endnote

[1]

This example is listed under the present document title in the i Examples dialog (use the menu command Help -> Examples… to access.