Simulating underground coal mining

Introduction

In an underground longwall coal mine, coal is mined in "panels", as shown in Figure 1. These panels are typically 3–4m in height, 150–400m in width and 1000–4000m long. Coal is extracted from one end, moving towards the other end. As mining progresses, a void is created behind the mining machinery where the coal used to be. The rock material above this void— the "overburden"—is not strong enough to support itself, so it collapses downwards. The resulting partial void and collapsed material is called the "goaf" (or "gob").

Figure 1: Pictorial representation of a single panel in longwall mining. The dimensions and process are described in the text. Figure sourced from citizensagainstlongwallmining.org who obtained it from the Energy Information Administration.

There are many geotechnical aspects that are of interest in this process, but this Example concentrates on the following.

  1. The vertical displacement of the ground surface due to the longwall mining. This is called the "subsidence". This obviously has effects on buildings and other man-made structures like roads, but it can also change surface water pathways, affect vegetation, and cause surface-water ponding, which all have effects on the local ecology.

  2. The vertical deformation within the overburden.

  3. The patterns of fracture in the overburden. Together with the previous item, this can have effects on operational aspects of longwall mining, such as the load patterns on the "chocks" that support the roof of the goaf next to the mining machinery. It also strongly effects the flow of fluids in the porous rocks that make up the overburden, as discussed next.

  4. The flow of water through the rock. If the fracturing and deformation is extensive then fluids can easily move through the rock. This can result in mines flooding with water, or causing excessive drawdown of groundwater aquifers which can effect other users of groundwater (such as irrigators for farming) or can effect the rates of groundwater baseflow to surrounding river systems which may effect local and not-so-local ecosystems. The deformation and fracture of the overburden can also lead to aquifer mixing, where an aquifer consisting of "dirty" water (with high salinity, for example) pollutes a nearby aquifer of clean water. Although not studied in this Example, the deformation and fracture can also lead to large releases of methane gas from overlying (and underlying) coal seams. This methane can then flow through the highly permeable fractured rock system to the atmosphere, resulting in large greenhouse gas emissions. Or it can flow to the mine workings which is extremely hazardous in terms of mine fires and explosions.

This Example does not seek to build a realistic model of a specific mine. Instead, it explains how the PorousFlow module can be used to build such a model, by exploring a simple 3D model.

The model

Two input files are provided with this example: coarse_with_fluid.i and fine_with_fluid.i. They essentially only differ in the input mesh. The results below are drawn from the fine model, but it takes a few hours to complete on a computer cluster, so the coarse model might be a better model to explore initially.

Geometry

Figure 2 shows the geometrical setup as well as the fine mesh. The model represents a single longwall panel of total width 300m and length 1000m. Only half of the width is modelled, and appropriate boundary conditions used to simulate the full panel. In the fine model, the panel is meshed with 10m square elements.

Figure 2: The geometry of the model and the mesh used. A single longwall panel of total width 300m and length 1000m is represented by the region with a fine mesh.

Fluid mechanics

Single-phase unsaturated fluid is used with the nonlinear Variable being the fluid porepressure. Full coupling with mechanical deformations is used, so the Kernels read

  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = porepressure
  []
  [flux]
    type = PorousFlowAdvectiveFlux
    use_displaced_mesh = false
    variable = porepressure
    gravity = '0 0 -10E-6'
    fluid_component = 0
  []
  [poro_vol_exp]
    type = PorousFlowMassVolumetricExpansion
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    variable = porepressure
    fluid_component = 0
  []
[]
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

All stresses are measured in MPa and the time unit is years. A Biot coefficient of 0.7 is used, and the fluid properties are:

