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]

average interpolation for advected quantity is unstable, you'll want to use upwind

I also want to ask a question, I used the boundary conditions set by the NSFVHeatFluxBC class, but I found that the temperature rise fed back to me was less than the value I expected, maybe it was about 10,000 times worse than the value to match the corresponding $W/m^2$ result, is the unit of heat flow given in this class what I thought?

[GiudGiud]

is the unit of heat flow given in this class what I thought?

The heat flux is in the units of the mesh and your fluid properties. If you want W/m2 make sure you are using SI units everywhere.

I used the boundary conditions set by the NSFVHeatFluxBC class

NSFVHeatFlux is meant for splitting the heat flux between the fluid and the solid. Make sure to use it with a porosity of 1 such that al l the heat flux goes to the fluid.
Also make sure to select the fluid phase.

If you already know the value of the heat flux, you might as well pick a FVNeumannBC

[shin-chan-cpu]

Thank you GiudGiud
I solved this problem by adjusting the boundary conditions, but there was a new problem when I computed the k-e model, Refer to https://github.com/idaholab/moose/blob/next/modules/navier_stokes/test/tests/finite_volume/ins/turbulence/channel/channel_ERCOFTAC.i input card and add the skewness correction option,What makes me strange is that when I use Cartesian grid calculation, the result is still normal, and the three equations converge to a good degree (even if the turbulent two equations do not converge to a good degree), but when I use my geometry, the result of the calculation becomes very bad, I tried to refer to the hypre tutorial on the official website to modify the petsc option but it feels not very useful, do you have any suggestions?

[GiudGiud]

I tried to refer to the hypre tutorial on the official website to modify the petsc option but it feels not very usefuf

hypre wont be great for advection unless there is so much turbulence the problem is more diffusive than advective.
You could:

• stick with LU for the whole system of 4 equations (you dont have energy right?)
• work on designing a field split for this 4 equation problem
  https://mooseframework.inl.gov/source/splits/Split.html
  https://mooseframework.inl.gov/source/preconditioners/FieldSplitPreconditioner.html
• wait for us to add k-epsilon to the linear solve and use that. The linear FV solver for fluid flow is seeing a lot of development and we are aiming to have this by the end of the month

[shin-chan-cpu]

you dont have energy right?

I've added the energy equation, but after joining, I found that my whole system can't converge, and if I don't open the energy equation, his convergence effect is still good.