Simulating excavation and refill of full blocks using CoupledVarThresholdElementSubdomainModifier

[jmeier]

Hi all,

For various tasks in the field of engineering and mining, the excavation (and possibly refilling) of defined volumes must be modelled. Examples are excavation pits or the excavation of a tunnel. The discussions here at https://github.com/idaholab/moose/discussions contain several corresponding questions and hints (e.g. #23677, #20601, #16387, #25264, and others), on how this can be achieved with Moose.

To the best of my current knowledge, there are three options for modelling the removal of model parts in Moose:

• Option 1: The element is not really removed from the model, but a reduced stiffness and weight is assigned. An example model that uses this approach can be found here: coal_mining
• Option 2: A restart analysis is performed. Strictly speaking, a series of Moose models is calculated in which the subsequent model has an adapted mesh and imports the initial state from the predecessor model. As this approach "breaks down" the model definition into individual models, I have not pursued this approach any further to date.
• Option 3: Elements are moved from an "active" subdomain to an "inactive" subdomain during a Moose calculation. I would like to pursue this option with this "show and tell".

First of all: Many thanks to @GiudGiud, @Kavan-Khaledi, and @AdrienWehrle for their valuable advices, without which this show and tell would not have been possible.

Screenshot: excavated block at t=0.4

The code below creates a rectangular volume ("Box1") in the centre of which a small cube ("Box2") is to be removed at t=0.4. To keep the model simple, both Box1 and Box2 are linearly elastic. The approach can be described very briefly as follows:

• In addition to the subdomains Box1 and Box2, in which the corresponding elements are located, an empty subdomain 3 is created.
• Subdomain 3 contains a greatly reduced "physics" so that the elements moved there are "inactive".
• The tensor mechanics master block is therefore only created for the subdomains of Box1 and Box2. Gravity also only works here.
• In subdomain 3, a diffusion kernel is only created as a "dummy" kernel.
• When Box2 is moved to subdomain 3 is controlled by a function (PiecewiseConstant). The value of this function is stored in an aux variable. This aux variable is in turn used by a CoupledVarThresholdElementSubdomainModifier to move the elements of Box2 to subdomain 3 at t=0.4.

In the model initially Box1 and Box2 are active. At t=0.4 Box2 gets "removed". As expected, as a result of the unloading, the material of Box1 lifts in the proximity of the removed Box2.

To run the input file, your Moose app needs the modules "MISC" (for the gravity) and "TENSOR_MECHANICS".

As a side-note: I presume that in productive models (especially with more complex material behaviour) a mixture of option 1 (reduction of stiffness and weight) and option 3 (subdomain modifier) is probably the method of choice. First make it soft and light and then move it to the inactive subdomain.

Code Block: Moose Input File

# testing construction stages
# primary model units (m | s | kg --> N/m²)

#some model-specific variables
material_density = 21.5E3         # kg/m^3
gravitational_acceleration = 9.81 # m/s^2
sig_top = 0.0                     # N/m^2
z_top = 10.0

[Problem]
  kernel_coverage_check = false
  material_coverage_check = false
[]

Box2_inactive = '3'
inactive_domain_blocks = ${Box2_inactive}

[Mesh]

  [BaseMesh]
    type = GeneratedMeshGenerator
    elem_type = "TET10"
    dim = 3
    nx = 20
    ny = 20
    nz = 2
    xmin = -10
    xmax = +10
    ymin = -10
    ymax = +10
    zmin = -2
    zmax = +2
  []

  [Box1]
    type = SubdomainBoundingBoxGenerator
    input = "BaseMesh"
    block_id = 1
    location = "INSIDE"
    bottom_left = "-20 -20 -2"
    top_right   = "+20 +20 +2"
  []

  [Box2]
    type = SubdomainBoundingBoxGenerator
    input = "Box1"
    block_id = 2
    location = "INSIDE"
    bottom_left = "-2 -2 +2"
    top_right   = "+2 +2 0"
  []

  add_subdomain_ids = ${inactive_domain_blocks}

[]

active_domain_blocks =  '1 2' #Box1 Box2

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  large_kinematics = false
[]

# ===== Tensor mechanics for all active domain blocks =====
[Modules/TensorMechanics/Master]
  [./all]
    new_system = true
    formulation = total
    strain = small

    block = ${active_domain_blocks}
    add_variables = true

    generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz '
  []
[]

# ===== Gravity =====
[Kernels]
  [./Gravity]
    type = Gravity
    block = ${active_domain_blocks}
    use_displaced_mesh = false
    variable = disp_z
    value = -${gravitational_acceleration}
  []
 [./donothing1]
    type = MatDiffusion
    block = ${inactive_domain_blocks}
    variable = disp_z
    diffusivity = 1e-7
  []
 [./donothing2]
    type = MatDiffusion
    block = ${inactive_domain_blocks}
    variable = disp_x
    diffusivity = 1e-7
  []
 [./donothing3]
    type = MatDiffusion
    block = ${inactive_domain_blocks}
    variable = disp_y
    diffusivity = 1e-7
  []
[]

