PorousFlowMassVolumetricExpansion is not working with Solid mechanics

[clariu25]

Hello, I'm trying to couple mechanics to a Hydraulic problem. The model has 3 layers and I'm trying to reach the steady state before I run a reservoir simulation. I ran the TH model without any problems, but when I added the M component, it did not converge anymore. Now I reduced the problem to HM.
Currently, if I remove the PorousFlowMassVolumetricExpansion kernel, the simulation runs but I don't understand why or if it makes sense. As far as I understand, this kernel lumps the mass component (that is in the quadrature point) into the nodes right? And that is necessary for the computation of the Time derivatives (mass in this case) as they are evaluated in the nodes, therefore, the mass calculations might not be correct.
Can somebody help me? Thank you very much

density_overburden = 2420

[Mesh]
  allow_renumbering = true
  [main]
    type = FileMeshGenerator
    file = mesh/3lyrs-PF-18km-Ldisks_sq1.msh
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  fp = the_simple_fluid
  displacements = 'disp_x disp_y disp_z'
[]

[Variables]
  [porepressure]
  []
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

[ICs]
  [Pressure]
    type = FunctionIC
    variable = porepressure
    function = 'hydrostaticP'
  []
[]

[BCs]
  [Pore_Pressure_BC]
    type = FunctionDirichletBC
    variable = porepressure
    boundary = 'Top Bottom Left Back Right Front'
    function = 'hydrostaticP'
  []
  [Disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'Back Front'
    value = 0
  []
  [Disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'Left Right'
    value = 0
  []
  [Disp_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'Bottom'
    value = 0
  []
  [litho_stress_zz]
    type = Pressure
    boundary = 'Top'
    variable = disp_z
    function = 'lithost_zz_BC' #Lithostatic pressure
  []
[]

[Functions]
  [hydrostaticP]
    type = ParsedFunction
    expression = '1000*9.8*(-z)+1e5'
  []
  [dts]
    type = PiecewiseLinear
    x = '0 3600 86400 864000 15768000
          15768001 15854400 16632000 31536000'
    y = '1 1800 43200 172800 864000
         3600 43200 86400 864000'
  []
  [lithost_zz_BC]
    type = ParsedFunction
    expression = '${density_overburden}*9.81*(-z)'
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    viscosity = 0.00026383
    thermal_expansion = 0
  []
[]

[Kernels]
  ##### Hydraulics #####
  [mass_dot]
    type = PorousFlowFullySaturatedMassTimeDerivative
    variable = porepressure
  []
  [advection]
    type = PorousFlowFullySaturatedAdvectiveFlux
    fluid_component = 0
    variable = porepressure
    gravity = '0 0 -9.81'
  []
  ####### Mechanics ########
  [gravity]
    type = Gravity
    use_displaced_mesh = false
    variable = disp_z
    value = -9.81 # MPa
    density = ${density_overburden}
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    variable = porepressure
  []
  [grad_stress_x]
    type = StressDivergenceTensors
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  []
  [grad_stress_z]
    type = StressDivergenceTensors
    variable = disp_z
    use_displaced_mesh = false
    component = 2
  []
  [poro_z]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_z
    use_displaced_mesh = false
    component = 2
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure disp_x disp_y disp_z'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [fluid_mass_prod]
    type = PorousFlowSumQuantity
  []
[]

[Materials]
  [undrained_density]
    type = PorousFlowTotalGravitationalDensityFullySaturatedFromPorosity
    rho_s = 3000 ## (por=15%, pf=1000, pb=2700)
    output_properties = density
    outputs = 'exodus'
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 1
    fluid_bulk_modulus = 2e9 #Keep it consistant with the FluidProperties
  []
  [saturation_calculator]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [water]
    type = PorousFlowSingleComponentFluid
    phase = 0
  []
  [porosity_const]
    type = PorousFlowPorosityConst
    porosity = 0.15
  []
  [permeability_res]
    type = PorousFlowPermeabilityConst
    permeability = '2.5e-14 0 0  0 2.5e-14 0  0 0 2.5e-14'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [eff_fluid_pressure]
    # Used by variant porosity
    type = PorousFlowEffectiveFluidPressure
  []
  [temperature]
    type = PorousFlowTemperature
    # temperature = temperature
  []
  ######## Mechanics ########
  [elasticity_tensor_res]
    type = ComputeIsotropicElasticityTensor
    bulk_modulus = 55e9
    shear_modulus = 32.59e9
  []
  [strain]
    type = ComputeSmallStrain
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Preconditioning]
  active = 'HM'
  [smp]
    type = SMP
    full = true
  []
  [superlu]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'gmres lu superlu_dist'
  []
  [HM]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_agg_nl  -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_max_levels  -pc_hypre_boomeramg_coarsen_type   -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_P_max -pc_hypre_boomeramg_truncfactor'
    petsc_options_value = 'hypre    boomeramg 0.7 4 5 25 HMIS  ext+i 2 0.3'
  []
[]

[Executioner]
  type = Transient
  end_time = 31536000 #1 year (365 days)
  [TimeStepper]
    type = FunctionDT
    function = dts
  []
  dtmin = 1
  nl_rel_tol = 1e-8
  nl_rel_step_tol = 1e-8
  l_tol = 1e-8
  automatic_scaling = true
  compute_scaling_once = true
  solve_type = NEWTON
  start_time = 0
  steady_state_detection = true
  steady_state_tolerance = 1e-09
[]

[Outputs]
  console = true
  csv = false
  exodus = true
  print_linear_residuals = false
[]

[Debug]
  show_var_residual_norms = true
[] 

[GiudGiud]

Hello

I ran the TH model without any problems

did you have the volumetric_expansion set to 0 then?

how big is the mesh?

There are instructions on troubleshooting convergence issues there
https://mooseframework.inl.gov//application_usage/failed_solves.html

[clariu25]

I did not have this kernel then. I added the expansion one only when coupling with mechanics

The mesh has 586.998 elements and covers 18km x 18km x 1km. It has very different sizes as my interest area is only 10m thick but over a huge area. Therefore, I need a lot of detail there but not so much in the other layers. I should have enough elements in the transition between fine and coarse. It is definitely not an ideal mesh so I will test again with a different mesh

Thanks for the link!

[GiudGiud]

But you have volumetric expansion set to 0 in this version of the input file?

[GiudGiud]

Does it solve if you are solving both TH and mechanics, but with 0 expansion so they are effectively decoupled?

[clariu25]

No, it still had the same convergence issue.
So far if I remove the Top BC for porepressure it converges. But if I change the fluid properties, I get the same problem so it is not really fixed

[josebastiase]

Hi,

I dont see any issues with the kernels... Have you tried a simple mesh? I think maybe the issue is with your BCs. For overburden, I typically use NeumannBC.

[clariu25]

Not really.. As the Hydraulic simulation worked on this mesh and the issue is with pressure I did not think about changing the mesh. Maybe I should, Guillaume also asked about the mesh. Thanks
Regarding the BC, I checked with both NeumannBC and Pressure and the results were the same

[clariu25]

Hello,
I did some tests with the BC, and so far removing the hydrostatic pressure at the Top boundary allows the simulation to converge. The results look fine, the pressure still reaches the expected values. But this solution only applies in the case where I use SimpleFluidProperties. For something more complex like Water97Properties it doesn't. Meaning that the problem is not solved yet

[Pore_Pressure_BC]
    type = FunctionDirichletBC
    variable = porepressure
    ## boundary = 'Top Bottom Left Back Right Front'
    boundary = 'Bottom Left Back Right Front'
    function = 'hydrostaticP'
  []

This is the current file: https://github.com/clariu25/PF_tests/blob/main/HM_5km_SF.i
And the mesh: https://github.com/clariu25/PF_tests/blob/main/5km_mesh.msh

[GiudGiud]

Seems odd to me to impose pressure boundary conditions everywhere?

Usually in flow simulations, you set a pressure BC somewhere (even a single node!) and let the equations dictate the pressure elsewhere.
Sometimes you set pressure in two locations, and the pressure gradient creates the flow. But pressure everywhere?

[clariu25]

Ooook. I always set them fully as, in the underground, the pressure should remain constant at the borders if I'm injecting in the center of my domain. But indeed, by applying the IC with a gradient, it will reach the actual pressure value, due to the gradient and the stress effect. Maybe I had a problem of concept. Thanks and sorry for my dumb mistake.

[cpgr]

It's kind of hard to know what is the culprit. Can you make a simple generated mesh that is similar so that we can try and run it?

[clariu25]

Hello, sure. I was trying to solve it on my own. Sorry, thanks for the help. Here I attached both the mesh and the input file
HM_problem.zip

I tested with this mesh but the convergence behavior was very similar as before