Action for Brown & Booker (1985) approach to simulate excavation

[nlosacco]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

TL;DR: successfully implemented a robust excavation procedure for geotechnical applications. Seeking for help in developing a corresponding MOOSE Action, for sake conciseness and simplicity.

---

Hi all,
as researchers in Geotechnical Engineering, we often deal with excavations (tunnels, open cuts, slopes...). I personally think the approach proposed by Brown & Booker (1985): https://www.sciencedirect.com/science/article/abs/pii/0266352X85900242 is the most reliable, and it is implemented in many FEM packages, either "geotechnics-oriented" or general purpose. The approach can be roughly summarised as follows:

1. elements representing the soil to be excavated are removed (i.e. internal variables vanish, global matrix reduces, and so forth...);
2. nodal forces are applied at the new excavation boundary to preserve the existing equilibrium;
3. the same nodal forces are reduced to zero (or possibly to a fraction 1 - $\lambda$, with $\lambda$ an "unloading ratio") linearly over a certain time.

Steps 1. and 2. occur simultaneously and instantaneously at the beginning of the excavation whereas step 3. develops in subsequent increments. If porewater is present we may add the following to 2 and 3:

2. porewater pressure values existing prior to the excavation are prescribed as a BC at the new excavation boundary;
3. porewater pressure from 2. are linearly changed over time to reach the appropriate value at the end of the excavation.

I finally managed to apply the whole procedure, please find pasted here a very simple example consisting of a fully saturated element column composed of 2 stacked 1x1 QUAD4 elements, the top element being excavated in $\Delta$ t = 0.5 starting at t = 0.5, with pore pressure = 0 at the top of the column. Element removal is simulated by moving the block to be excavated to an inactive block (i.e. with no active *Kernels or *Variables) as suggested in other posts, using MeshModifiers/TimedSubdomainModifier.

Few comments on the input file, that you might want to address:

1. I had to use Mesh/SideSetsBetweenSubdomainsGenerator instead of the newer MeshModifiers/SidesetAroundSubdomainUpdater, because I was stuck with an old MOOSE version while running the tests;
2. I couldn't use a PorousFlowBasicTHM Action as it probably tries and create AuxVariables also in the inactive blocks, thus throwing an error;
3. AuxKernels/FreezeForceAux is a simple object I created to store the value of a variable at a certain time; the same could probably be obtained through a built-in ParsedAux;
4. BCs/LockBoundaryBC does what its name suggests: it applies a DirichletBC with the value of a variable fixed at a certain time (surprisingly I could not find a built-in way to do it, but I might have missed it);
5. I had to use checkpoint = true in [Outputs] before the exodus block, otherwise, I'd get a weird ExodusII_IO_Helper: block id 1 not found in block_ids. error when using PorousFlow. No clue why it happens.

Now, as you can see, even for this tiny example the input file becomes bloated with loads of AuxVariables and AuxKernels for the excavation procedure. Furthermore, I am forced to create a dummy step (in my input from 0.5 to 1.0) to remove elements / block nodes / extract nodal forces.

Hence, I think a proper Action would be required. Ideally, it should take the following parameters:

• block to be excavated;
• new excavation boundary (it would be nice if this was automatically calculated, though);
• start time;
• duration;
• displacement variables (i.e. those for which a NodalKernel is applied during excavation);
• porewater pressure variable (and any other variables, e.g temperature, for which a DirichletBC would be appropriate during excavation);
• excavation percentage $\lambda$ (for the time being, $\lambda$ = 1.0 is ok, though);

I would still be happy if element deletion with SubdomainModifier and the relevant [Controls] blocks needed to be kept in the input file, but at least the Action should enforce instantaneous calculation and storage of values of nodal reactions and variables at the very beginning of the excavation step, simultaneously with element deletion, and subsequent change of those values.

I think such a functionality would be greatly appreciated by the ever increasing geotechnical community using MOOSE.
Hence, I'm seeking for guidance to develop the Action, as thoroughly explained.

Thanks for any help.

[GlobalParams]
    displacements = 'disp_x disp_y'
    use_displaced_mesh = false
    PorousFlowDictator = dictator
    biot_coefficient = 1
    block = '0 1'
