PINSFV convergence issue

[lucuo]

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 team, I'm using the porous media flow module in the NS module to simulate a 3D heat exchanger. I currently have two issues:
① When using the Darcy drag coefficient, it requires extremely large values to converge. I suspect small coefficients might cause turbulent flow instability, and though I've added the Mixlength turbulence model, what other methods can improve convergence?
② When using the Forchheimer drag coefficient, it not converge - what causes this？ and What can I do to improve convergence
Attached are my input file,

mu = 0.001003
rho = 998.2
cp = 4182
u_inlet = 3
T_inlet = 373
advected_interp_method = 'upwind'#'average'#'upwind'
velocity_interp_method = 'rc'
von_karman_const = 0.21
[Mesh]
  [gen]
    type = FileMeshGenerator
    file = exchange_3.e
  []
[]
[GlobalParams]
  rhie_chow_user_object = 'rc'
[]
[Problem]
  kernel_coverage_check = false
  fv_bcs_integrity_check = true
[]
[UserObjects]
  [rc]
    type = PINSFVRhieChowInterpolator
    u = superficial_vel_x
    v = superficial_vel_y
    w = superficial_vel_z
    pressure = pressure
    porosity = porosity
  []
[]
[Variables]
  [superficial_vel_x]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1
  []
  [superficial_vel_y]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 0.5
  []
  [superficial_vel_z]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1
  []  
  [pressure]
    type = INSFVPressureVariable
    initial_condition = 1e5
  []
  [T_fluid]
    type = INSFVEnergyVariable
  []
[]
[Functions]
  [vel_xx]
    type = ParsedFunction
    expression = 'if(-0.4<z<0.4,if(0.62<x<1.22,0.5,1),0.3)'
  []
  [vel_zz]
    type = ParsedFunction
    expression = 'if(0.5<x<1.35,if(0.75<x<1.15,0.2,-0.5),-2)'
  []
[]
[ICs]
  [vel_xx]
    type = FunctionIC
    variable = 'superficial_vel_x'
    function = vel_xx
  []
  [vel_zz]
    type = FunctionIC
    variable = 'superficial_vel_z'
    function = vel_zz
  []
[]
[AuxVariables]
  [T_solid]
    family = 'MONOMIAL'
    order = 'CONSTANT'
    fv = true
    initial_condition = 293
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 0.5
  []
  [mixing_len]
    order = CONSTANT
    family = MONOMIAL
    fv = true
  []
  [Reaux]
    type = MooseVariableFVReal
  []
  [Refunctor]
    type = MooseVariableFVReal
  []
