Use the NS module to calculate fluid occurrence

[shin-chan-cpu]

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

I used Newton method to calculate the fluid flow in the pipeline. Since no three-dimensional fluid calculation example with energy conservation equation was given on the official website, I could only modify it through two-dimensional calculation example, but the convergence was poor. That is, only the physical properties given by the example can be convergent, and once they are changed, errors or non-convergence will be reported, and they will also exist.I wonder why this is, is it the iteration itself?

Additional information

Mesh size and type:Hexahedral grid, about 20,000 meshes
Reynolds number:2.2e8
Discretization (finite element CG/DG, finite volume, etc):FV
Models (turbulence, porous media, etc):noTurb
Solver method (fully coupled, segregated, multiapps, etc):Newton
Base input you started from:

# mu = 1.1
# rho = 1.1
# k = 1e-3
# cp = 1.1
#alpha = 150
mu = 1e-6
rho = 1000
k = 0.5
cp = 4222
# Operating conditions
w_inlet = 5
T_inlet = 300
#T_solid = 300
p_outlet = 15.5e6
#h_fs = 1


advected_interp_method = 'upwind'
velocity_interp_method = 'average'

#pressure_tag = "pressure_grad"

[Mesh]
  # [mesh]
  #   type = CartesianMeshGenerator
  #   dim = 3
  #   dx = '2'
  #   dy = '2'
  #   dz = '8'
  #   ix = '15'
  #   iy = '15'
  #   iz = '20'
  # []
  [mesh]
    type = FileMeshGenerator
    show_info = true
    file = 1h_sff.e
[]
[]

[GlobalParams]
  rhie_chow_user_object = 'rc'
[]

# [Problem]
#   nl_sys_names = 'u_system v_system w_system pressure_system energy_system'
#   previous_nl_solution_required = true
#   error_on_jacobian_nonzero_reallocation = true
# []

