Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Mass flow rates not respected in WCNSFV #29709

Open
tanoret opened this issue Jan 20, 2025 · 0 comments · May be fixed by #29710
Open

Mass flow rates not respected in WCNSFV #29709

tanoret opened this issue Jan 20, 2025 · 0 comments · May be fixed by #29710
Labels
P: normal A defect affecting operation with a low possibility of significantly affects. T: defect An anomaly, which is anything that deviates from expectations.

Comments

@tanoret
Copy link
Contributor

tanoret commented Jan 20, 2025

Bug Description

The issue appears when using variable density materials in Navier Stokes, e.g., density dependent on temperature.

When specifying inlet mass flow rates using the WCNSFVMassFluxBC and WCNSFVMomentumFluxBC, the specified mass flow rates are reconstructed via second order interpolation to the face and are not necessarily reproducing the values specified by the user. This could be fine for the outlet, where the density have to be reconstructed to the outlet face, but not for the inlet where the user is specifying the mass flow rate.

Steps to Reproduce

This is an example input reproducing the error. The mass flow rate should be 1,400 kg/s. However, the post-processors are showing the following values:

+----------------+----------------+----------------+----------------+----------------+----------------+
| time | area_pp_left | inlet_mdot | inlet_mfr | outlet_mfr | velocity_pp |
+----------------+----------------+----------------+----------------+----------------+----------------+
: : : : : : :
| 7.167000e+00 | 1.000000e+00 | 1.400000e+03 | -1.391791e+03 | 1.635018e+03 | 1.000000e+00 |
| 8.191000e+00 | 1.000000e+00 | 1.400000e+03 | -1.390112e+03 | 1.575618e+03 | 1.000000e+00 |
| 1.023900e+01 | 1.000000e+00 | 1.400000e+03 | -1.388269e+03 | 1.511694e+03 | 1.000000e+00 |
| 1.433500e+01 | 1.000000e+00 | 1.400000e+03 | -1.385957e+03 | 1.478058e+03 | 1.000000e+00 |
| 2.252700e+01 | 1.000000e+00 | 1.400000e+03 | -1.384653e+03 | 1.471348e+03 | 1.000000e+00 |
| 3.891100e+01 | 1.000000e+00 | 1.400000e+03 | -1.384485e+03 | 1.471105e+03 | 1.000000e+00 |
| 7.167900e+01 | 1.000000e+00 | 1.400000e+03 | -1.384581e+03 | 1.471155e+03 | 1.000000e+00 |
| 1.372150e+02 | 1.000000e+00 | 1.400000e+03 | -1.384645e+03 | 1.471159e+03 | 1.000000e+00 |
| 2.682870e+02 | 1.000000e+00 | 1.400000e+03 | -1.384676e+03 | 1.471159e+03 | 1.000000e+00 |
| 5.304310e+02 | 1.000000e+00 | 1.400000e+03 | -1.384691e+03 | 1.471159e+03 | 1.000000e+00 |
| 1.054719e+03 | 1.000000e+00 | 1.400000e+03 | -1.384699e+03 | 1.471159e+03 | 1.000000e+00 |
| 2.103295e+03 | 1.000000e+00 | 1.400000e+03 | -1.384703e+03 | 1.471159e+03 | 1.000000e+00 |
| 4.200447e+03 | 1.000000e+00 | 1.400000e+03 | -1.384705e+03 | 1.471159e+03 | 1.000000e+00 |
| 8.394751e+03 | 1.000000e+00 | 1.400000e+03 | -1.384705e+03 | 1.471159e+03 | 1.000000e+00 |
| 1.000000e+04 | 1.000000e+00 | 1.400000e+03 | -1.384704e+03 | 1.471159e+03 | 1.000000e+00 |
+----------------+----------------+----------------+----------------+----------------+----------------+

The computed inlet mass flow rate is 1,384.7 kg/s at the inlet and 1,4711.5 at the outlet. This is wrong for both cases. The errors reduce as we refine the mesh, of course, but we would generally like to use coarse meshes.

rho = 'rho'
l = 10
velocity_interp_method = 'rc'
advected_interp_method = 'average'

# Artificial fluid properties
# For a real case, use a GeneralFluidFunctorProperties and a viscosity rampdown
k = 1
cp = 1000
mu = 1e2

# Operating conditions
inlet_temp = 100
outlet_pressure = 1e5
inlet_v = 1.0

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${l}
    ymin = 0
    ymax = 1
    nx = 10
    ny = 4
  []
[]

[GlobalParams]
  rhie_chow_user_object = 'rc'
[]

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

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = ${inlet_v}
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 1e-15
  []
  [pressure]
    type = INSFVPressureVariable
    initial_condition = ${outlet_pressure}
  []
  [T_fluid]
    type = INSFVEnergyVariable
    initial_condition = ${inlet_temp}
  []
[]

[AuxVariables]
  [mixing_length]
    type = MooseVariableFVReal
  []
  [power_density]
    type = MooseVariableFVReal
    initial_condition = 1e8
  []
[]