# ===== Boundary Conditions: Fixies at XMin, XMax, YMin, YMax, ZMin =====
[BCs]
  
  [./BoundaryXMin_DispX]
    type = DirichletBC
    boundary = 'left'
    variable = disp_x
    value = 0
  []

  [./BoundaryXMax_DispX]
    type = DirichletBC
    boundary = 'right'
    variable = disp_x
    value = 0
  []
  
  [./BoundaryYMin_DispY]
    type = DirichletBC
    boundary = 'bottom'
    variable = disp_y
    value = 0
  []
  
  [./BoundaryYMax_DispY]
    type = DirichletBC
    boundary = 'top'
    variable = disp_y
    value = 0
  []
  
  [./BoundaryZMin_DispX]
    type = DirichletBC
    boundary = 'back'
    variable = disp_x
    value = 0
  []

  [./BoundaryZMin_DispY]
    type = DirichletBC
    boundary = 'back'
    variable = disp_y
    value = 0
  []

  [./BoundaryZMin_DispZ]
    type = DirichletBC
    boundary = 'back'
    variable = disp_z
    value = 0
  []

[]

# ===== Initial Conditions: Initial Stress Field =====
[Functions]
  [ini_xx]
    type = ParsedFunction
    expression = '(-sig_top - rho * g * (z_top - z)) * K0'
    symbol_names  = 'sig_top     z_top      rho                   g                               K0  '
    symbol_values = '${sig_top}  ${z_top}   ${material_density}   ${gravitational_acceleration}   0.3 '
  []
  [ini_yy]
    type = ParsedFunction
    expression = '-sig_top - rho * g * (z_top - z)'
    symbol_names  = 'sig_top     z_top      rho                   g                               K0  '
    symbol_values = '${sig_top}  ${z_top}   ${material_density}   ${gravitational_acceleration}   0.3 '
  []
  [ini_zz]
    type = ParsedFunction
    expression = '(-sig_top - rho * g * (z_top - z)) * K0'
    symbol_names  = 'sig_top     z_top      rho                   g                               K0  '
    symbol_values = '${sig_top}  ${z_top}   ${material_density}   ${gravitational_acceleration}   0.3 '
  []
[]

# ===== Materials (linear-elastic to keep it simple) =====
[Materials]

  [./volAll_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = ${active_domain_blocks}
    youngs_modulus = 450E6 # 450 MN/m^2
    poissons_ratio = 0.15
  []

  [./volAll_stress]
    type = ComputeLagrangianWrappedStress
    block = ${active_domain_blocks}
  []

  [./volAll_stress_base]
    type = ComputeLinearElasticStress
    block = ${active_domain_blocks}
  []

  [./volAll_ini_stress]
    type = ComputeEigenstrainFromInitialStress
    eigenstrain_name = ini_stress
    initial_stress = 'ini_xx 0 0  0 ini_yy 0  0 0 ini_zz'
    block = ${active_domain_blocks}
  []

  [./density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = ${material_density}
    block = ${active_domain_blocks}
  []

[]

# ===== Block Activation/Deactivation =====
[Functions]
  # this function at which time (x-value) the block is active (y=1) or inactive (y=0)
  [Box2_ActivationStateFunction]
    type = PiecewiseConstant
    direction = left
    xy_data = '0.0  0
               0.2  0
               0.4  1
               1.0  1'
  []
[]
# aux variable to store the value of Box2_ActivationStateFunction
[AuxVariables]
  [Box2_Activation_AuxVariable]
    family = MONOMIAL
    order = CONSTANT
  []
[]
# link Box2_ActivationStateFunction and Box2_Activation_AuxVariable
[AuxKernels]
  [Box2_Activation_FunctionAux]
    type = FunctionAux
    variable = 'Box2_Activation_AuxVariable'
    function = 'Box2_ActivationStateFunction'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]
# move Box2 between subdomains
[UserObjects]
  [Box2_SubdomainModifier]
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'Box2_Activation_AuxVariable'
    criterion_type = ABOVE
    threshold = 0.5
    block = 2 #'Box2'
    subdomain_id = 2 #Box2
    complement_subdomain_id = ${Box2_inactive}
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

# ===== Executioner =====
[Executioner]
  type = Transient

  end_time = 1.0
  [TimeSteppers]
    [BlockEventTimeStepper]
      type = TimeSequenceStepper
      time_sequence = '0.0 0.2 0.4 1.0'
    []
  []
  
  solve_type = 'PJFNK'
  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'

  nl_abs_tol = 1E-3
  nl_max_its = 400
  
  l_tol = 1E-3
  l_max_its = 200
[]

[Outputs]
  perf_graph = true
  exodus = true
[]

[PEI0214]

Hello, I would like to ask how to realize the filling process? Thank you for your help.

[PEI0214]

When simulating the fill process, a non-convergence occurs when moving elements from a subdomain “without physics” to a subdomain “with physics”. Is this due to mesh distortion? Do you have any good solutions?

[jmeier]

When moving elements from one subdomain to another, you may very well want to ramp in or ramp out various material parameters over time. Which parameters these are depends very much on what you want to model. In solid mechanics, however, density is often one such parameter. As a rule of thumb: try to model sufficiently well what happens in reality (e.g. a water dam is usually not put down "all at once" - rather in many small steps).

[PEI0214]

Do you have examples of simulated fill construction? If so, can you provide it for learning?Thank you.

[PEI0214]

How do you implement the "Option 2"?

[jmeier]

Do you have examples of simulated fill construction? If so, can you provide it for learning?

Currently I have no example showing this.

How do you implement the "Option 2"?

Option 1 and 2 are not shown in this example.