Unexpected high pressure point and convergence issue

[JuliaHLee]

Check these boxes if you have followed the posting rules.

• Q&A Navier Stokes is the most appropriate category 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 (avoid screenshots if possible, triple back quotes before/after pasted text, etc)

Question

Dear MOOSE experts,

I have been developing a 3D thermo-hydraulic model with Navier-Stokes module using incompressible finite volume. I run a transient simulation. The problems are that the model shows unexpected pressure development and slowing down with convergence. I have a cube filled with spheres (solid phase). Inside the cube, water flows in z direction from the top to the bottom through the spheres, while the side is set with no flow boundary. Please see the model showing the mesh below:

The issues are:

1. I have noticed a problem, that there is unexpected high pressure at one point near a side boundary surface, while I expected a pressure gradient in z direction. As the geometry has the symmetric structure, it should not happen. Could it be a problem with the mesh, choice of solver, time-step or convergence criteria?

2. After around 200 time steps, the model does not converge with the given dt = 1 and continues with a reduced time step dt = 0.5. And it repeats that, slowing down the simulation. Can it be related to the first issue shown by the high pressure point?

3. When I reduce time step to dt = 0.01, hoping for a better convergence, the high pressure at one point even grows over time. And the pressure value gets higher than the pressure gradient range. Is it showing that the model is not stable?

I am not sure about the cause of the problems. So it would be great if I could get some tips on the issue.

Many thanks in advance!

Best regards,

Julia

Additional information

Mesh size and type: Num DOFs: 370877
Reynolds number: Re < 2.5
Discretization (finite element CG/DG, finite volume, etc): Finite volume
Models (turbulence, porous media, etc): porous media
Solver method (fully coupled, segregated, multiapps, etc): Newton
Base input you started from:

Details

Input_file

T_0 = 283.15
T_in = 303.15
p_outlet = 0
lambda_sol = 1.5
cs = 840
rho_s = 2558
lambda_fl = 0.59
cf = 4180
rho_f = 1000
mu = 0.001002
vf = 5e-05
hsf = 1e+08


[Mesh]
  [mesh_file]
    type = FileMeshGenerator
    file = mesh_spheres_in_box.msh
  []
  [sideset]
    input = mesh_file
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 'box'
    paired_block = 'solid'
    new_boundary = sideset
  []
[]


[GlobalParams]
  two_term_boundary_expansion = false
  rhie_chow_user_object = rc
