NS with FE not converging

[adigc]

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 am having difficulties with convergence of a 3D CFD simulation with INSAD at higher Reynolds Number.
I am providing the input I am using below.

NOTE: The mesh contains only linear tetra elements.

With the above settings, that is with

inlet_velocity = 1
mu = 1
rho = 1

Reynolds number = ~60

the simulation takes 8 nonlinear iterations to converge. Everything is fine. But as soon as I increase the inlet velocity to 10 (Reynolds number = 600) the simulation takes forever to converge.

Any hints what I am missing would be helpful.
Thanks in advance.

Additional information

Mesh size and type: 250K nodes, 4 noded linear tet elements.
Reynolds number: Target Reynolds number is around ~4000
Discretization (finite element CG/DG, finite volume, etc): Finite element CG
Models (turbulence, porous media, etc): Not including any turbulence models.
Solver method (fully coupled, segregated, multiapps, etc): fully coupled (input provided.)
Base input you started from:

# input for cfd sim

inlet_boundary = 'inlet'

wall_boundary = 'wall_bottom_top wall_sides '

pressure_boundary = 'outlet'

inlet_velocity = 1

mu = 1

rho = 1.0

[GlobalParams]
  integrate_p_by_parts = true
  viscous_form = traction
[]


[Mesh]
  file = mesh.inp
  construct_side_list_from_node_list=true
[]

[Variables]
  [./velocity]
    order = SECOND
    family = LAGRANGE_VEC
  [../]
  [./p]
  [../]
[]


[AuxVariables]
  [rvel_x]
    # Velocity residual in (x) directions
    order = SECOND
  [] 
  [rvel_y]
    # Velocity residual in (y) directions
    order = SECOND
  [] 
  [rvel_z]
    # Velocity residual in (z) directions
    order = SECOND
  []   
[]


[AuxKernels] 
  [reaction_vx]
    type = DebugResidualAux
    variable = rvel_x
    debug_variable = velocity
    component = 0
  []
  [reaction_vy]
    type = DebugResidualAux
    variable = rvel_y
    debug_variable = velocity
    component = 1
  []
  [reaction_vz]
    type = DebugResidualAux
    variable = rvel_z
    debug_variable = velocity
    component = 2
  []  
[]



[Kernels]
  [./mass]
    type = INSADMass
    variable = p
  [../]
  [./momentum_convection]
    type = INSADMomentumAdvection
    variable = velocity
  [../]
  [./momentum_viscous]
    type = INSADMomentumViscous
    variable = velocity
  [../]
  [./momentum_pressure]
    type = INSADMomentumPressure
    variable = velocity
    pressure = p
  [../]
  [momentum_supg]
    type = INSADMomentumSUPG
    variable = velocity
    velocity = velocity
  []  
[]

[BCs]
  [./inlet]
    type = VectorDirichletBC
    variable = velocity
    boundary = ${inlet_boundary}
    values = '0.0 ${inlet_velocity} 0.0'
  [../]
  [./wall]
    type = VectorDirichletBC
    variable = velocity
    boundary = ${wall_boundary}
    values = '0.0 0.0 0.0'
  [../]

  [./pressure_pin]
    type = DirichletBC
    variable = p
    boundary = ${pressure_boundary}
    value = 0.0
  [../]
[]

[Materials]
  [./const]
    type = ADGenericConstantMaterial
    prop_names = 'rho mu'
    prop_values = '${rho} ${mu}'
  [../]
  [ins_mat]
    type = INSADTauMaterial
    velocity = velocity
    pressure = p
  []
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
    solve_type = 'NEWTON'
  [../]
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
 
  snesmf_reuse_base = false
  petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'preonly   lu       mumps'  
  nl_rel_tol = 1e-6
  nl_max_its = 100
  l_max_its = 1000
  automatic_scaling = false
[]


[Postprocessors]
  [reaction_x]
    type = NodalSum
    variable = rvel_x
    boundary = ${wall_boundary}
  []
  [reaction_y]
    type = NodalSum
    variable = rvel_y
    boundary = ${wall_boundary}
  []
  [reaction_z]
    type = NodalSum
    variable = rvel_z
    boundary = ${wall_boundary}
  []  
[]


[Outputs]
  [fe_pattern]
    type = Exodus
    execute_on = 'nonlinear'
  []
  perf_graph = false
  csv = true
  print_nonlinear_residuals = true
  print_linear_residuals = true  
[]

[GiudGiud]

Hello

Adding a turbulence model will facilitate convergence.
A few things you could try:

• initializing the Re=600 case with the Re=60 solution
• ramping up the velocity slowly using a Function and turning this simulation into a transient

[adigc]

You'll need to code one, we never merged a turbulence model in INSAD. The mixing length we have in FV would be simple enough to add. You can look at that one too for an example https://mooseframework.inl.gov/source/kernels/INSADSmagorinskyEddyViscosity.html

well I used INSAD only because I saw that example. If just INSFE or INSFV makes it easy I will go for it !

[adigc]

Also I though the SUPG will help dealing with high Re simulations ! Isn't that the case ?

[GiudGiud]

Yes it helps. And you can try PSPG as well if you think that will be enough.
We are developing INSFV at the moment

[adigc]

Thank you for the suggestions.
I believe I have progress.

So I tried

ramping up the velocity slowly using a Function and turning this simulation into a transient

together with https://mooseframework.inl.gov/source/kernels/INSADSmagorinskyEddyViscosity.html

the simulation is more stable now. What puzzles me now is that simulation now hangs after some 200 time steps.

    340 Linear |R| = 8.285654e-09
    341 Linear |R| = 8.043995e-09
    342 Linear |R| = 7.831560e-09
 1 Nonlinear |R| = 2.301015e-04

    Computing Jacobian                                                                   [  5.24 s] [ 1273 MB]

Just like that. It isn't proceeding further !

[GiudGiud]

Can you use a debugger to attach to the process, then print a backtrace.

you should see where it is hanging.
My guess is that the linear solve is failing. It would only report a failure at 10000 linear iterations

  petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'preonly   lu       mumps'  

this should have been solving in only 1 linear iteration if you did not use preonly can you try that?