Two-phase porous media flow: Strange pressure field

[zhang-rite]

Hi MOOSE Porous media team

I am now testing a simple 2D CO2 injection problem (coupled with mechanics and with gravity effect on) based on the "modules/porous_flow/examples/thm_example/2D.i" example.

The model simply has an aquifer with a high permeability (1e-12 m2) between upper caprock and bottom both with low permeability (1e-19 m2).

As expected, while CO2 injection, under bouncy forces, the CO2 plume gradually moves to upper regions on in the aquifer. However, the strange thing is in the change of fluid pressure field, which is much large in the caprock than expected.

Could you help me point out the problem? Thank you.

The input file link

(The pressure change (to the hydrostatic pressure) is denoted "pgas2". )

[WilkAndy]

No answer, just some questions

1. why don't you have a PorousFlowEffectiveStressCoupling Kernel for disp_z ?
2. Have you checked your sidesets in your mesh? Viz, the 'left' and 'injection_area' sidesets?
3. why are your bulk moduli so high?
4. Have you checked your initial conditions?
5. Are the displacements as expected?

a

[zhang-rite]

Andy, thank you very much for quick reply. Following your questions I have revised the input file and done modeling again.

why don't you have a PorousFlowEffectiveStressCoupling Kernel for disp_z ?

• Now I have added, however no apparent effect for above mentioned pressure field.

Have you checked your sidesets in your mesh? Viz, the 'left' and 'injection_area' sidesets?

• It seems the 'injection area' functions. This figure shows that the pressure field during initial pressure build up (no CO2 enter into the aquifer domain). We see the pressure build up well in vertical range of the aquifer.
  The above mentioned pressure anomaly appeared after CO2 plume migrated to upper aquifer.

why are your bulk moduli so high?

• Oh, I Just followed the example "modules/porous_flow/examples/thm_example/2D.i" without checking 😂. Now I have changed they to low values; however, the pressure problem remains.

Have you checked your initial conditions?

• I specially use another run to get an hydrostatic initial conditions. Link here

The problem. Here I do not know how to inherit the initialized strain field to the following run.

• In the initial run

  [./strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'ini_stress'
  [../]
  [./ini_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '-12.8E6 0 0  0 -12.8E6 0  0 0 -12.8E6'
    eigenstrain_name = ini_stress
  [../]

• In the subsequent CO2 injection run, I commented out the "ini_strain" block.

  [./strain]
    type = ComputeAxisymmetricRZSmallStrain
    # eigenstrain_names = 'ini_stress'
  [../]
  # [./ini_strain]
  #   type = ComputeEigenstrainFromInitialStress
  #   # initial_stress = '-12.8E6 0 0  0 -12.8E6 0  0 0 -12.8E6'
  #   # eigenstrain_name = ini_stress
  # [../]

Now I have changed the setting in the subsequent CO2 injection run the same as the initial run. The result seems better.
Could tell me what is the right way for inheriting strain stress field from the initial run?

Are the displacements as expected?

• The strain field has the same problem as the pressure field.

[zhang-rite]

Else is there a way to read "ini_strain" from output file of an initialization run? I want to get more realistic vertical stress by an initialization run.
Like SolutionFunction? initial_from_file_var?

[zhang-rite]

In https://mooseframework.inl.gov/source/materials/ComputeEigenstrainFromInitialStress.html, there writes "Using AuxVariables to set the initial stress may be advantageous when reading from a solution file using a SolutionAux."
So the idea is through "SolutionUserObject" -> "AuxVariables" -> "initial_stress_aux"?
I wish there is an example using "initial_stress_aux" with restarting run.

[WilkAndy]

You're correct. The initial run has to save its stress field into AuxVariables, then, in the second run, you have to use a SolutionUserObject to read those into initial_stress_aux. I thought i had included an example of that in the MOOSE repo as i do it all the time in my work, but it's not there in full. You can combine the ideas from the following, however:

modules/tensor_mechanics/test/tests/initial_stress/gravity_with_aux.i
test/tests/auxkernels/solution_aux/solution_aux_exodus.i

[zhang-rite]

Andy, thank you again for the help.

I find, as your suggested way, after load and initialize strain using the output of a previous run (already reached steady state), there are still gradual changes in displacement (and strain and stress) always ( since run even after a very long time run) while without fluid injection.
I do keep all the boundary conditions and mesh the same as the initialization run.
I just want to the get gravity equilibrium fields of stress-strain and pressure before injection.

Could you help me point out the problem. The file link.

Best,

Yi Zhang

[WilkAndy]

There are two possibilities here:

1. You've made a mistake
2. You have not made a mistake, and the result is expected, and annoying, and i'm not sure how to eliminate the problem.

To determine whether (1) or (2) is the case, look at your stresses and displacements - are they changing by very much (without fluid injection)?

The result is expected because when writing/reading exodus files you are averaging stress over an element. However, internally within MOOSE, stress is represented at the quadpoints of the element (there are 4 in each element in your case). For the initial configuration, MOOSE finds the stresses at those quadpoints, then averages them to write to the exodus file. Then you read those averaged versions and assign those averaged versions to each quadpoint. Then MOOSE sees the result is not in exact equilibrium, and changes the stresses (and displacements) slightly to achieve equilibrium (so you're back at the true equilibrium situation MOOSE computed in the initial simulation).

What i've done in cases like this is to create some "change_disp" AuxVariables that record the change in displacement after the first timestep. But that's a bit of a hack. Hopefully those disps and stress changes are so small you can ignore them.

[zhang-rite]

Andy, thanks!
Currently I may get the differential change to first timestep by using ParsedAux in AuxKernels with execute_on = initial and doing subtraction in another ParsedAux. But how could I check the differential changes in displacement after first timestep (relative to the step 2)?
Is there something like execute_on = step 2?

[WilkAndy]

You could use Controls. Turn the AuxKernel off after a time. Controls are a really useful feature of MOOSE.

[WilkAndy]

Here's a better example: https://mooseframework.inl.gov/source/controls/TimePeriod.html

[zhang-rite]

By viewing using Paraview, I find the gradual changes in displacement (or strain, pressure) occur at the transition between the high-permeability aquifer zone and low-permeability zones. The reason is still unknow. Perhaps the error comes from the physics setting some where. I have not found the improper setting.

I will try Controls system. Thanks!

[zhang-rite]

I had tested another whole day today, but still cannot work out the problem. Is there a two-phase flow mechanical coupled problem with first gravity equilibrium? (I wish). Thank you.

[zhang-rite]

Is there any recommended strategy for setting an initial gravity equilibrium for hydromechanical problems?

[WilkAndy]

Sorry for not answering sooner. I'm pretty sure you haven't made a mistake, @zhang-rite , and i'm pretty sure that your problem isn't related to 2-phase flow. (You could try a single-phase version just to check - there are examples in the tests of porous flow. If the problem is 2-phase flow related it's probably a bug in your input file, but i can't see one.) Alternatives to what you're doing already are recover/restart (https://mooseframework.inl.gov/application_usage/restart_recover.html) or: (a) read from your steady-state; (b) do a very long-time transient without injecting to get the quadpoint stresses almost-perfect; (c) then start injecting.

[zhang-rite]

Andy, thank you very much for your help!
Yes, the initialization can be viewed as a 1-phase problem. Now I am trying to do as the suggested way--starting from 1-phase steady state run. However, I again encountered problem for the convergence, which I have posted in another discussion.
#17091
I will also try the "recover/restart" way.