INSAD failing to converge

[j-bowhay]

Hi
I am trying to model fluid flow through a straight pipe as a simple benchmark. I have had good success using INSFV but I am now trying to use INSAD and struggling to get the solve to converge. Have tried making the mesh finer and messing around with various options to no avail.

Any suggestions welcome.
Jake

Input File

# 2d siumulation of a water through a pipe.

mu=1e-3 # Nsm^-2
rho=997.0 # kgm^-3
Re=1000.0
pipe_length=1 # m
pipe_radius=0.1 # m
u_inlet=${fparse (mu * Re)/(2 * pipe_radius * rho)} # ms^-1

[GlobalParams]
  integrate_p_by_parts = true
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${pipe_length}
    ymin = 0
    ymax = ${pipe_radius}
    nx = 1000
    ny = 100
  []
[]

[Problem]
  coord_type = 'RZ'
  rz_coord_axis = x
[]

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

[ICs]
  [velocity]
    type = VectorConstantIC
    x_value = 1e-15
    y_value = 1e-15
    variable = velocity
  []
[]

[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
  [../]
[]

[BCs]
  [inlet]
    type = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'left'
    function_x = ${u_inlet}
    function_y = 0
  [../]
  [wall]
    type = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'top'
    function_x = 0
    function_y = 0
  []
  [axis]
    type = ADVectorFunctionDirichletBC
    variable = velocity
    boundary = 'bottom'
    set_x_comp = false
    function_y = 0
  []
  [outlet]
    type = INSADMomentumNoBCBC
    variable = velocity
    pressure = p
    boundary = 'right'
  []
  # When the NoBCBC is applied on the outlet boundary then there is nothing
  # constraining the pressure. Thus we must pin the pressure somewhere to ensure
  # that the problem is not singular. If the below BC is not applied then
  # -pc_type svd -pc_svd_monitor reveals a singular value
  [p_corner]
    type = DirichletBC
    boundary = 'right'
    value = 0
    variable = p
  []
[]

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

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

[Executioner]
  type = Steady

  petsc_options_iname = '-pc_type -sub_pc_factor_levels -ksp_gmres_restart'
  petsc_options_value = 'asm      6                     200'
  nl_rel_tol = 1e-12
  nl_abs_tol = 1e-12
  l_max_its = 200
  line_search = 'none'
[]

[Outputs]
  exodus = true
[]

[GiudGiud]

Hello

What is your Reynolds number? Do you expect the flow to be turbulent?
INSAD, which is based on finite element, requires numerical stabilization techniques such as PSPG or SUPG.

Guillaume

[lindsayad]

I have not, but I think I have some time right now

[lindsayad]

There was a bug in our RZ coordinate handling. See #21070 for the fix. Your input also needed a couple of tweaks. Here is the tweaked input:

# 2d siumulation of a water through a pipe.

mu=1e-3 # Nsm^-2
rho=997.0 # kgm^-3
Re=1000.0
pipe_length=1 # m
pipe_radius=0.1 # m
u_inlet=${fparse (mu * Re)/(2 * pipe_radius * rho)} # ms^-1

[GlobalParams]
  integrate_p_by_parts = false
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${pipe_length}
    ymin = 0
    ymax = ${pipe_radius}
    nx = 1000
    ny = 100
    # elem_type = 'QUAD9'
  []
[]

[Problem]
  coord_type = 'RZ'
  rz_coord_axis = x
[]

[Variables]
  [velocity]
    family = LAGRANGE_VEC
    # order = second
  []
  [p]
  []
[]

[ICs]
  [velocity]
    type = VectorConstantIC
    x_value = 1e-15
    y_value = 1e-15
    variable = velocity
  []
[]

[Kernels]
  [mass]
    type = INSADMass
    variable = p
  []
  [mass_pspg]
    type = INSADMassPSPG
    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 = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'left'
    function_x = ${u_inlet}
    function_y = 0
  [../]
  [wall]
    type = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'top'
    function_x = 0
    function_y = 0
  []
  [axis]
    type = ADVectorFunctionDirichletBC
    variable = velocity
    boundary = 'bottom'
    set_x_comp = false
    function_y = 0
  []
  # [outlet]
  #   type = INSADMomentumNoBCBC
  #   variable = velocity
  #   pressure = p
  #   boundary = 'right'
  # []
  # # When the NoBCBC is applied on the outlet boundary then there is nothing
  # # constraining the pressure. Thus we must pin the pressure somewhere to ensure
  # # that the problem is not singular. If the below BC is not applied then
  # # -pc_type svd -pc_svd_monitor reveals a singular value
  [p_corner]
    type = DirichletBC
    boundary = 'right'
    value = 0
    variable = p
  []
[]

[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

  petsc_options_iname = '-pc_type -sub_pc_factor_levels -ksp_gmres_restart'
  petsc_options_value = 'asm      6                     200'
  nl_rel_tol = 1e-12
  nl_abs_tol = 1e-12
  l_max_its = 200
  line_search = 'none'
[]

[Outputs]
  exodus = true
[]

[lindsayad]

The important input change was to set integrate_p_by_parts = false and then to turn off the outlet BC. These changes force the viscous flux for the x-velocity to be 0 at the axis-of-symmetry and the viscous fluxes at the outlet to be 0 for both x-velocity and y-velocity. We definitely want the x-velocity viscous flux to be 0 along the axis-of-symmetry; that's critical for the accuracy of the solution. The outlet BC change is probably not that important as long as your channel is sufficiently long for the flow to fully develop

[j-bowhay]

Thank you for the very helpful response. Will hopefully have a play with this over the few days.

[j-bowhay]

Finally got a chance to try this and the solution qualitatively looks good. Planning on doing some qualitative validation in the near future.