Single-phase multi-component boundary conditions and kernel settings

[beidongW]

I am implementing single-phase multi-component gas flow and found that my .i file fails to converge if I do not set boundary conditions. Components only stop accumulating at the outlet when I apply DirichletBC to fix the mass fractions. I would like to request assistance with three questions regarding my input file:
1、How should I set the boundary conditions to achieve free flow?
2、Is my kernel usage correct?
3、I have implemented a custom startup pressure gradient kernel that public PorousFlowFullySaturatedDarcyBase named StartFlowPressureDarcyFlow. Can this be used here? I still do not fully understand the relationship between PorousFlowAdvectiveFlux and Darcy flow.

# ----------------------
# name：3component.i
# SINGLE-PHASE MULTI-COMPONENT POROUSFLOW EXAMPLE
# Darcy flow with 3 gas components (CO2, CH4, N2)
# ----------------------

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 30
  ny = 30
  xmin = 0
  xmax = 5
  ymin = 0
  ymax = 5
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]


[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'p massfrac0 massfrac1'
    number_fluid_phases = 1
    number_fluid_components = 3
  []
[]


[Variables]
  [p]
    order = FIRST
    family = LAGRANGE
  []

  [massfrac0]
    order = FIRST
    family = LAGRANGE
  []

  [massfrac1]
    order = FIRST
    family = LAGRANGE
  []
[]


[AuxVariables]
  [viscosity_field]
    family = MONOMIAL
    order = CONSTANT
  []
  [massfrac2]
    family = MONOMIAL
    order = CONSTANT
  []
  [./porosity_field]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [output_viscosity]
    type = PorousFlowPropertyAux
    variable = viscosity_field
    property = viscosity
  []
  [massfrac2_aux]
    type = ParsedAux
    variable = massfrac2
    coupled_variables = 'massfrac0 massfrac1'
    expression = 'max(0.0, min(1.0, 1.0 - massfrac0 - massfrac1))'
  []
  [./porosity_field_aux]
    type = ParsedAux
    variable = porosity_field
    expression = '0.15 + 0.05*sin(2*3.14*x/5)'
    use_xyzt = true
    execute_on = INITIAL
  [../]

[]

[ICs]
  [p_init]
    type = FunctionIC
    variable = p
    function = '(5-x)*2e6+1e7'
  []

  [mf0_init]
    # type = ConstantIC
    # value =0.3
    type = RandomIC
    variable = massfrac0
    min = 0
    max =0.3
  []

  [mf1_init]
    # type = ConstantIC
    # value = 0.4
    type = RandomIC
    variable = massfrac1
    min = 0
    max = 0.4
  []
[]


[BCs]
  [p_right]
    type = DirichletBC
    variable = p
    value = 1e7
    boundary = right
  []
  [mf0_left]   
    type =  DirichletBC
    variable = massfrac0
    value = 0.0
    boundary = right
  []

  [mf1_left]   
    type =  DirichletBC
    variable = massfrac1
    value = 0.0
    boundary = right
  []

[]


[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [ch4]
    type = MethaneFluidProperties
  []
  [n2]
    type = NitrogenFluidProperties
  []
[]


[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 300
  []

  [phase]
    type = PorousFlow1PhaseFullySaturated
    porepressure = p
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'massfrac0 massfrac1'
    #  at_nodes = true
  []
  [fluid]
    type = SimpleMixtureFluidProperties
    PorousFlowDictator = dictator
    phase = 0
    fp0 = n2
    fp1 = ch4
    fp2 = co2
  []
  [diff]
    type = PorousFlowDiffusivityConst
    diffusion_coeff = '1e-5 1e-5 1e-5'
    tortuosity = 0.1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = porosity_field
  []

  [relperm]
    type = PorousFlowRelativePermeabilityConst
    PorousFlowDictator = dictator
    phase = 0
    kr = 1.0
  []

  [perm]
    type = PorousFlowPermeabilityConst
    permeability = '1e-14 0 0  0 1e-14 0  0 0 1e-14'
  []
[]

[Kernels]
  [./darcy_velocity]
    type = StartFlowPressureDarcyFlow
    variable = p
    gravity = '0 9.8 0'
    startup_grad = '1e4 1e4 1e4'      
    transition_width = '5e2 5e2 5e2'  
  [../] 
  [mass_p]
    type = PorousFlowMassTimeDerivative
    variable = p
    fluid_component = 2
  []
  [adv_p]
    type = PorousFlowAdvectiveFlux
    variable = p
    gravity = '0 9.8 0'
    fluid_component = 2
  []

  [mass_mf0]
    type = PorousFlowMassTimeDerivative
    variable = massfrac0
    fluid_component = 0
  []

  [adv_mf0]
    type = PorousFlowAdvectiveFlux
    variable = massfrac0
    fluid_component = 0
    gravity = '0 9.8 0'
  []

  [mass_mf1]
    type = PorousFlowMassTimeDerivative
    variable = massfrac1
    fluid_component = 1
  []

  [adv_mf1]
    type = PorousFlowAdvectiveFlux
    variable = massfrac1
    fluid_component = 1
    gravity = '0 9.8 0'
  []

  [diff_mf0]
    type = PorousFlowDispersiveFlux
    variable = massfrac0
    disp_long = 0
    disp_trans = 0
    gravity = '0 9.8 0'
    fluid_component = 0    
  []
  [diff_mf1]
    type = PorousFlowDispersiveFlux
    variable = massfrac1
    disp_long = 0
    disp_trans = 0
    gravity = '0 9.8 0'
    fluid_component = 0  
  []
[]


[Executioner]
  type = Transient
  solve_type = Newton
  start_time = 0
  end_time = 4e4
  dtmax =1000
  dt = 1
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-5

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 20
    optimal_iterations = 7       
    iteration_window = 3
    growth_factor = 1.5
    cutback_factor = 0.5
   [../] 
[]

[Preconditioning]
  [petsc_lu]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu mumps'
  []
[]


[Outputs]
  [exodus_out]
    type = Exodus
    execute_on = 'initial timestep_end'
    file_base = problem/test1/personalmodel/testdata/testdata/3component_1
  []
[]

[GiudGiud]

@beidongW Is this other post relevant? As in you are adding the additional component flow / scalar species advection to this free flow problem?
#30311

[beidongW]

My target to simulate gas reservoir depletion development. The physical meaning of the model is that a Dirichlet constant pressure condition is applied at the outlet, and components are allowed to flow out freely. The current issue is that if I use Dirichlet constant pressure at the outlet, components will accumulate at the outlet boundary, with values exceeding 10, as shown in the figure. However, I can run the official case (porous_flow/test/tests/sinks/injection_production_eg_outflowBC.i) successfully with PorousFlowOutflowBC, but it doesn't work in my own three-component model. So now I suspect there might be an issue with my kernel implementation, as I don't fully understand the operating principles of PorousFlowMassTimeDerivative and PorousFlowAdvectiveFlux. I wonder if using component mass fractions as variables and calculating them with these two kernels is fundamentally wrong.

phase = 0 is liquid phase
phase = 1 is gas phase
fluid_component = 0 is water
fluid_component = 1 is CO2

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 30
  ny = 30
  xmin = 0
  xmax = 5
  ymin = 0
  ymax = 5
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [frac_water_in_liquid]
  []
  [frac_water_in_gas]
  []
  [water_kg_per_s]
  []
  [co2_kg_per_s]
  []
[]

[ICs]
  [pwater_init]
    type = FunctionIC
    variable = pwater
    function = '(5-x)*2e6+1e7'
  []
  [pgas_init]
    type = FunctionIC
    variable = pgas
    function = '(5-x)*2e6+1.1e7'
  []
  [./saturation_gas]
    type = RandomIC
    variable = saturation_gas
    min= 0
    max= 0.15
  [../] 
  [./frac_water_in_liquid]
    type = ConstantIC
    variable = frac_water_in_liquid
    value = 1
  [../]
  [./frac_water_in_gas]
    type = ConstantIC
    variable = frac_water_in_gas
    value = 0     
  [../]
[]  

[BCs]
  [pwater_right]
    type = DirichletBC
    boundary = right
    value = 1e7
    variable = pwater
  []
  [pgas_right]
    type = DirichletBC
    boundary = right
    value = 1.1e7
    variable = pgas
  []
  [right_water_component0]
    type = PorousFlowOutflowBC
    boundary = right
    variable = pwater
    mass_fraction_component = 0
    save_in = water_kg_per_s
  []
  [right_co2_component1]
    type = PorousFlowOutflowBC
    boundary = right
    variable = pgas
    mass_fraction_component = 1
    save_in = co2_kg_per_s
  []
[]
[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
[]

[Variables]
  [pwater]
  []
  [pgas]
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = pgas
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas pwater'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedBicubicFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 293.15
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.2
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
[]



[Preconditioning]
  active = 'basic'
  [basic]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu NONZERO 2'
  []
  [preferred]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  start_time = 0
  end_time = 50000
  dtmax =500
  dt = 1
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-5
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 20
    optimal_iterations = 7       # Newton迭代次数理想值 5~10
    iteration_window = 3
    growth_factor = 1.5
    cutback_factor = 0.5
   [../] 
[]


[Outputs]
  [exodus_out]
    type = Exodus
    execute_on = 'initial timestep_end'
    file_base = problem/test1/personalmodel/testdata/testdata/ceshi2_1
  []
[]

[GiudGiud]

non-expert here, what does it mean to have this

  [pwater_right]
    type = DirichletBC
    boundary = right
    value = 1e7
    variable = pwater
  []
  [pgas_right]
    type = DirichletBC
    boundary = right
    value = 1.1e7
    variable = pgas
  []

if you specify different value, do you effectively dictate the concentration of each phase?

[WilkAndy]

Yes, @GiudGiud , that is correct, so you're not a complete non-expert !! :-). Usually we don't use the term "concentration" (because that's reserved for solutes) but instead use "saturation", since we're talking about water + something else (in this case, gas). There is a thing called capillary pressure, which is pgas - pwater, (so is 0.1e7 at the boundary) and in this case is a PorousFlowCapillaryPressureVG UserObject. Given a value of pgas - pwater, the capillary pressure function gives the saturation.

[WilkAndy]

Hi @beidongW ,

Before i write my suggestions, I just want to reiterate that BCs are very difficult to understand in multiphase, multicomponent models (whether they be in MOOSE or other codes) so please don't feel sad if you struggle to get them correct. Even experts frequently make mistakes. When you've eventually got them correct, you'll feel very happy !

1. I suggest starting with a simple single-component, single-phase example, using SimpleFluidProperties, and only PorousFlow objects (not your modifications).
2. Then do 2-components, single-phase, still using SimpleFluidProperties and PorousFlow objects
3. Then do multiphase, still using SimpleFluidProperties and PorousFlow objects
4. Then do your real case.

These steps might take some days of experimentation. You can find (1)-(3) in the PorousFlow examples and tests. But you'll learn a lot by trying stuff yourself.

In my opinion, you've taken the wrong approach with your BCs. Even though it might work in this one situation, it's just luck. Please look again at https://mooseframework.inl.gov/modules/porous_flow/boundaries.html , eg, there's a few paragraphs on why DirichletBCs are usually bad to use, and you should use PorousFlowSink-type instead. In the future, and maybe here too, you're very likely to want https://mooseframework.inl.gov/modules/porous_flow/sinks.html as well, since i think you'll be simulating wells.

  [pwater_right]
    type = DirichletBC
    boundary = right
    value = 1e7
    variable = pwater
  []

This is saying: "fix the variable pwater to 1e7 by providing/removing fluid component 0". In this case fluid component 0 is water (fluid component 0 is not always water: that depends on the rest of your input file). Similarly, for the other DirichletBC. What if MOOSE gets to a situation with high pressure, BUT, there's no fluid component 0 at the node - it'll try to remove a non-existent fluid component, and will either crash or not converge.

  [right_water_component0]
    type = PorousFlowOutflowBC
    boundary = right
    variable = pwater
    mass_fraction_component = 0
    save_in = water_kg_per_s
  []

This is saying: "to the right of the right boundary, there is another domain that i'm choosing not to model. That domain has permeability and porosity similar to my model. It can act as a source of water (if the porepressure gradient in my model is negative) or a sink of water (if the porepressure gradient in my model is positive), and the speed of water flow into/out-of my model will be dictated by my model's permeability." I think this isn't what you want, conceptually (but i could be misunderstanding you, sorry): use a PorousFlowSink , or PorousFlowPolyLineSink or Peaceman borehole, as at the aforementioned webpages.

[beidongW]

Thank you very much. I will gradually follow the suggestions to experiment.