[]

[Problem]
    solve = true
    material_coverage_check = ONLY_LIST
    material_coverage_block_list = '0 1'
    kernel_coverage_check = ONLY_LIST
    kernel_coverage_block_list = '0 1'
[]

[Mesh]
    [gen]
        type = GeneratedMeshGenerator
        elem_type = QUAD4
        dim = 2
        xmin = 0.
        xmax = 1.
        ymin = 0.
        ymax = 2.
        nx = 1
        ny = 2
        subdomain_ids = 0
        subdomain_name = 'column'
    []
    [exc]
        type = SubdomainBoundingBoxGenerator
        input = 'gen'
        block_id = 1
        block_name = 'exc'
        bottom_left = '0. 1. 0.'
        top_right = '1. 2. 1.'
    []
    [top_exc]
        input = 'exc'
        type = SideSetsBetweenSubdomainsGenerator
        primary_block = '0'
        paired_block = '1'
        new_boundary = 'top_exc'
    []

    add_subdomain_ids = '10'
    add_subdomain_names = 'excavated'
[]

[MeshModifiers]
  [GlobalSubdomainModifier]
    type = TimedSubdomainModifier
    times = '0.5'
    blocks_from = '1'
    blocks_to = '10' 
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
[]

[Physics]
    [SolidMechanics]
        [QuasiStatic]
            [all]
                strain = SMALL
                incremental = true
                add_variables = false
                eigenstrain_names = 'ini_stress'
                material_output_order = FIRST
                save_in = 'react_x react_y'
            []
        []
    []
[]

[Variables]
    [disp_x]
        family = LAGRANGE
        order = FIRST
    []
    [disp_y]
        family = LAGRANGE
        order = FIRST
    []
    [porepressure]
        family = LAGRANGE
        order = FIRST
    []
[]

# ===== Kernels: Gravity =====
[Kernels]
    [gravity]
        type = Gravity
        variable = disp_y
        value = -10.
        save_in = 'react_y'
        use_displaced_mesh = false
        enable = true
    []
[]

# ===== Kernels: PorousFlow =====
[Kernels]
    [effective_stress_x]
        type = PorousFlowEffectiveStressCoupling
        variable = 'disp_x'
        component = 0
        save_in = 'react_x'
    []

    [effective_stress_y]
        type = PorousFlowEffectiveStressCoupling
        variable = 'disp_y'
        component = 1
        save_in = 'react_y'
    []

    [flux]
        type = PorousFlowFullySaturatedDarcyFlow
        variable = 'porepressure'
        fluid_component = 0
        gravity = '0 -10 0'
    []
[]

# ===== UserObjects for porous flow analysis =====
[UserObjects]
    [dictator]
        type = PorousFlowDictator
        porous_flow_vars = 'porepressure disp_x disp_y'
        number_fluid_phases = 1
        number_fluid_components = 1
    []
[]

# ===== AuxVariable & AuxKernel: stress =====
[AuxVariables]
    [effective_mean_pressure]
        family = MONOMIAL
        order = FIRST
    []
    [stress_xx]
        order = FIRST
        family = MONOMIAL
    []
    [stress_yy]
        order = FIRST
        family = MONOMIAL
    []
    [stress_zz]
        order = FIRST
        family = MONOMIAL
    []
    [stress_xy]
        order = FIRST
        family = MONOMIAL
    []
    [stress_xz]
        order = FIRST
        family = MONOMIAL
    []
    [stress_yz]
        order = FIRST
        family = MONOMIAL
    []
    [stress_maxprincipal]
        order = FIRST
        family = MONOMIAL
    []
    [stress_minprincipal]
        order = FIRST
        family = MONOMIAL
    []
