INSFVMomentumViscousSourceRZ

INSFVMomentumDiffusion provides the Laplacian portion of the viscous stress. In a cylindrical, axisymmetric coordinate system (no swirl) the radial momentum equation also contains an additional source term that comes from $var element = document.getElementById("moose-equation-1625f14e-2fcc-4287-afe2-20d6d081ef44");katex.render("\\nabla^2 \\boldsymbol{u}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ :

$var element = document.getElementById("moose-equation-6cbcbb87-984b-46d7-871b-f8ec0636915c");katex.render("-\\mu \\frac{u_r}{r^2}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

INSFVMomentumViscousSourceRZ adds this volumetric contribution so that both the residual and the Rhie-Chow $var element = document.getElementById("moose-equation-2ec21ea2-3ff6-487f-b8cc-b2d85521726e");katex.render("a", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ -coefficients include the effect of the $var element = document.getElementById("moose-equation-701ce472-1443-4dd2-a857-3292c7b756c4");katex.render("- \\mu u_r / r^2", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ term. The radius is evaluated from the element vertex average (cell centroid) to match the cell-centered discretization. This kernel must be restricted to blocks that use an RZ coordinate system and to the radial momentum component (the component selected with "momentum_component" must be the axisymmetric radial index).

note:Automatically added for WCNSFVFlowPhysics

When using Navier Stokes Flow / WCNSFVFlowPhysics, this kernel is added automatically on blocks that use an RZ coordinate system. Set "add_rz_viscous_source" to false if the default behavior needs to be disabled.

note

The kernel expects the same viscosity functor that is passed to INSFVMomentumDiffusion, and it only needs to be added for the radial momentum equation. Axial components should not include this kernel.

note

When the user enables turbulence, in some cases, kernels can be split into contributions from the molecular and turbulent viscosity. In these cases, the user has to make sure to include one of these kernels for each INSFVMomentumDiffusion!

note:Use only for RZ viscous sources

INSFVMomentumViscousSourceRZ complements the Laplacian viscous term implemented by INSFVMomentumDiffusion. Do not add this kernel in Cartesian problems; only include it for the radial momentum equation of axisymmetric $var element = document.getElementById("moose-equation-d29f57cf-08f7-442e-a9bc-95870e7aea29");katex.render("RZ", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ setups. When the INSFVMomentumDiffusion kernel uses the complete-expansion form, set "complete_expansion" so the same factor of 2 is applied here.

Input Parameters

momentum_componentThe component of the momentum equation that this kernel applies to.
C++ Type:MooseEnum
Options:x, y, z
Controllable:No
Description:The component of the momentum equation that this kernel applies to.
muThe dynamic viscosity that multiplies the viscous source term. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The dynamic viscosity that multiplies the viscous source term. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
rhie_chow_user_objectThe rhie-chow user-object
C++ Type:UserObjectName
Controllable:No
Description:The rhie-chow user-object
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
complete_expansionFalseMirror of the INSFVMomentumDiffusion 'complete_expansion' switch. When true, the viscous contribution is multiplied by 2.
Default:False
C++ Type:bool
Controllable:No
Description:Mirror of the INSFVMomentumDiffusion 'complete_expansion' switch. When true, the viscous contribution is multiplied by 2.

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

ghost_layers1The number of layers of elements to ghost.
Default:1
C++ Type:unsigned short
Controllable:No
Description:The number of layers of elements to ghost.
use_point_neighborsFalseWhether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.

Parallel Ghosting Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/cylindrical/no-slip-tris.i)
(modules/navier_stokes/test/tests/finite_volume/ins/mms/cylindrical/2d-rc-rz-viscous-source.i)

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/cylindrical/no-slip-tris.i)

mu = 1
rho = 1
advected_interp_method = 'average'
velocity_interp_method = 'rc'

[GlobalParams]
  two_term_boundary_expansion = true
  rhie_chow_user_object = 'rc'
[]

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

[Mesh]
  type = GeneratedMesh
  nx = 4
  ny = 4
  xmax = 3.9
  ymax = 4.1
  elem_type = TRI3
  dim = 2
  coord_type = 'RZ'
[]

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = 1e-15
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 1e-15
  []
  [pressure]
    type = INSFVPressureVariable
  []
[]

[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_viscous_source]
    type = INSFVMomentumViscousSourceRZ
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
    pressure = pressure
  []

  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'y'
    # we can think of the axis as a slip wall boundary, no normal velocity and no viscous shear
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  []
[]

[FVBCs]
  active = 'inlet-u inlet-v free-slip-wall-u free-slip-wall-v outlet-p axis-u axis-v'

  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = 'bottom'
    variable = vel_x
    functor = 0
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = 'bottom'
    variable = vel_y
    functor = 1
  []
  [free-slip-wall-u]
    type = INSFVNaturalFreeSlipBC
    boundary = 'right'
    variable = vel_x
    momentum_component = 'x'
  []
  [free-slip-wall-v]
    type = INSFVNaturalFreeSlipBC
    boundary = 'right'
    variable = vel_y
    momentum_component = 'y'
  []
  [no-slip-wall-u]
    type = INSFVNoSlipWallBC
    boundary = 'right'
    variable = vel_x
    function = 0
  []
  [no-slip-wall-v]
    type = INSFVNoSlipWallBC
    boundary = 'right'
    variable = vel_y
    function = 0
  []
  [outlet-p]
    type = INSFVOutletPressureBC
    boundary = 'top'
    variable = pressure
    function = 0
  []
  [axis-u]
    type = INSFVSymmetryVelocityBC
    boundary = 'left'
    variable = vel_x
    u = vel_x
    v = vel_y
    mu = ${mu}
    momentum_component = x
  []
  [axis-v]
    type = INSFVSymmetryVelocityBC
    boundary = 'left'
    variable = vel_y
    u = vel_x
    v = vel_y
    mu = ${mu}
    momentum_component = y
  []
  [axis-p]
    type = INSFVSymmetryPressureBC
    boundary = 'left'
    variable = pressure
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
  nl_rel_tol = 1e-12
[]

[Debug]
  show_var_residual_norms = true
[]

[Postprocessors]
  [in]
    type = SideIntegralVariablePostprocessor
    variable = vel_y
    boundary = 'bottom'
  []
  [out]
    type = SideIntegralVariablePostprocessor
    variable = vel_y
    boundary = 'top'
  []
  [num_lin]
    type = NumLinearIterations
    outputs = 'console'
  []
  [num_nl]
    type = NumNonlinearIterations
    outputs = 'console'
  []
  [cum_lin]
    type = CumulativeValuePostprocessor
    outputs = 'console'
    postprocessor = 'num_lin'
  []
  [cum_nl]
    type = CumulativeValuePostprocessor
    outputs = 'console'
    postprocessor = 'num_nl'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

(modules/navier_stokes/test/tests/finite_volume/ins/mms/cylindrical/2d-rc-rz-viscous-source.i)

mu_ref = 1.2
mu_r = 0.1
mu_x = 0.2
rho_ref = 1.4
rho_r = 0.15
rho_x = 0.25
advected_interp_method = 'average'
velocity_interp_method = 'rc'

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  []
  coord_type = 'RZ'
  rz_coord_axis = x
[]

[GlobalParams]
  rhie_chow_user_object = rc
[]

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = 0.0
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 0.0
  []
  [pressure]
    type = INSFVPressureVariable
    initial_condition = 0.0
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
[]

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

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    rho = rho_fun
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
  []
  [pressure_pin]
    type = FVPointValueConstraint
    lambda = lambda
    phi0 = 0.125
    point = '0.5 0.5 0'
    variable = pressure
  []

  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = rho_fun
    momentum_component = 'x'
  []
  [u_diffusion]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = mu_fun
    momentum_component = 'x'
    complete_expansion = true
    u = vel_x
    v = vel_y
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
    pressure = pressure
  []
  [u_body]
    type = INSFVBodyForce
    variable = vel_x
    functor = forcing_u
    momentum_component = 'x'
  []

  [v_diffusion]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = mu_fun
    momentum_component = 'y'
    complete_expansion = true
    u = vel_x
    v = vel_y
  []
  [v_viscous_source]
    type = INSFVMomentumViscousSourceRZ
    variable = vel_y
    mu = mu_fun
    momentum_component = 'y'
    complete_expansion = true
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = rho_fun
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  []
  [v_body]
    type = INSFVBodyForce
    variable = vel_y
    functor = forcing_v
    momentum_component = 'y'
  []
[]

[FVBCs]
  [axis_u]
    type = INSFVSymmetryVelocityBC
    boundary = 'bottom'
    variable = vel_x
    u = vel_x
    v = vel_y
    mu = mu_fun
    momentum_component = x
  []
  [axis_v]
    type = INSFVSymmetryVelocityBC
    boundary = 'bottom'
    variable = vel_y
    u = vel_x
    v = vel_y
    mu = mu_fun
    momentum_component = y
  []
  [axis_p]
    type = INSFVSymmetryPressureBC
    boundary = 'bottom'
    variable = pressure
  []
  [dirichlet_u_wall]
    type = INSFVNoSlipWallBC
    boundary = 'left right top'
    variable = vel_x
    function = exact_u
  []
  [dirichlet_v_wall]
    type = INSFVNoSlipWallBC
    boundary = 'left right top'
    variable = vel_y
    function = exact_v
  []
[]

[Functions]
  # --------------------------------------------------
  # Material property functions
  # --------------------------------------------------
  [rho_fun]
    type = ParsedFunction
    symbol_names = 'rho_ref rho_r rho_x'
    symbol_values = '${rho_ref} ${rho_r} ${rho_x}'
    expression = 'rho_ref + rho_x*x + rho_r*y'
  []
  [mu_fun]
    type = ParsedFunction
    symbol_names = 'mu_ref mu_r mu_x'
    symbol_values = '${mu_ref} ${mu_r} ${mu_x}'
    expression = 'mu_ref + mu_x*x + mu_r*y'
  []

  # --------------------------------------------------
  # Convenience reciprocals
  # --------------------------------------------------
  [S]
    type=ParsedFunction
    symbol_names = 'rho_fun'
    symbol_values = 'rho_fun'
    expression='1.0/rho_fun'
  []
  [S2]
    type=ParsedFunction
    symbol_names = 'rho_fun'
    symbol_values = 'rho_fun'
    expression='(1.0/rho_fun)^2'
  []
  [S3]
    type=ParsedFunction
    symbol_names = 'rho_fun'
    symbol_values = 'rho_fun'
    expression='(1.0/rho_fun)^3'
  []

  # --------------------------------------------------
  # Building blocks
  # --------------------------------------------------
  [A]
    type=ParsedFunction
    expression='x^2*(1-x)^2'
  []
  [Ap]
    type=ParsedFunction
    expression='2*x - 6*x^2 + 4*x^3'
  []
  [App]
    type=ParsedFunction
    expression='2 - 12*x + 12*x^2'
  []
  [Appp]
    type=ParsedFunction
    expression='-12 + 24*x'
  []
  [R]
    type=ParsedFunction
    expression='y^3 - 2*y^4 + y^5'
  []
  [Rp]
    type=ParsedFunction
    expression='3*y^2 - 8*y^3 + 5*y^4'
  []
  [Rpp]
    type=ParsedFunction
    expression='6*y - 24*y^2 + 20*y^3'
  []
  [Rppp]
    type=ParsedFunction
    expression='6 - 48*y + 60*y^2'
  []

  # --------------------------------------------------
  # Exact solutions (from mass streamfunction)
  # --------------------------------------------------
  [exact_u]
    type=ParsedFunction
    symbol_names = 'A Rp S'
    symbol_values = 'A Rp S'
    expression='A*Rp*S / y'
  []
  [exact_v]
    type=ParsedFunction
    symbol_names = 'Ap R S'
    symbol_values = 'Ap R S'
    expression='-Ap*R*S / y'
  []
  [exact_p]
    type=ParsedFunction
    expression='x*y^2'
  []

  # --------------------------------------------------
  # First derivatives
  # --------------------------------------------------
  [ux]
    type=ParsedFunction
    symbol_names = 'Ap Rp A S rho_x S2'
    symbol_values = 'Ap Rp A S ${rho_x} S2'
    expression='(Ap*Rp*S - A*Rp*rho_x*S2)/y'
  []
  [ur]
    type=ParsedFunction
    symbol_names = 'A Rpp Rp S rho_r S2'
    symbol_values = 'A Rpp Rp S ${rho_r} S2'
    expression='(A*Rpp*S - A*Rp*rho_r*S2)/y - (A*Rp*S)/y^2'
  []
  [vx]
    type=ParsedFunction
    symbol_names = 'App R Ap S rho_x S2'
    symbol_values = 'App R Ap S ${rho_x} S2'
    expression='-(App*R*S - Ap*R*rho_x*S2)/y'
  []
  [vr]
    type=ParsedFunction
    symbol_names = 'Ap Rp R S rho_r S2'
    symbol_values = 'Ap Rp R S ${rho_r} S2'
    expression='(-Ap*Rp*S + Ap*R*rho_r*S2)/y + (Ap*R*S)/y^2'
  []

  # --------------------------------------------------
  # Divergence (for stress)
  # --------------------------------------------------
  [theta]
    type=ParsedFunction
    symbol_names = 'A Rp Ap R rho_x rho_r S2'
    symbol_values = 'A Rp Ap R ${rho_x} ${rho_r} S2'
    expression='( -A*Rp*rho_x + Ap*R*rho_r )*S2 / y'
  []

  [T]
    type=ParsedFunction
    symbol_names = 'A Rp Ap R rho_x rho_r'
    symbol_values = 'A Rp Ap R ${rho_x} ${rho_r}'
    expression='-A*Rp*rho_x + Ap*R*rho_r'
  []
  [thetax]
    type=ParsedFunction
    symbol_names = 'Ap Rp App R rho_x rho_r S2 S3 T'
    symbol_values = 'Ap Rp App R ${rho_x} ${rho_r} S2 S3 T'
    expression='( (-Ap*Rp*rho_x + App*R*rho_r)*S2 - 2*rho_x*T*S3 )/y'
  []
  [thetar]
    type=ParsedFunction
    symbol_names = 'A Rpp Ap Rp rho_x rho_r S2 S3 T'
    symbol_values = 'A Rpp Ap Rp ${rho_x} ${rho_r} S2 S3 T'
    expression='( (-A*Rpp*rho_x + Ap*Rp*rho_r)*S2 - 2*rho_r*T*S3 )/y - (T*S2)/y^2'
  []

  # --------------------------------------------------
  # Second derivatives
  # --------------------------------------------------
  [uxx]
    type=ParsedFunction
    symbol_names = 'App Rp Ap rho_x A S S2 S3'
    symbol_values = 'App Rp Ap ${rho_x} A S S2 S3'
    expression='( App*Rp*S - 2*Ap*Rp*rho_x*S2 + 2*A*Rp*(rho_x^2)*S3 )/y'
  []
  [uxr]
    type=ParsedFunction
    symbol_names = 'Ap Rpp Rp rho_r rho_x A S S2 S3'
    symbol_values = 'Ap Rpp Rp ${rho_r} ${rho_x} A S S2 S3'
    expression='( Ap*Rpp*S - Ap*Rp*rho_r*S2 - A*Rpp*rho_x*S2 + 2*A*Rp*(rho_x*rho_r)*S3 )/y - (Ap*Rp*S - A*Rp*rho_x*S2)/y^2'
  []
  [urr]
    type=ParsedFunction
    symbol_names = 'A Rppp Rpp Rp rho_r S S2 S3'
    symbol_values = 'A Rppp Rpp Rp ${rho_r} S S2 S3'
    expression='( A*Rppp*S - 2*A*Rpp*rho_r*S2 + 2*A*Rp*(rho_r^2)*S3 )/y - 2*(A*Rpp*S - A*Rp*rho_r*S2)/y^2 + 2*(A*Rp*S)/y^3'
  []
  [vxx]
    type=ParsedFunction
    symbol_names = 'Appp R App rho_x Ap S S2 S3'
    symbol_values = 'Appp R App ${rho_x} Ap S S2 S3'
    expression='( -Appp*R*S + 2*App*R*rho_x*S2 - 2*Ap*R*(rho_x^2)*S3 )/y'
  []
  [vxr]
    type=ParsedFunction
    symbol_names = 'App Rp R rho_r rho_x Ap S S2 S3'
    symbol_values = 'App Rp R ${rho_r} ${rho_x} Ap S S2 S3'
    expression='( -App*Rp*S + App*R*rho_r*S2 + Ap*Rp*rho_x*S2 - 2*Ap*R*(rho_x*rho_r)*S3 )/y + (App*R*S - Ap*R*rho_x*S2)/y^2'
  []
  [vrr]
    type=ParsedFunction
    symbol_names = 'Ap Rpp R Rp rho_r S S2 S3'
    symbol_values = 'Ap Rpp R Rp ${rho_r} S S2 S3'
    expression='( -Ap*Rpp*S + 2*Ap*Rp*rho_r*S2 - 2*Ap*R*(rho_r^2)*S3 )/y + 2*(Ap*Rp*S - Ap*R*rho_r*S2)/y^2 - 2*(Ap*R*S)/y^3'
  []

  # --------------------------------------------------
  # Stresses
  # --------------------------------------------------
  [tau_xx]
    type=ParsedFunction
    symbol_names = 'mu_fun ux'
    symbol_values = 'mu_fun ux'
    expression='2*mu_fun*ux'
  []
  [tau_rr]
    type=ParsedFunction
    symbol_names = 'mu_fun vr'
    symbol_values = 'mu_fun vr'
    expression='2*mu_fun*vr'
  []
  [tau_xr]
    type=ParsedFunction
    symbol_names = 'mu_fun ur vx'
    symbol_values = 'mu_fun ur vx'
    expression='mu_fun*(ur + vx)'
  []
  [tau_tt]
    type=ParsedFunction
    symbol_names = 'mu_fun exact_v'
    symbol_values = 'mu_fun exact_v'
    expression='2*mu_fun*(exact_v/y)'
  []

  # --------------------------------------------------
  # Stress derivatives (for symmetric-only stress)
  # --------------------------------------------------
  [tau_xx_x]
    type = ParsedFunction
    symbol_names  = 'mu_fun ux uxx mu_x'
    symbol_values = 'mu_fun ux uxx ${mu_x}'
    expression = '2*(mu_x*ux + mu_fun*uxx)'
  []

  [tau_xr_x]
    type = ParsedFunction
    symbol_names  = 'mu_fun ur vx uxr vxx mu_x'
    symbol_values = 'mu_fun ur vx uxr vxx ${mu_x}'
    expression = 'mu_x*(ur + vx) + mu_fun*(uxr + vxx)'
  []

  [tau_xr_r]
    type = ParsedFunction
    symbol_names  = 'mu_fun ur vx urr vxr mu_r'
    symbol_values = 'mu_fun ur vx urr vxr ${mu_r}'
    expression = 'mu_r*(ur + vx) + mu_fun*(urr + vxr)'
  []

  [tau_rr_r]
    type = ParsedFunction
    symbol_names  = 'mu_fun vr vrr mu_r'
    symbol_values = 'mu_fun vr vrr ${mu_r}'
    expression = '2*(mu_r*vr + mu_fun*vrr)'
  []

  # --------------------------------------------------
  # Convective fluxes
  # --------------------------------------------------
  [Fxx]
    type=ParsedFunction
    symbol_names = 'rho_fun exact_u'
    symbol_values = 'rho_fun exact_u'
    expression='rho_fun*exact_u*exact_u'
  []
  [Fxr]
    type=ParsedFunction
    symbol_names = 'rho_fun exact_u exact_v'
    symbol_values = 'rho_fun exact_u exact_v'
    expression='rho_fun*exact_u*exact_v'
  []
  [Frx]
    type=ParsedFunction
    symbol_names = 'rho_fun exact_u exact_v'
    symbol_values = 'rho_fun exact_u exact_v'
    expression='rho_fun*exact_v*exact_u'
  []
  [Frr]
    type=ParsedFunction
    symbol_names = 'rho_fun exact_v'
    symbol_values = 'rho_fun exact_v'
    expression='rho_fun*exact_v*exact_v'
  []

  # --------------------------------------------------
  # Conservative flux derivatives (product rule)
  # --------------------------------------------------
  [Fxx_x]
    type = ParsedFunction
    symbol_names  = 'rho_fun exact_u ux rho_x'
    symbol_values = 'rho_fun exact_u ux ${rho_x}'
    expression = 'rho_x*exact_u*exact_u + rho_fun*(2*exact_u*ux)'
  []

  [Frx_r]
    type = ParsedFunction
    symbol_names  = 'rho_fun exact_u exact_v ur vr rho_r'
    symbol_values = 'rho_fun exact_u exact_v ur vr ${rho_r}'
    expression = 'rho_r*exact_v*exact_u + rho_fun*(vr*exact_u + exact_v*ur)'
  []

  [Fxr_x]
    type = ParsedFunction
    symbol_names  = 'rho_fun exact_u exact_v ux vx rho_x'
    symbol_values = 'rho_fun exact_u exact_v ux vx ${rho_x}'
    expression = 'rho_x*exact_u*exact_v + rho_fun*(ux*exact_v + exact_u*vx)'
  []

  [Frr_r]
    type = ParsedFunction
    symbol_names  = 'rho_fun exact_v vr rho_r'
    symbol_values = 'rho_fun exact_v vr ${rho_r}'
    expression = 'rho_r*exact_v*exact_v + rho_fun*(2*exact_v*vr)'
  []

  # --------------------------------------------------
  # Pressure gradients
  # --------------------------------------------------
  [px]
    type=ParsedFunction
    expression='y^2'
  []
  [pr]
    type=ParsedFunction
    expression='2*x*y'
  []

  # --------------------------------------------------
  # Forcing terms (momentum residuals)
  # --------------------------------------------------
  [forcing_u]
    type=ParsedFunction
    symbol_names = 'Fxx_x Frx Frx_r px tau_xx_x tau_xr tau_xr_r'
    symbol_values = 'Fxx_x Frx Frx_r px tau_xx_x tau_xr tau_xr_r'
    expression='Fxx_x + (Frx/y) + Frx_r + px - ( tau_xx_x + (tau_xr/y) + tau_xr_r )'
  []
  [forcing_v]
    type=ParsedFunction
    symbol_names = 'Fxr_x Frr Frr_r pr tau_xr_x tau_rr tau_rr_r tau_tt'
    symbol_values = 'Fxr_x Frr Frr_r pr tau_xr_x tau_rr tau_rr_r tau_tt'
    expression='Fxr_x + (Frr/y) + Frr_r + pr - ( tau_xr_x + (tau_rr/y) + tau_rr_r - (tau_tt)/y )'
  []
[]

[AuxVariables]
  [vel_x_aux]
    type = MooseLinearVariableFVReal
    initial_condition = 0.0
  []
  [vel_y_aux]
    type = MooseLinearVariableFVReal
    initial_condition = 0.0
  []
  [pressure_aux]
    type = MooseLinearVariableFVReal
    initial_condition = 0.0
  []
[]

[AuxKernels]
  [vel_x_function_aux]
    type = FunctionAux
    variable = vel_x_aux
    function = exact_u
    execute_on = TIMESTEP_END
  []
  [vel_y_function_aux]
    type = FunctionAux
    variable = vel_y_aux
    function = exact_v
    execute_on = TIMESTEP_END
  []
  [pressure_function_aux]
    type = FunctionAux
    variable = pressure_aux
    function = exact_p
    execute_on = TIMESTEP_END
  []
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -ksp_gmres_restart'
  petsc_options_value = 'lu 200'
  line_search = 'none'
  nl_abs_tol = 1e-12
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = 'timestep_end'
  []
  [L2u]
    type = ElementL2Error
    variable = vel_x
    function = exact_u
    execute_on = 'timestep_end'
  []
  [L2v]
    type = ElementL2Error
    variable = vel_y
    function = exact_v
    execute_on = 'timestep_end'
  []
  [L2p]
    type = ElementL2Error
    variable = pressure
    function = exact_p
    execute_on = 'timestep_end'
  []
[]

[Outputs]
  csv = true
  exodus = true
  execute_on = timestep_end
[]

Input Parameters
Input Files