[Modules]
  [FluidProperties]
    [simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E3
      density0 = 1000
      thermal_expansion = 0
      viscosity = 3.5E-17
    []
  []
[]
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

A van Genuchten capillary pressure relationship is used

  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1 # MPa^-1
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

A Corey relative permeability is used

  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.4
    sum_s_res = 0.4
    phase = 0
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

The fluid porepressure is fixed at the outer boundaries, and a PorousFlowPiecewiseLinearSink with a spatially and temporally varying flux_function is used to withdraw water from the rock into the goaf:

  [fix_porepressure]
    type = FunctionDirichletBC
    variable = porepressure
    boundary = 'ymin ymax xmax'
    function = ini_pp
  []
  [roof_porepressure]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    pt_vals = '-1E3 1E3'
    multipliers = '-1 1'
    fluid_phase = 0
    flux_function = roof_conductance
    boundary = roof
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

In multi-phase scenarios, such a simple representation of the goaf boundary condition is not valid.

Solid mechanics

Layered Cosserat solid mechanics is employed, meaning there are 3 translational and 2 rotational degrees of freedom. A quasi-static approximation is used, so the solid-mechanical Kernels are

  [cx_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [cy_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [cz_elastic]
    type = CosseratStressDivergenceTensors
    use_displaced_mesh = false
    variable = disp_z
    component = 2
  []
  [x_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_x
    displacements = 'wc_x wc_y wc_z'
    component = 0
    base_name = couple
  []
  [y_couple]
    type = StressDivergenceTensors
    use_displaced_mesh = false
    variable = wc_y
    displacements = 'wc_x wc_y wc_z'
    component = 1
    base_name = couple
  []
  [x_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_x
    component = 0
  []
  [y_moment]
    type = MomentBalancing
    use_displaced_mesh = false
    variable = wc_y
    component = 1
  []
  [gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -10E-6 # remember this is in MPa
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_x
    component = 0
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    variable = disp_y
    component = 1
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    use_displaced_mesh = false
    component = 2
    variable = disp_z
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Quite a complicated elasto-plastic model is used and this accounts for the majority of nonlinear iterations of this model, as well as the relative slowness of computing the residual and Jacobian. The elasticity tensors (one for standard, one for Cosserat) are computed via:

  [elasticity_tensor_0]
    type = ComputeLayeredCosseratElasticityTensor
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
  []
  [elasticity_tensor_1]
    type = ComputeLayeredCosseratElasticityTensor
    block = 1
    young = 8E3 # MPa
    poisson = 0.25
    layer_thickness = 1.0
    joint_normal_stiffness = 1E9 # huge
    joint_shear_stiffness = 1E3 # MPa
    elasticity_tensor_prefactor = excav_sideways
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

where the excav_sideways Function simulates the excavation of coal (it is 1 ahead of the coal face and zero behind it). Small strain is used, and an insitu stress that increases with depth is assumed:

  [strain]
    type = ComputeCosseratIncrementalSmallStrain
    eigenstrain_names = ini_stress
  []
  [ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_xx 0  0 0 ini_zz'
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Capped Mohr-Coulomb plus weak-plane plasticity is used and the plastic models are cycled

  [stress_0]
    type = ComputeMultipleInelasticCosseratStress
     block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    inelastic_models = 'mc wp'
    cycle_models = true
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [stress_1]
    type = ComputeMultipleInelasticCosseratStress
    block = 1
    inelastic_models = ''
    relative_tolerance = 2.0
    absolute_tolerance = 1E6
    max_iterations = 1
    tangent_operator = nonlinear
    perform_finite_strain_rotations = false
  []
  [mc]
    type = CappedMohrCoulombCosseratStressUpdate
    warn_about_precision_loss = false
    host_youngs_modulus = 8E3
    host_poissons_ratio = 0.25
    base_name = mc
    tensile_strength = mc_tensile_str_strong_harden
    compressive_strength = mc_compressive_str
    cohesion = mc_coh_strong_harden
    friction_angle = mc_fric
    dilation_angle = mc_dil
    max_NR_iterations = 100000
    smoothing_tol = 0.1 # MPa  # Must be linked to cohesion
    yield_function_tol = 1E-9 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0
  []
  [wp]
    type = CappedWeakPlaneCosseratStressUpdate
    warn_about_precision_loss = false
    base_name = wp
    cohesion = wp_coh_harden
    tan_friction_angle = wp_tan_fric
    tan_dilation_angle = wp_tan_dil
    tensile_strength = wp_tensile_str_harden
    compressive_strength = wp_compressive_str_soften
    max_NR_iterations = 10000
    tip_smoother = 0.05
    smoothing_tol = 0.05 # MPa  # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
    yield_function_tol = 1E-11 # MPa.  this is essentially the lowest possible without lots of precision loss
    perfect_guess = true
    min_step_size = 1.0E-3
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

The plastic moduli are:

  [mc_coh_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.99 # MPa
    value_residual = 2.01 # MPa
    rate = 1.0
  []
  [mc_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.61 # 35deg
  []
  [mc_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.15 # 8deg
  []

  [mc_tensile_str_strong_harden]
    type = TensorMechanicsHardeningExponential
    value_0 = 1.0 # MPa
    value_residual = 1.0 # MPa
    rate = 1.0
  []
  [mc_compressive_str]
    type = TensorMechanicsHardeningCubic
    value_0 = 100 # Large!
    value_residual = 100
    internal_limit = 0.1
  []

  [wp_coh_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_tan_fric]
    type = TensorMechanicsHardeningConstant
    value = 0.26 # 15deg
  []
  [wp_tan_dil]
    type = TensorMechanicsHardeningConstant
    value = 0.18 # 10deg
  []

  [wp_tensile_str_harden]
    type = TensorMechanicsHardeningCubic
    value_0 = 0.05
    value_residual = 0.05
    internal_limit = 10
  []
  [wp_compressive_str_soften]
    type = TensorMechanicsHardeningCubic
    value_0 = 100
    value_residual = 1
    internal_limit = 1.0
  []
[]
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Roller boundary conditions are prescribed at the boundaries. At the roof of the excavation, a StickyBC is employed to prevent the roof from collapsing further than 3m:

  [roof_bcs]
    type = StickyBC
    variable = disp_z
    min_value = -3.0
    boundary = roof
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Coupling

Full coupling between the fluid and solid mechanics is evident in the Kernels. The PorousFlowDictator ensures all nonzero Jacobian entries are computed:

  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Other aspects of the fluid-solid coupling are the porosity:

  [porosity_bulk]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    ensure_positive = true
    porosity_zero = 0.02
    solid_bulk = 5.3333E3
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

and the permeability:

  [permeability_bulk]
    type = PorousFlowPermeabilityKozenyCarman
    block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
    poroperm_function = kozeny_carman_phi0
    k0 = 1E-15
    phi0 = 0.02
    n = 2
    m = 2
  []
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)

Results

Despite being a very simple model with no interesting lithology that usually induces characteristic fracture and displacement patterns, the results display many interesting features.

Figure 3: The evolution of the destressed zone where the vertical stress has less magnitude than the initial vertical stress. The zone is colored by the magnitude of de-stressing. The black rectangle shows the coal that will be excavated.

Figure 4: The evolution of the zone where the vertical displacement is cm. The zone is colored by the magnitude of vertical displacement. The black rectangle shows the coal that will be excavated.

Figure 5: The evolution of the zone where the porepressure is reduced by more than MPa. The zone is colored by the magnitude of the porepressure reduction. The green arrows show the Darcy velocity. The black rectangle shows the coal that will be excavated.

Some of interesting features are:

  • Initially the roof holds up before suddenly collapsing at around 100m of excavation.

  • After this point, there is some periodicity in the plastic strains and vertical stresses, although this may be an artificat of the mesh rather than anything physical.

  • The minimum vertical displacement is actually m in the region where the sudden collapse occurs towards the start of the panel. The StickyBC are not sufficient to prevent this. Therefore, the results towards the panel's beginning shouldn't be treated too seriously. To avoid this, MOOSE's Constraint system, or the mining-specific features CSIRO's private MOOSE-based app, could be employed.

  • There is significant rock fracture, indicated by the Mohr-Coulomb plastic strains of greater than % in the immediate 10m of the roof.

  • Mohr-Coulomb shear failure is much more prevelant than Mohr-Coulomb tensile failure throughout the model.

  • Mohr-Coulomb shear failure also occurs in the unexcavated material close to the panel (the ribs and roadways).

  • Tensile opening of the joints occurs in large regions of the roof, indicated by weak-plane tensile plasticity of greater than %.

  • Shear failure of the joints occurs throughout all regions of the model above the goaf. It is greater than 1% in the immediate m of the roof.

  • The surrounding unexcavated material experiences large vertical loads, while there is an annular region on the outside and above the excavation that is de-stressed. The material above the goaf therefore experiences de-stressing and then recompaction.

  • Porosity and permeability enhancement typically occurs in regions of highest weak-plane failure. The permeability is increased by approximately a factor of 2 in the goaf region. This is actually a tiny enhancement compared to the expected value of around and is due to the Kozeny-Carman permeability relationship which is not really valid in this setting. If readers are interested in simulating realistic mining-induced permeability enhancements they should enquire about using CSIRO's private MOOSE-based app.

  • Mining-induced drawdown of the water porepressure at the ground surface starts to occur after around m (2 weeks) of excavation. In this model, the region above the panel rapidly becomes unsaturated, with saturation reaching approximately 80%, as water flows to the goaf region.

  • Most of the high Darcy fluid velocities occur around the periphery of the mining panel.

Results along the centre-line of the panel.

The figures below show results contoured on a vertical section along the centre-line of the panel, after excavation has completed. The excavation is shown as a black line.

Mohr-Coulomb shear plastic strain on a vertical section along the centre-line of the panel

Mohr-Coulomb tensile plastic strain on a vertical section along the centre-line of the panel

Weak-plane shear plastic strain on a vertical section along the centre-line of the panel

Weak-plane tensile plastic strain on a vertical section along the centre-line of the panel

Vertical permeability on a vertical section along the centre-line of the panel

Porepressure (color) and Darcy flow vectors on a vertical section along the centre-line of the panel

Results across the panel.

The figures below show results contoured on a vertical section across the half-panel that is modelled, after excavation has completed. The excavation is shown as a black line.

Mohr-Coulomb shear plastic strain on a vertical section across the panel

Mohr-Coulomb tensile plastic strain on a vertical section across the panel

Weak-plane shear plastic strain on a vertical section across the panel

Weak-plane tensile plastic strain on a vertical section across the panel

Vertical permeability on a vertical section across the panel

Porepressure (color) and Darcy flow vectors on a vertical section across the panel