[]
[AuxKernels]
    [effective_mean_pressure]
        type = RankTwoScalarAux
        rank_two_tensor = stress
        variable = effective_mean_pressure
        scalar_type = hydrostatic
        execute_on = 'TIMESTEP_BEGIN'
    []

    [stress_xx]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 0
        index_j = 0
        variable = stress_xx
    []
    [stress_yy]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 1
        index_j = 1
        variable = stress_yy
    []
    [stress_zz]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 2
        index_j = 2
        variable = stress_zz
    []
    [stress_xy]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 0
        index_j = 1
        variable = stress_xy
    []
    [stress_yz]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 1
        index_j = 2
        variable = stress_yz
    []
    [stress_xz]
        type = RankTwoAux
        rank_two_tensor = stress
        index_i = 2
        index_j = 0
        variable = stress_xz
    []
    [stress_maxprincipal]
        type = RankTwoScalarAux
        rank_two_tensor = stress
        variable = stress_maxprincipal
        scalar_type = MaxPrincipal
    []
    [stress_minprincipal]
        type = RankTwoScalarAux
        rank_two_tensor = stress
        variable = stress_minprincipal
        scalar_type = MinPrincipal
    []
[]

# ===== AuxVariable & AuxKernel: excavation  =====
[AuxVariables]
    [react_x]
        family = LAGRANGE
        order = FIRST
    []
    [react_y]
        family = LAGRANGE
        order = FIRST
    []
    [react_base_y]
        family = LAGRANGE
        order = FIRST
    []
    [frozen_react_y]
        family = LAGRANGE
        order = FIRST
    []
    [frozen_porepressure]
        family = LAGRANGE
        order = FIRST
    []
    [relaxing_react_y]
        family = LAGRANGE
        order = FIRST
    []
    [relaxing_porepressure]
        family = LAGRANGE
        order = FIRST
    []
[]
[AuxKernels]
  [freeze_react_y]
    type = FreezeForceAux
    variable = frozen_react_y
    source_variable = react_y
    freeze_time = 0.5
    execute_on = 'timestep_end'
  []
  [freeze_pwp]
    type = FreezeForceAux
    variable = frozen_porepressure
    source_variable = porepressure
    freeze_time = 0.5
    execute_on = 'timestep_end'
  []
  [frozen_react_y]
    type = ParsedAux
    variable = frozen_react_y
    coupled_variables = 'frozen_react_y'
    expression = 'frozen_react_y'
    execute_on = 'timestep_end'
  []
  [scale_excavation_force]
    type = ParsedAux
    variable = relaxing_react_y
    coupled_variables = frozen_react_y
    use_xyzt = true
    expression = 'frozen_react_y * max(1.0 - (t - 1.0)/0.5, 0.0)'
    execute_on = 'timestep_begin'
  []
  [scale_pwp]
    type = ParsedAux
    variable = relaxing_porepressure
    coupled_variables = 'frozen_porepressure'
    use_xyzt = true
    expression = 'frozen_porepressure * max(1.0 - (t - 1.0)/0.5, 0.0)'
    execute_on = 'timestep_begin'
  []
[]

# ===== Initial Conditions: Pore-Pressure =====
# due to gravity, the initial pore pressure is hydrostatic
[Functions]
    [func_ini_porepressure]
        type = 'ParsedFunction'
        expression = '10. * (2. - y) '
    []
  [weight]
    type = ParsedFunction
    expression = '10.* ( y - 2.)' # initial EFFECTIVE stress that should result from the weight force
  []
[]

[ICs]
    [porepressure]
        type = FunctionIC
        variable = 'porepressure'
        function = 'func_ini_porepressure'
    []
[]