[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    w = vel_z
    pressure = pressure
  []
[]

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = 1e-12
  #  solver_sys = u_system
  #  two_term_boundary_expansion = false
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 1e-12
  #  solver_sys = v_system
  #  two_term_boundary_expansion = false
  []
  [vel_z]
    type = INSFVVelocityVariable
    initial_condition = ${w_inlet}
  #  solver_sys = w_system
  #  two_term_boundary_expansion = false
  []
  [pressure]
    type = INSFVPressureVariable
  #  solver_sys = pressure_system
    initial_condition = ${p_outlet}
  #  two_term_boundary_expansion = false
  []
  [T_fluid]
    type = INSFVEnergyVariable
    initial_condition = ${T_inlet}
  #  solver_sys = energy_system
  #  two_term_boundary_expansion = false
  []
  [T_solid]
#    type = INSFVEnergyVariable
    order = CONSTANT
    family = MONOMIAL
    fv = true
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
  []

  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'x'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
    pressure = pressure
  #  extra_vector_tags = ${pressure_tag}
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  #  extra_vector_tags = ${pressure_tag}
  []
  [w_advection]
    type = INSFVMomentumAdvection
    variable = vel_z
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'z'
  []
  [w_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_z
    mu = ${mu}
    momentum_component = 'z'
  []
  [w_pressure]
    type = INSFVMomentumPressure
    variable = vel_z
    momentum_component = 'z'
    pressure = pressure
  #  extra_vector_tags = ${pressure_tag}
  []
  # [p_diffusion]
  #   type = FVAnisotropicDiffusion
  #   variable = pressure
  #   coeff = "Ainv"
  #   coeff_interp_method = 'average'
  # []
  # [p_source]
  #   type = FVDivergence
  #   variable = pressure
  #   vector_field = "HbyA"
  #   force_boundary_execution = true
  # []
  [energy_advection]
    type = INSFVEnergyAdvection
    variable = T_fluid
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_interp_method}
  []
  [energy_diffusion]
    type = FVDiffusion
    coeff = ${k}
    variable = T_fluid
  []
  # [ambient_convection]
  #   type = NSFVEnergyAmbientConvection
  #   variable = T_fluid
  #   T_ambient = 350
  #   alpha = 'alpha'
  # []
  # [energy_convection]
  #   type = PINSFVEnergyAmbientConvection
  #   variable = T_fluid
  #   is_solid = false
  #   T_fluid = 'T_fluid'
  #   T_solid = 'T_solid'
  #   h_solid_fluid = 'h_cv'
  # []  
  [solid_temp_conduction]
    type = FVDiffusion
    coeff = 1
    variable = T_solid
  []
[]

[FVBCs]
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = vel_x
    function = '0'
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = vel_y
    function = '0'
  []
  [inlet-w]
    type = INSFVInletVelocityBC
    boundary = 'inlet'
    variable = vel_z
    function = ${w_inlet}
  []
  [walls-u]
    type = INSFVNoSlipWallBC
    boundary = 'walls   '
    variable = vel_x
    function = 0.0
  []
  [walls-v]
    type = INSFVNoSlipWallBC
    boundary = 'walls '
    variable = vel_y
    function = 0.0
  []
  [walls-w]
    type = INSFVNoSlipWallBC
    boundary = 'walls '
    variable = vel_z
    function = 0.0
  []

  [symmetry-u]
    type = INSFVSymmetryVelocityBC
    boundary = 'symmetry'
    variable = vel_x
    u = vel_x
    v = vel_y
    w = vel_z
    mu = ${mu}
    momentum_component = 'x'
  []
  [symmetry-v]
    type = INSFVSymmetryVelocityBC
    boundary = 'symmetry'
    variable = vel_y
    u = vel_x
    v = vel_y
    w = vel_z
    mu = ${mu}
    momentum_component = 'y'
  []
  [symmetry-w]
    type = INSFVSymmetryVelocityBC
    boundary = 'symmetry'
    variable = vel_z
    u = vel_x
    v = vel_y
    w = vel_z
    mu = ${mu}
    momentum_component = 'z'
  []  
  [symmetry-p]
    type = INSFVSymmetryPressureBC
    boundary = 'symmetry'
    variable = pressure
  []

  [outlet_p]
    type = INSFVOutletPressureBC
    boundary = 'outlet'
    variable = pressure
    function = '${p_outlet}'
  []
  # [zero-grad-pressure]
  #   type = FVFunctionNeumannBC
  #   variable = pressure
  #   boundary = 'back left right top bottom'
  #   function = 0.0
  # []
  [inlet_t]
    type = FVDirichletBC
    boundary = 'inlet'
    variable = T_fluid
    value = ${T_inlet}
  []

 # [left_Ts]
 #   type = NSFVHeatFluxBC
 #   variable = T_fluid
 #   boundary = 'walls '
 #   phase = 'fluid'
 #   value = 100000
 #   splitting = thermal_conductivity
  #  locality = global
 #   average_k_fluid='0.04'
  #  average_k_solid='500'
 #   average_kappa='0.04'  # because of vector matprop, should be kappa
 #   average_kappa_solid='500'   
 # []
  # [heater]
  #   type = FVDirichletBC
  #   variable = 'T_solid'
  #   boundary = 'left right top bottom'
  #   value = 300
  # []
[]

# [FVInterfaceKernels]
#   [convection]
#     type = FVConvectionCorrelationInterface
#     variable1 = T_fluid
#     variable2 = T_solid
#     boundary = 'left right top bottom'
#     h = 0.1
#     T_solid = T_solid
#     T_fluid = T_fluid
#     subdomain1 = 0
#     subdomain2 = 1
#     wall_cell_is_bulk = true
#   []
# []

[FunctorMaterials]
  # [constants]
  #   type = ADGenericFunctorMaterial
  #   prop_names = 'h_cv T_solid'
  #   prop_values = '${h_fs} ${T_solid}'
  # [] 
  [const_functor]
    type = ADGenericFunctorMaterial
    prop_names = 'cp'
    prop_values = '${cp}'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    rho = ${rho}
    temperature = 'T_fluid'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
#  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
#  petsc_options_value = '201    hypre   boomeramg    4'  
  nl_rel_tol = 7e-5
  nl_abs_tol = 1e-9
  nl_max_its = 50
  line_search = 'none'
  automatic_scaling = true
  off_diagonals_in_auto_scaling = true
  compute_scaling_once = false
[]

[Outputs]
  exodus = true
  csv = true
  file_base = output/f_constant
  # perf_graph = false
  # print_nonlinear_residuals = false
  # print_linear_residuals = true
[]

[GiudGiud]

Hello

With a Reynolds number of 2.2e8 you need a turbulence model to converge. The additional viscosity from the turbulence will help lower the effective Reynolds number

[shin-chan-cpu]

Under theoretical calculation, the pressure drop should be within 1 MPa

[GiudGiud]

can you plot the pressure axial profile? See if it's oscillating maybe

[shin-chan-cpu]

I have outputted the pressure distribution along the axis, which is a straight line sloping downward without any oscillations.

[shin-chan-cpu]

I also encountered a problem where such a warning appears during the calculation. After this, my calculation results have a very small residual, but the deviation is significant. Do you know why this happens?

[GiudGiud]

The pressure and temperature must be so out of range that the correlation for viscosity evaluates to negative values.

This makes no sense because you set viscosity to a constant in the input you pasted.

This

velocity_interp_method = 'average'

should be set to

velocity_interp_method = 'rc'

I would also switch to using a regular cartesian mesh and get correct results. Once that works you can use your mesh file.