[FVKernels]
  inactive = 'u_turb v_turb temp_turb'
  [mass_time]
    type = WCNSFVMassTimeDerivative
    variable = pressure
    drho_dt = drho_dt
  []
  [mass]
    type = WCNSFVMassAdvection
    variable = pressure
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
  []

  [u_time]
    type = WCNSFVMomentumTimeDerivative
    variable = vel_x
    drho_dt = drho_dt
    rho = rho
    momentum_component = 'x'
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_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
  []
  [u_turb]
    type = INSFVMixingLengthReynoldsStress
    variable = vel_x
    rho = ${rho}
    mixing_length = 'mixing_length'
    momentum_component = 'x'
    u = vel_x
    v = vel_y
  []

  [v_time]
    type = WCNSFVMomentumTimeDerivative
    variable = vel_y
    drho_dt = drho_dt
    rho = rho
    momentum_component = 'y'
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_interp_method}
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    momentum_component = 'y'
    mu = ${mu}
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  []
  [v_turb]
    type = INSFVMixingLengthReynoldsStress
    variable = vel_y
    rho = ${rho}
    mixing_length = 'mixing_length'
    momentum_component = 'y'
    u = vel_x
    v = vel_y
  []

  [temp_time]
    type = WCNSFVEnergyTimeDerivative
    variable = T_fluid
    rho = rho
    drho_dt = drho_dt
    h = h
    dh_dt = dh_dt
  []
  [temp_conduction]
    type = FVDiffusion
    coeff = 'k'
    variable = T_fluid
  []
  [temp_advection]
    type = INSFVEnergyAdvection
    variable = T_fluid
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_interp_method}
  []
  [heat_source]
    type = FVCoupledForce
    variable = T_fluid
    v = power_density
  []
  [temp_turb]
    type = WCNSFVMixingLengthEnergyDiffusion
    variable = T_fluid
    rho = rho
    cp = cp
    mixing_length = 'mixing_length'
    schmidt_number = 1
    u = vel_x
    v = vel_y
  []
[]

[FVBCs]
  [slip_x]
    type = INSFVNaturalFreeSlipBC
    variable = vel_x
    boundary = 'top bottom'
    momentum_component = 'x'
  []
  [slip_y]
    type = INSFVNaturalFreeSlipBC
    variable = vel_y
    boundary = 'top bottom'
    momentum_component = 'y'
  []

  # Inlet
  [inlet_mass]
    type = WCNSFVMassFluxBC
    variable = pressure
    boundary = 'left'
    mdot_pp = 'inlet_mdot'
    # velocity_pp = 'velocity_pp'
    area_pp = 'area_pp_left'
    rho = 'rho'
    vel_x = vel_x
    vel_y = vel_y
  []
  [inlet_u]
    type = WCNSFVMomentumFluxBC
    variable = vel_x
    boundary = 'left'
    mdot_pp = 'inlet_mdot'
    # velocity_pp = 'velocity_pp'
    area_pp = 'area_pp_left'
    rho = 'rho'
    momentum_component = 'x'
    vel_x = vel_x
    vel_y = vel_y
  []
  [inlet_v]
    type = WCNSFVMomentumFluxBC
    variable = vel_y
    boundary = 'left'
    mdot_pp = 0
    # velocity_pp = 0
    area_pp = 'area_pp_left'
    rho = 'rho'
    momentum_component = 'y'
    vel_x = vel_x
    vel_y = vel_y
  []
  [inlet_T]
    type = FVDirichletBC
    variable = T_fluid
    boundary = 'left'
    value = ${inlet_temp}
  []

  [outlet_p]
    type = INSFVOutletPressureBC
    variable = pressure
    boundary = 'right'
    function = ${outlet_pressure}
  []
[]

[Functions]
  [k_func]
    type = ParsedFunction
    expression = '${k}'
  []
  [rho_func]
    type = ParsedFunction
    expression = '1.5e3 - x'
  []
  [mu_func]
    type = ParsedFunction
    expression = '${mu}'
  []
[]

[FluidProperties]
  # [fp]
  #   type = FlibeFluidProperties
  # []
  [fp]
    type = TemperaturePressureFunctionFluidProperties
    cv = ${cp}
    k = k_func
    rho = rho_func
    mu = mu_func
  []
[]

[FunctorMaterials]
  [const_functor]
    type = ADGenericFunctorMaterial
    prop_names = 'cp k'
    prop_values = '${cp} ${k}'
  []
  [rho]
    type = RhoFromPTFunctorMaterial
    fp = fp
    temperature = T_fluid
    pressure = pressure
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    temperature = 'T_fluid'
    rho = ${rho}
  []
[]

[AuxKernels]
  inactive = 'mixing_len'
  [mixing_len]
    type = WallDistanceMixingLengthAux
    walls = 'top'
    variable = mixing_length
    execute_on = 'initial'
    delta = 0.5
  []
[]

[Postprocessors]
  [inlet_mdot]
    type = Receiver
    default = 1.4e3
  []
  [velocity_pp]
    type = Receiver
    default = ${inlet_v}
  []
  [area_pp_left]
    type = AreaPostprocessor
    boundary = 'left'
    execute_on = 'INITIAL'
  []
  [inlet_mfr]
    type = VolumetricFlowRate
    boundary = left
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 'rho'
  []
  [outlet_mfr]
    type = VolumetricFlowRate
    boundary = right
    vel_x = vel_x
    vel_y = vel_y
    advected_quantity = 'rho'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-3
    optimal_iterations = 6
  []
  end_time = 1e4

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

Impact

This will improve the domain-overlap coupling action and enable the user to actually see in the postprocessors the values that were actually specified for the inlet mass flow rates.

[Optional] Diagnostics

No response
@tanoret tanoret added P: normal A defect affecting operation with a low possibility of significantly affects. T: defect An anomaly, which is anything that deviates from expectations. labels Jan 20, 2025
tanoret pushed a commit to tanoret/moose that referenced this issue Jan 20, 2025
adding dirichlet BC for inlet mass flow rate Refs idaholab#29709
acf1c1d
@tanoret tanoret linked a pull request Jan 20, 2025 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
P: normal A defect affecting operation with a low possibility of significantly affects. T: defect An anomaly, which is anything that deviates from expectations.
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Adding dirichlet BC for inlet mass flow rate tanoret/moose
1 participant
@tanoret