# ===== BCs: Mechanical =====
[BCs]
  [left]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []
  [right]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0.0
  []
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0.0
  []
  [bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

# ===== BCs: Hydraulic =====
[BCs]
  [top_pwp]
    type = DirichletBC
    variable = porepressure
    boundary = top
    value = 0.
    enable = true 
  []
[]

# ===== BCs: Excavation =====
[BCs]
    [lock_exc_top_disp]
        type = LockBoundaryBC
         variable = disp_y
         locked_variable = disp_y
         boundary = 'top_exc'
         value = 0.0
         enable = false
    []
    [lock_exc_top_pwp]
        type = LockBoundaryBC
        variable = porepressure
        locked_variable = porepressure
        boundary = 'top_exc'
        value = 0.0
        enable = false
    []
    [exc_top_pwp]
        type = FunctorDirichletBC
        variable = porepressure
        boundary = 'top_exc'
        functor = 'relaxing_porepressure'
        enable = false
    []
[]
[NodalKernels]
    [force_exc_top]
        type = CoupledForceNodalKernel
        variable = disp_y
        v = relaxing_react_y
        boundary = 'top_exc'
        enable = false
    []
[]

# ===== Material: Fluid Properties =====
[FluidProperties]
    [simple_fluid]
        type = SimpleFluidProperties
        bulk_modulus = 2.2E+6
        density0 = 1.
        thermal_expansion = 0
        viscosity = 1.E-6
    []
[]

# ===== Controls: Excavation  =====
[Controls]
  [hold]
    type = TimePeriod
    enable_objects = 'BCs::lock_exc_top_disp BCs::lock_exc_top_pwp'
    start_time = 0.5
    end_time = 1.0
    execute_on = 'timestep_begin'
  []
  [force]
    type = TimePeriod
    enable_objects = 'NodalKernels::force_exc_top BCs::exc_top_pwp'
    start_time = 1.0
    execute_on = 'timestep_begin'
  []
[]

# ===== Material: Volume Elements =====
[Materials]
    # initial stresses
    [ini_stress]
        type = ComputeEigenstrainFromInitialStress
        initial_stress = '0. 0. 0. 0. weight 0. 0. 0. 0.'
        eigenstrain_name = 'ini_stress'
    []
    [elasticity_tensor]
        type = ComputeIsotropicElasticityTensor
        youngs_modulus = 1.E+4
        poissons_ratio = 0.
    []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
    [temperature]
        type = PorousFlowTemperature
    []
    [density]
        type = GenericConstantMaterial
        prop_names = density
        prop_values = 2.
    []

    # porous flow

    [eff_fluid_pressure] #as in example undrained_oedometer.i
        type = PorousFlowEffectiveFluidPressure
    []
    [ppss] 
        type = PorousFlow1PhaseFullySaturated
        porepressure = 'porepressure'
    []
    [massfrac] 
        type = PorousFlowMassFraction
    []
    [simple_fluid] 
        type = PorousFlowSingleComponentFluid
        fp = simple_fluid
        phase = 0
    []
    [porosity_bulk] 
        type = PorousFlowPorosityConst
        PorousFlowDictator = dictator
        porosity = 0.5
    []
    [permeability]
        type = PorousFlowPermeabilityConst
        permeability = '1e-13 0 0  0 1e-13 0  0 0 1e-13' 
    []
[]

[Preconditioning]
    [SMP]
        type = SMP
        full = true

        petsc_options = '-ksp_snes_ew'
        petsc_options_iname = '-ksp_type -pc_type -pc_hypre_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels -ksp_gmres_restart'
        petsc_options_value = ' gmres     hypre    boomeramg      lu           NONZERO                   4                     301'
    []
[]

[Executioner]
    type = Transient
    solve_type = 'PJFNK'

    petsc_options = '-snes_converged_reason'

    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'

    end_time = 1.5
    dt = .1

# controls for linear iterations
  l_max_its = 100
  l_tol = 1e-8

# controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-8
[]

[Postprocessors]
    [stress_base]
        type = PointValue
        point = '0. 0. 0.'
        variable = stress_yy
    []
    [stress_exc_top]
        type = PointValue
        point = '0. 1. 0.'
        variable = stress_yy
    []
    [disp_top]
        type = PointValue
        point = '0. 2. 0.'
        variable = disp_y
    []
    [disp_mid]
        type = PointValue
        point = '0. 1. 0.'
        variable = disp_y
    []
    [react_y_base]
        type = PointValue
        point = '0. 0. 0.'
        variable = 'react_y'
    []
    [react_y_exc_top]
        type = PointValue
        point = '0. 1. 0.'
        variable = 'react_y'
    []
    [frozen_react_y_exc_top]
        type = PointValue
        point = '0. 1. 0.'
        variable = 'frozen_react_y'
    []
    [pwp_base]
        type = PointValue
        point = '0. 0. 0.'
        variable = porepressure
    []
    [pwp_exc_top]
        type = PointValue
        point = '0. 1. 0.'
        variable = porepressure
    []
[]

[Outputs]
    checkpoint = true
    [exodus]
    type = Exodus
    []
[]

@vabufano

[lindsayad]

Do you want this action to live in the MOOSE modules or in your own application? Are you looking for coding advice on how to write the Action? Your parameter desires all sound reasonable. There is some guidance on developing actions at https://mooseframework.inl.gov/source/actions/Action.html

[nlosacco]

Do you want this action to live in the MOOSE modules or in your own application?

Does it matter? I'd keep it in my own app for the time being but, is it a big deal?

Are you looking for coding advice on how to write the Action?

Yes! I would definitely need that. I've read the page you pointed at but it's still a bit obscure to me. If you have any other reference you can share...

Before I dive into this–quite scary—territory, I think there is one point I need to address: I think the most cumbersome bit of the approach, as implemented in the input file I shared earlier, is the stepwise procedure:

1. block displacements
2. calculate reactions
3. apply reactions
4. reduce reactions

Can an action be conceived in such a way that 1, 2 and 3 occur simultaneously at the beginning of a step (i.e. Dt = 0)?

[GiudGiud]

Can an action be conceived in such a way that 1, 2 and 3 occur simultaneously at the beginning of a step (i.e. Dt = 0)?

Are 1, 2 and 3 all performed by auxkernels and user objects?
If you list which ones we might be able to give you a way to do that from the input file using UO parameters for executioner ordering

Note that the creation of the objects is handled by the Action. the execution logic is coming from the objects' parameters. The action can set up the right execution logic when creating objects (the same a regular input does).

Creating an Action does not mean we can have a different execution logic than an input file would have been able to do

[nlosacco]

Hi @GiudGiud ,

for sake of completeness I'll add porewater pressure to the picture:

Are 1, 2 and 3 all performed by auxkernels and user objects?

1. is actually performed through custom BCs:

[BCs]
    [lock_exc_top_disp]
        type = LockBoundaryBC
         variable = disp_y
         locked_variable = disp_y
         boundary = 'top_exc'
         value = 0.0
         enable = false
    []
    [lock_exc_top_pwp]
        type = LockBoundaryBC
        variable = porepressure
        locked_variable = porepressure
        boundary = 'top_exc'
        value = 0.0
        enable = false

2. Computed nodal reactions and porepressures are stored into AuxVariables. I'm afraid the calculation of nodal reactions due to application of the previous BCs::lock_exc_top_disp needs a timestep to be performed, which I want to avoid. These are the blocks:

[AuxKernels]
  [freeze_react_y]
    type = FreezeForceAux
    variable = frozen_react_y
    source_variable = react_y
    freeze_time = 0.5
    execute_on = 'timestep_end'
  []
  [freeze_pwp]
    type = FreezeForceAux
    variable = frozen_porepressure
    source_variable = porepressure
    freeze_time = 0.5
    execute_on = 'timestep_end'
  []

3. Nodal reactions are then applied via a NodalKernels, porepressures through a BCs:

[NodalKernels]
    [force_exc_top]
        type = CoupledForceNodalKernel
        variable = disp_y
        v = relaxing_react_y
        boundary = 'top_exc'
        enable = false
    []
[]
[BCs]
    [exc_top_pwp]
        type = FunctorDirichletBC
        variable = porepressure
        boundary = 'top_exc'
        functor = 'relaxing_porepressure'
        enable = false
    []
[]

Creating an Action does not mean we can have a different execution logic than an input file would have been able to do

Look, honestly I don't even know if an Actions is the way to go, or if this should be pursued in a different way. This is actually what I'd like to hear from you

[GiudGiud]

Look, honestly I don't even know if an Actions is the way to go, or if this should be pursued in a different way. This is actually what I'd like to hear from you

ok let's figure it out in an input file first, and once it works it can be hardcoded into an Action

1. BCs
2. AuxKernels
3. BCs

I'm afraid the calculation of nodal reactions due to application of the previous BCs::lock_exc_top_disp needs a timestep to be performed,

Exactly.
But do you need to use time steps for time dependence (like in dynamics) ?
Otherwise it's fine to use time steps as load steps instead.