[]
[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    w = vel_z
    block = 'box'
    pressure = pressure
  []
[]
[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 0
  []
  [vel_y]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 0
  []
  [vel_z]
    type = INSFVVelocityVariable
    block = box
    initial_condition = ${fparse -1*vf}
  []
  [pressure]
    type = INSFVPressureVariable
    block = 'box'
  []
  [T]
    type = INSFVEnergyVariable
    block = 'box'
    initial_condition = ${T_0}
  []
  [Ts]
    type = INSFVEnergyVariable
    initial_condition = ${T_0}
    block = 'solid'
  []
[]
[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    rho = ${rho_f}
  []
  [u_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
    block = 'box'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = x
    block = 'box'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = x
    pressure = pressure
    block = 'box'
  []
  [v_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
    block = 'box'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = y
    block = 'box'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = y
    pressure = pressure
    block = 'box'
  []
  [w_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
  []
  [w_advection]
    type = INSFVMomentumAdvection
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
    block = box
  []
  [w_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_z
    mu = ${mu}
    momentum_component = z
    block = box
  []
  [w_pressure]
    type = INSFVMomentumPressure
    variable = vel_z
    momentum_component = z
    pressure = pressure
    block = box
  []
  [energy_time]
    type = INSFVEnergyTimeDerivative
    variable = T
    rho = ${rho_f}
    dh_dt = dh_dt
    block = 'box'
  []
  [temp_advection]
    type = INSFVEnergyAdvection
    variable = T
    block = 'box'
  []
  [temp_conduction]
    type = FVDiffusion
    coeff = ${lambda_fl}
    variable = T
    block = 'box'
  []
  [grain_energy_time]
    type = INSFVEnergyTimeDerivative
    variable = Ts
    rho = ${rho_s}
    dh_dt = dh_solid_dt
    block = 'solid'
  []
  [grain_conduction]
    type = FVDiffusion
    coeff = ${lambda_sol}
    variable = Ts
    block = 'solid'
  []
[]
[FVInterfaceKernels]
  [convection]
    type = FVConvectionCorrelationInterface
    variable1 = T
    variable2 = Ts
    boundary = sideset
    h = ${hsf}
    T_solid = Ts
    T_fluid = T
    subdomain1 = 'box'
    subdomain2 = 'solid'
    wall_cell_is_bulk = true
  []
[]
[FVBCs]
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_x
    functor = 0
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_y
    functor = 0
  []
  [inlet-w]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_z
    functor = ${fparse -1*vf}
  []
  [inlet_T]
   type = FVDirichletBC
    variable = T
    boundary = top
    value = ${T_in}
  []
  [no-slip-u]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_x
    function = 0
  []
  [no-slip-v]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_y
    function = 0
  []
  [no-slip-w]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_z
    function = 0
  []
  [outlet_p]
    type = INSFVOutletPressureBC
    boundary = bottom
    variable = pressure
    function = ${p_outlet}
  []
[]

[FunctorMaterials]
  [fluid_functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = cp
    prop_values = ${cf}
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = T
    rho = ${rho_f}
    block = 'box'
  []
  [solid_functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = cp_s
    prop_values = ${cs}
  []
  [solid_ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = Ts
    cp = cp_s
    rho = ${rho_s}
    h = h_solid
    rho_h = rho_h_solid
    block = 'solid'
  []
[]
[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_pc_type -sub_pc_factor_shift_type'
  petsc_options_value = 'asm      100                lu           NONZERO'
  line_search = 'none'
  l_tol = 1e-8
  nl_rel_tol = 1e-8
  dt = 1
  end_time = 1200
[]
[Outputs]
  exodus = true
  csv = true
[]

[GiudGiud]

Hello

Do you have skewness-correction turned on? Your mesh is made of tets right?

[GiudGiud]

hello

If you use a GeneralFunctorFluidProps you dont need to define cp as a functor like this

  [fluid_functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = cp
    prop_values = ${cf}
  []

So it's either you define all the properties (including mu) using separate materials,
or you use the GeneralFunctorFluidProps with a fluid properties object

[JuliaHLee]

Thank you for the detailed explanation. But it still does not converge showing increasing Nonlinear residual from the first time step:

Time Step 1, time = 0.01, dt = 0.01
 0 Nonlinear |R| = 3.973301e+01
      0 Linear |R| = 3.973301e+01
      1 Linear |R| = 4.118636e-10
 1 Nonlinear |R| = 6.406181e+01
      0 Linear |R| = 6.406181e+01
      1 Linear |R| = 1.707751e-10
 2 Nonlinear |R| = 4.559898e+02
      0 Linear |R| = 4.559898e+02
      1 Linear |R| = 1.475171e-11
 3 Nonlinear |R| = 2.769738e+03
      0 Linear |R| = 2.769738e+03
      1 Linear |R| = 6.053305e-09
 4 Nonlinear |R| = 8.219541e+04

Updated_inputfile

T_0 = 283.15
T_in = 303.15
p_outlet = 0
lambda_sol = 1.5
cs = 840
rho_s = 2558
lambda_fl = 0.59
cf = 4180
rho_f = 1000
mu = 0.001002
vf = 5e-05
hsf = 1e+08


[Mesh]
  [mesh_file]
    type = FileMeshGenerator
    file = mesh_spheres_in_box.msh
  []
  [sideset]
    input = mesh_file
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 'box'
    paired_block = 'solid'
    new_boundary = sideset
  []
[]


[GlobalParams]
  rhie_chow_user_object = rc
[]
[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    w = vel_z
    block = 'box'
    pressure = pressure
  []
[]
[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 1e-10
    face_interp_method = 'skewness-corrected'
  []
  [vel_y]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 1e-10
    face_interp_method = 'skewness-corrected'
  []
  [vel_z]
    type = INSFVVelocityVariable
    block = box
    initial_condition = ${fparse -1*vf}
    face_interp_method = 'skewness-corrected'
  []
  [pressure]
    type = INSFVPressureVariable
    block = 'box'
    face_interp_method = 'skewness-corrected'
  []
  [T]
    type = INSFVEnergyVariable
    block = 'box'
    initial_condition = ${T_0}
    face_interp_method = 'skewness-corrected'
  []
  [Ts]
    type = INSFVEnergyVariable
    initial_condition = ${T_0}
    block = 'solid'
    face_interp_method = 'skewness-corrected'
  []
[]

[AuxVariables]
  [velocity_norm]
    type = MooseVariableFVReal
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    rho = ${rho_f}
    advected_interp_method = 'skewness-corrected'
  []
  [u_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
    block = 'box'
    advected_interp_method = 'skewness-corrected'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = x
    block = 'box'
    variable_interp_method = 'skewness-corrected'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = x
    pressure = pressure
    block = 'box'
    correct_skewness = 'true'
  []
  [v_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
    block = 'box'
    advected_interp_method = 'skewness-corrected'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = y
    block = 'box'
    variable_interp_method = 'skewness-corrected'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = y
    pressure = pressure
    block = 'box'
    correct_skewness = 'true'
  []
  [w_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
  []
  [w_advection]
    type = INSFVMomentumAdvection
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
    block = box
    advected_interp_method = 'skewness-corrected'
  []
  [w_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_z
    mu = ${mu}
    momentum_component = z
    block = box
    variable_interp_method = 'skewness-corrected'
  []
  [w_pressure]
    type = INSFVMomentumPressure
    variable = vel_z
    momentum_component = z
    pressure = pressure
    block = box
    correct_skewness = 'true'
  []
  [energy_time]
    type = INSFVEnergyTimeDerivative
    variable = T
    rho = ${rho_f}
    dh_dt = dh_dt
    block = 'box'
  []
  [temp_advection]
    type = INSFVEnergyAdvection
    variable = T
    block = 'box'
    advected_interp_method = 'skewness-corrected'
  []
  [temp_conduction]
    type = FVDiffusion
    coeff = ${lambda_fl}
    variable = T
    block = 'box'
    variable_interp_method = 'skewness-corrected'
  []
  [grain_energy_time]
    type = INSFVEnergyTimeDerivative
    variable = Ts
    rho = ${rho_s}
    dh_dt = dh_solid_dt
    block = 'solid'
  []
  [grain_conduction]
    type = FVDiffusion
    coeff = ${lambda_sol}
    variable = Ts
    block = 'solid'
    variable_interp_method = 'skewness-corrected'
  []
[]
[FVInterfaceKernels]
  [convection]
    type = FVConvectionCorrelationInterface
    variable1 = T
    variable2 = Ts
    boundary = sideset
    h = ${hsf}
    T_solid = Ts
    T_fluid = T
    subdomain1 = 'box'
    subdomain2 = 'solid'
    wall_cell_is_bulk = true
  []
[]
[FVBCs]
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_x
    functor = 1e-10
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_y
    functor = 1e-10
  []
  [inlet-w]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_z
    functor = ${fparse -1*vf}
  []
  [inlet_T]
   type = FVDirichletBC
    variable = T
    boundary = top
    value = ${T_in}
  []
  [no-slip-u]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_x
    function = 0
  []
  [no-slip-v]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_y
    function = 0
  []
  [no-slip-w]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_z
    function = 0
  []
  [outlet_p]
    type = INSFVOutletPressureBC
    boundary = bottom
    variable = pressure
    function = ${p_outlet}
  []
[]

[FunctorMaterials]
  [fluidprops]
    type = GeneralFunctorFluidProps
    mu_rampdown = mu_rampdown
    fp = 'simple_fluid'
    porosity = 1.0
    characteristic_length = 1.0
    T_fluid = 'T'
    pressure = 'pressure'
    speed = 1.0
  []
  [scalar_props]
    type = ADGenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '0.47'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = T
    rho = ${rho_f}
    block = 'box'
  []
  [solid_functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = cp_s
    prop_values = ${cs}
  []
  [solid_ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = Ts
    cp = cp_s
    rho = ${rho_s}
    h = h_solid
    rho_h = rho_h_solid
    block = 'solid'
  []
[]

[Functions]
  [mu_rampdown]
    type = ParsedFunction
    expression = 0.001002*(1+exp(-3*t))
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e6
    viscosity = ${mu}
    density0 =  ${rho_f}
    cp = ${cf}
    thermal_expansion = 0
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
  line_search = 'none'
  l_tol = 1e-8
  nl_rel_tol = 1e-8
  dt = 0.01
  end_time = 1200
[]

[Outputs]
  exodus = true
  csv = true
[]

Also when I quickly tested it with a very high viscosity, it seemed not converging.

[GiudGiud]

ok so I think the way forward is going to be turning off the energy equations (by commenting out the kernels and making the variables AuxVariables instead).
And working on that reduced problem until it converges.

[JuliaHLee]

Thanks for suggesting the way forward. But I could not manage to make it converged (increasing non linear residual). I have commented out the energy equation kernels. But I kept "FVDiffusion" for the 'solid' blocks, because I got an error message that I should include an active kernel for each block. Is there a way to avoid the error and remove the kernel for the solid block? My input file is:

T_0 = 283.15
#T_in = 303.15
p_outlet = 0
lambda_sol = 1.5
cs = 840
rho_s = 2558
#lambda_fl = 0.59
cf = 4180
rho_f = 1000
mu = 0.001002
vf = 5e-05
#hsf = 1e+08


[Mesh]
  [mesh_file]
    type = FileMeshGenerator
    file = mesh_spheres_in_box.msh
  []
  [sideset]
    input = mesh_file
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 'box'
    paired_block = 'solid'
    new_boundary = sideset
  []
[]


[GlobalParams]
  rhie_chow_user_object = rc
[]
[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    w = vel_z
    block = 'box'
    pressure = pressure
  []
[]
[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 1e-10
    face_interp_method = 'skewness-corrected'
  []
  [vel_y]
    type = INSFVVelocityVariable
    block = 'box'
    initial_condition = 1e-10
    face_interp_method = 'skewness-corrected'
  []
  [vel_z]
    type = INSFVVelocityVariable
    block = box
    initial_condition = ${fparse -1*vf}
    face_interp_method = 'skewness-corrected'
  []
  [pressure]
    type = INSFVPressureVariable
    block = 'box'
    face_interp_method = 'skewness-corrected'
  []
  [Ts]
    type = INSFVEnergyVariable
    initial_condition = ${T_0}
    block = 'solid'
    face_interp_method = 'skewness-corrected'
  []
[]

[AuxVariables]
  [velocity_norm]
    type = MooseVariableFVReal
  []
  [T]
    type = INSFVEnergyVariable
    block = 'box'
    initial_condition = ${T_0}
    face_interp_method = 'skewness-corrected'
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    rho = ${rho_f}
    advected_interp_method = 'skewness-corrected'
  []
  [u_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    rho = ${rho_f}
    momentum_component = x
    block = 'box'
    advected_interp_method = 'skewness-corrected'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = x
    block = 'box'
    variable_interp_method = 'skewness-corrected'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = x
    pressure = pressure
    block = 'box'
    correct_skewness = 'true'
  []
  [v_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    rho = ${rho_f}
    momentum_component = y
    block = 'box'
    advected_interp_method = 'skewness-corrected'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = y
    block = 'box'
    variable_interp_method = 'skewness-corrected'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = y
    pressure = pressure
    block = 'box'
    correct_skewness = 'true'
  []
  [w_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
  []
  [w_advection]
    type = INSFVMomentumAdvection
    variable = vel_z
    rho = ${rho_f}
    momentum_component = z
    block = box
    advected_interp_method = 'skewness-corrected'
  []
  [w_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_z
    mu = ${mu}
    momentum_component = z
    block = box
    variable_interp_method = 'skewness-corrected'
  []
  [w_pressure]
    type = INSFVMomentumPressure
    variable = vel_z
    momentum_component = z
    pressure = pressure
    block = box
    correct_skewness = 'true'
  []
#  [energy_time]
#    type = INSFVEnergyTimeDerivative
#    variable = T
#    rho = ${rho_f}
#    dh_dt = dh_dt
#    block = 'box'
#  []
#  [temp_advection]
#    type = INSFVEnergyAdvection
#    variable = T
#    block = 'box'
#    advected_interp_method = 'skewness-corrected'
#  []
#  [temp_conduction]
#    type = FVDiffusion
#    coeff = ${lambda_fl}
#    variable = T
#    block = 'box'
#    variable_interp_method = 'skewness-corrected'
#  []
#  [grain_energy_time]
#    type = INSFVEnergyTimeDerivative
#    variable = Ts
#    rho = ${rho_s}
#    dh_dt = dh_solid_dt
#    block = 'solid'
#  []
  [grain_conduction]
    type = FVDiffusion
    coeff = ${lambda_sol}
    variable = Ts
    block = 'solid'
    variable_interp_method = 'skewness-corrected'
  []
[]
#[FVInterfaceKernels]
#  [convection]
#    type = FVConvectionCorrelationInterface
#    variable1 = T
#    variable2 = Ts
#    boundary = sideset
#    h = ${hsf}
#    T_solid = Ts
#    T_fluid = T
#    subdomain1 = 'box'
#    subdomain2 = 'solid'
#    wall_cell_is_bulk = true
#  []
#[]
[FVBCs]
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_x
    functor = 1e-10
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_y
    functor = 1e-10
  []
  [inlet-w]
    type = INSFVInletVelocityBC
    boundary = top
    variable = vel_z
    functor = ${fparse -1*vf}
  []
#  [inlet_T]
#   type = FVDirichletBC
#    variable = T
#    boundary = top
#    value = ${T_in}
#  []
  [no-slip-u]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_x
    function = 0
  []
  [no-slip-v]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_y
    function = 0
  []
  [no-slip-w]
    type = INSFVNoSlipWallBC
    boundary = 'left right front back sideset'
    variable = vel_z
    function = 0
  []
  [outlet_p]
    type = INSFVOutletPressureBC
    boundary = bottom
    variable = pressure
    function = ${p_outlet}
  []
[]

[FunctorMaterials]
  [fluidprops]
    type = GeneralFunctorFluidProps
    mu_rampdown = mu_rampdown
    fp = 'simple_fluid'
    porosity = 1.0
    characteristic_length = 1.0
    T_fluid = 'T'
    pressure = 'pressure'
    speed = 1.0
    block = 'box'
  []
  [scalar_props]
    type = ADGenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '0.47'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = T
    rho = ${rho_f}
    block = 'box'
  []
  [solid_functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = cp_s
    prop_values = ${cs}
  []
  [solid_ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = Ts
    cp = cp_s
    rho = ${rho_s}
    h = h_solid
    rho_h = rho_h_solid
    block = 'solid'
  []
[]

[Functions]
  [mu_rampdown]
    type = ParsedFunction
    expression = 0.001002*(1+exp(-3*t))
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1e6
    viscosity = ${mu}
    density0 =  ${rho_f}
    cp = ${cf}
    thermal_expansion = 0
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
  line_search = 'none'
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  dt = 0.01
  end_time = 1200
[]

[Outputs]
  exodus = true
  csv = true
[]

[GiudGiud]

But I kept "FVDiffusion" for the 'solid' blocks, because I got an error message that I should include an active kernel for each block. Is there a way to avoid the error and remove the kernel for the solid block?

You can turn off that check with Problem/kernel_coverage_check = false

Ideally you would be down to the minimum number of variables defined on the least number of blocks, get convergence, then ramp up the size of the system until you identify the breaking point

[JuliaHLee]

Hello,

Thank you for the suggestions. I am trying to make the flow model work without energy equations. But I noticed that it converges when I turn off the two_term_boundary_expansion and do not include "skewness correction". So when I add one of these two to fix the problem of high pressure concentration (the original issue), it does not converge. And when I tested the input file with a different mesh with a simpler structure (27 spheres), it converges without the high pressure issue.

Could you give me advice what I could test to get forward? Could the mesh quality be issue or should the "skewness correction" be included only in few blocks? Or should the "derivative-size" in the configuration be adjusted (changing from 120, which I used)?

[GiudGiud]

Just a note that when this happens

Time Step 49, time = 49, dt = 1
 0 Nonlinear |R| = 8.655340e-08
 Solve Converged!

the time step is essentially skipped

Would you suggest me to keep the setting with skewness correction and upwind and then adjust mesh for the next step?

Let's try to figure out mesh settings and interpolation settings that get you the desired result at steady state
Having a very high velocity in 1 time step is not necessarily wrong, because the initialization of the problem is unphysical so the time evolution back to the correct steady solution could include a spike in pressure or velocity.

For my problem, I wanted to have packed spheres, not floating spheres in the 3D space.

I think you can try to make the walls closer to the spheres first, then move the inlet, then the outlet. To see what makes things not converge

[JuliaHLee]

Thanks a lot for the suggestions and your help! I think I will need time to figure out the appropriate mesh setting and interpolation setting. I guess my problem requires fine adjustments due to the internal boundary and small space between boundaries. I wonder if there is another method in MOOSE, which could fit better to this problem.

[GiudGiud]

You could try a porous media approach?
A packing of pebble spheres sounds great!

[JuliaHLee]

I would like to keep the spheres defined separately from the pore space though. You mean the porous media approach with "homogenized" porous media (one phase), right?

[GiudGiud]

yes. The homogenized porous media could represent a pebble bed of packed sphere. This is a common configuration in the nuclear engineering world so there are correlations for it.
If you are located in the US or are collaborating with a US institution you can ask for access to Pronghorn that has some of these correlations already implemented