[]
[FVKernels]
  [mass]
    type = PINSFVMassAdvection
    variable = pressure
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
  []
  [u_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    porosity = porosity
    momentum_component = 'x'
  []
  [u_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_x
    mu = ${mu}
    porosity = porosity
    momentum_component = 'x'
  []
  [u_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_x
    momentum_component = 'x'
    pressure = pressure
    porosity = porosity
  []

  [v_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    porosity = porosity
    momentum_component = 'y'
  []
  [v_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_y
    mu = ${mu}
    porosity = porosity
    momentum_component = 'y'
  []
  [v_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_y
    momentum_component = 'y'
    pressure = pressure
    porosity = porosity
  []

  [w_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_z
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    porosity = porosity
    momentum_component = 'z'
  []
  [w_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_z
    mu = ${mu}
    porosity = porosity
    momentum_component = 'z'
  []
  [w_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_z
    momentum_component = 'z'
    pressure = pressure
    porosity = porosity
  []
  [energy_advection]
    type = PINSFVEnergyAdvection
    variable = T_fluid
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_interp_method}
  []
  [energy_diffusion_1]
    type = PINSFVEnergyAnisotropicDiffusion
    kappa = 'kappa'
    variable = T_fluid
    porosity = porosity
  []
  [energy_convection]
    type = PINSFVEnergyAmbientConvection
    variable = T_fluid
    is_solid = false
    T_fluid = T_fluid
    T_solid = T_solid
    h_solid_fluid = 'h_cv'
  []
###friction
  [u_friction]
    type = PINSFVMomentumFriction
    variable = superficial_vel_x
    momentum_component = 'x'
    Darcy_name = 'Darcy_coefficient'
    Forchheimer_name = 'Forchheimer_coefficient'
#    porosity = 0.5
    mu = ${mu}
    rho = ${rho}
    speed = speed
  []
  [v_friction]
    type = PINSFVMomentumFriction
    variable = superficial_vel_y
    momentum_component = 'y'
    Darcy_name = 'Darcy_coefficient'
    Forchheimer_name = 'Forchheimer_coefficient'
    rho = ${rho}
    speed = speed
    mu = ${mu}
  []
  [w_friction]
    type = PINSFVMomentumFriction
    variable = superficial_vel_z
    momentum_component = 'z'
    Darcy_name = 'Darcy_coefficient'
    Forchheimer_name = 'Forchheimer_coefficient'
    rho = ${rho}
    speed = speed
    mu = ${mu}
  []
###turb
  [u_viscosity_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_x
    rho = ${rho}
    mixing_length = mixing_len
    momentum_component = 'x'
    u = superficial_vel_x
    v = superficial_vel_y
    w = superficial_vel_z
  []
  [v_viscosity_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_y
    rho = ${rho}
    mixing_length = mixing_len
    momentum_component = 'y'
    u = superficial_vel_x
    v = superficial_vel_y
    w = superficial_vel_z
  []
  [w_viscosity_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_z
    rho = ${rho}
    mixing_length = mixing_len
    momentum_component = 'z'
    u = superficial_vel_x
    v = superficial_vel_y
    w = superficial_vel_z
  []
[]
[AuxKernels]
  [mixing_len]
    type = WallDistanceMixingLengthAux
    walls = 'wall_shell wall_baffle wall_baffle_top'
    variable = mixing_len
    execute_on = 'initial'
    von_karman_const = ${von_karman_const}
  []
  [Reynolds]
    type = ReynoldsNumberFunctorAux
    variable = Reaux
    speed = speed
    rho = ${rho}
    mu = ${mu}
  []
  [Refunctor]
    type = FunctorAux
    functor = 'Re'
    variable = Refunctor
    execute_on = 'timestep_begin'
  []
[]
[FVBCs]
  inactive ='inlet-T' #'inlet-T-dirichlet' #'inlet-T'#'inlet-T-dirichlet'
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = superficial_vel_x
    function = 0
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = superficial_vel_y
    function = 0
  []
  [inlet-w]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = superficial_vel_z
    function = -${u_inlet}
  []  
  [inlet-T]
    type = FVNeumannBC
    variable = T_fluid
    value = '${fparse u_inlet * rho * cp * T_inlet}'
    boundary = 'inlet'
  []
  [inlet-T-dirichlet]
    type = FVDirichletBC
    variable = T_fluid
    value = '${T_inlet}'
    boundary = 'inlet'
  []
  [free-slip-u]
    type = INSFVNaturalFreeSlipBC
    boundary = 'wall_shell wall_baffle wall_baffle_top'
    variable = superficial_vel_x
    momentum_component = 'x'
  []
  [free-slip-v]
    type = INSFVNaturalFreeSlipBC
    boundary = 'wall_shell wall_baffle wall_baffle_top'
    variable = superficial_vel_y
    momentum_component = 'y'
  []
  [free-slip-w]
    type = INSFVNaturalFreeSlipBC
    boundary = 'wall_shell wall_baffle wall_baffle_top'
    variable = superficial_vel_x
    momentum_component = 'z'
  []
  [outlet-p]
    type = INSFVOutletPressureBC
    boundary = 'outlet'
    variable = pressure
    function = 100000
  []
[]
[FunctorMaterials]
  [constants]
    type = ADGenericFunctorMaterial
    prop_names = 'h_cv'
    prop_values = '1000000'
  []
  [functor_constants]
    type = ADGenericFunctorMaterial
    prop_names = 'cp'
    prop_values = '${cp}'
  []
  [kappa]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'kappa'
    prop_values = '0.6 0.6 0.6'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    rho = ${rho}
    temperature = 'T_fluid'
  []
  [darcy]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'Darcy_coefficient Forchheimer_coefficient'
    prop_values = '6e7 6e7 6e7 0 0 0'
  []
  [speed]
    type = PINSFVSpeedFunctorMaterial
    superficial_vel_x = superficial_vel_x
    superficial_vel_y = superficial_vel_y
    superficial_vel_z = superficial_vel_z
    porosity = porosity
  []
  [Re_material]
    type = ReynoldsNumberFunctorMaterial
    speed = speed
    characteristic_length = 2
    rho = ${rho}
    mu = ${mu}
  []
[]

[Preconditioning]  
 [./SMP_NEWTON]                            
    type = SMP
    full = true
    solve_type = 'NEWTON'
petsc_options_iname = '-pc_mat_solver_package -pc_type -pc_factor_shift_type'
petsc_options_value ='mumps asm NONZERO'
[]
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-ksp_gmres_restart -pc_hypre_boomeramg_max_iter -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 = '201  4 hypre    boomeramg    0.7    4    5    25    HMIS    ext+i    2    0.3'
  line_search = 'none'
  nl_rel_tol = 1e-8                            
  nl_abs_tol = 1e-8
  nl_max_its = 300
  nl_forced_its = 5
  l_tol = 1e-8                                   
  l_max_its = 300
[]
# Some basic Postprocessors to examine the solution
[Postprocessors]
  [inlet-p]
    type = SideAverageValue
    variable = pressure
    boundary = 'inlet'
  []
  [outlet-u]
    type = SideAverageValue
    variable = superficial_vel_y
    boundary = 'outlet'
  []
  [outlet-temp]
    type = SideAverageValue
    variable = T_fluid
    boundary = 'outlet'
  []
  [solid-temp]
    type = ElementAverageValue
    variable = T_solid
  []
[]

[Outputs]
  file_base = output/ex_03_water_turb
  exodus = true
  csv = false
[]

mesh file,
Desktop.zip
mesh schematic."

Additional information

Mesh size and type:
Reynolds number: 1.9e3~3e5
Discretization (finite element CG/DG, finite volume, etc):finite volume
Models (turbulence, porous media, etc): porous media and mixinglength
Solver method (fully coupled, segregated, multiapps, etc): newton
Base input you started from:

mu = 0.001003
rho = 998.2
cp = 4182
u_inlet = 3
T_inlet = 373

[GiudGiud]

Hello

I recommend you set the preconditioning parameters only once in the Preconditioning block, and you should start with a direct solve to get a solution first.
I would also recommend solving just the flow equations first, and adding the energy equation later once it solves without the energy equation.

You can also explore using the SIMPLENonLinearAssembly executioner rather than the Newton solver

Guillaume