idaholab / moose

Multiphysics Object Oriented Simulation Environment
https://www.mooseframework.org
GNU Lesser General Public License v2.1
1.71k stars 1.04k forks source link

Navier Stokes Compressible Flow HLLC with RZ #24753

Open snschune opened 1 year ago

snschune commented 1 year ago

Bug Description

Incorrect pressure at R=0 boundary when using Navier-Stokes compressible HLLC implementation.

Steps to Reproduce

length = '${units 2 ft -> m}'

D_wick_i = '${units 0.65 in -> m}'
R_wick_i = '${fparse D_wick_i/2}'
A_core = '${fparse 0.25 * pi * D_wick_i^2}'
# P_core = '${fparse pi * D_wick_i}'

# T_initial = '${units 540.8 degF -> K}'
# T_source = '${units 540.8 degF -> K}'
# p_source = '${units 38.9 psi -> Pa}'
# p_sink = '${units 24.5 psi -> Pa}'

# m_dot_source = -0.068039
# m_dot_sink = 0.068039
# vel = 135

#rho_in = 2.35
#H_in = 2504225
#gamma = 1.4
#R = 8.3145
#molar_mass = 29e-3
#R_specific = 1718
#cv = 4290
#cp = '${fparse cv * gamma}'

T_in = '${units 540.8 degF -> K}'
m_dot = 0.068039

#mass_flux = ${fparse inlet_vel * rho_in}
mass_flux = '${fparse m_dot / A_core}'

#outlet_pressure = 168921.6

[Debug]
  show_material_props = true
[]

[Mesh]
  coord_type = 'RZ'
  rz_coord_axis = 'X'
  [tubo_2d]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${length}
    nx = 150
    ymin = 0
    ymax = ${R_wick_i}
    ny = 25
    bias_y = '${fparse 1 / 1.1}'
  []
  [rename]
    type = RenameBlockGenerator
    input = 'tubo_2d'
    old_block = '0'
    new_block = 'evap'
  []
  [mesh_modified1]
    type = SubdomainBoundingBoxGenerator
    input = 'rename'
    block_id = '2'
    block_name = 'cond'
    bottom_left = '${fparse 0.5 * length} 0 0'
    top_right = '${length} ${R_wick_i} 0'
  []
  [new_boundary]
    type = BreakBoundaryOnSubdomainGenerator
    input = 'mesh_modified1'
    boundaries = 'top bottom'
  []
[]

[FluidProperties]
  [fp]
    type = IdealGasFluidProperties
    gamma = 1.4
    molar_mass = 0.029
    allow_imperfect_jacobians = true
  []
[]
[GlobalParams]
  fp = fp
[]

[Variables]
  [rho]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 2.35
  []

  [rho_u]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 1e-15
    outputs = none
  []

  [rho_v]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 317
    outputs = none
  []

  [rho_E]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 5616723
  []
[]

[ICs]
[]

[FVKernels]
  # Mass conservation
  [mass_time]
    type = FVTimeKernel
    variable = rho
  []

  [mass_advection]
    type = CNSFVMassHLLC
    variable = rho
    fp = fp
  []

  # Momentum x conservation
  [momentum_x_time]
    type = FVTimeKernel
    variable = rho_u
  []

  [momentum_x_advection]
    type = CNSFVMomentumHLLC
    variable = rho_u
    momentum_component = x
    fp = fp
  []

  # Momentum y conservation
  [momentum_y_time]
    type = FVTimeKernel
    variable = rho_v
  []

  [momentum_y_advection]
    type = CNSFVMomentumHLLC
    variable = rho_v
    momentum_component = y
  []

  # Fluid energy conservation
  [fluid_energy_time]
    type = FVTimeKernel
    variable = rho_E
  []

  [fluid_energy_advection]
    type = CNSFVFluidEnergyHLLC
    variable = rho_E
    fp = fp
  []
[]

[FVBCs]
  [mass_inflow]
    type = CNSFVHLLCSpecifiedMassFluxAndTemperatureMassBC
    variable = rho
    boundary = 'top_to_evap'
    rhou = 0
    rhov = ${mass_flux}
    temperature = ${T_in}
  []

  [momentum_x_inflow]
    type = CNSFVHLLCSpecifiedMassFluxAndTemperatureMomentumBC
    variable = rho_u
    boundary = 'top_to_evap'
    rhou = 0
    rhov = ${mass_flux}
    temperature = ${T_in}
    momentum_component = x
  []
  [momentum_y_inflow]
    type = CNSFVHLLCSpecifiedMassFluxAndTemperatureMomentumBC
    variable = rho_v
    boundary = 'top_to_evap'
    rhou = 0
    rhov = ${mass_flux}
    temperature = ${T_in}
    momentum_component = y
  []
  [fluid_energy_inflow]
    type = CNSFVHLLCSpecifiedMassFluxAndTemperatureFluidEnergyBC
    variable = rho_E
    boundary = 'top_to_evap'
    rhou = 0
    rhov = ${mass_flux}
    temperature = ${T_in}
  []

  ## outflow implicit conditions
  [mass_outflow]
    type = CNSFVHLLCMassImplicitBC
    variable = rho
    fp = fp
    boundary = 'top_to_cond'
  []

  [momentum_x_outflow]
    type = CNSFVHLLCMomentumImplicitBC
    variable = rho_u
    momentum_component = x
    fp = fp
    boundary = 'top_to_cond'
  []

  [momentum_y_outflow]
    type = CNSFVHLLCMomentumImplicitBC
    variable = rho_v
    momentum_component = y
    fp = fp
    boundary = 'top_to_cond'
  []

  [fluid_energy_outflow]
    type = CNSFVHLLCFluidEnergyImplicitBC
    variable = rho_E
    fp = fp
    boundary = 'top_to_cond'
  []
  # wall conditions
  [momentum_x_pressure_wall]
    type = CNSFVMomImplicitPressureBC
    variable = rho_u
    momentum_component = x
    boundary = 'left right'
  []

  [momentum_y_pressure_wall]
    type = CNSFVMomImplicitPressureBC
    variable = rho_v
    momentum_component = y
    boundary = 'left right'
  []
[]

[AuxVariables]
  [Ma]
    family = MONOMIAL
    order = CONSTANT
  []

  [p]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]

  [p_aux]
    type = ADMaterialRealAux
    variable = p
    property = pressure
  []
[]

[Materials]
  [var_mat]
    type = ConservedVarValuesMaterial
    rho = rho
    rhou = rho_u
    rhov = rho_v
    rho_et = rho_E
  []
  [sound_speed]
    type = SoundspeedMat
    fp = fp
  []
[]

[Postprocessors]
  [cfl_dt]
    type = ADCFLTimeStepSize
    c_names = 'sound_speed'
    vel_names = 'speed'
    CFL = 1.25
  []
  # [should_be_one]
  #   type = ElementAverageValue
  #   block = 1
  #   variable = constantVar
  #   execute_on = 'initial timestep_end'
  # []
  [volume1]
    type = VolumePostprocessor
    block = 0
    execute_on = 'initial timestep_end'
  []
  [volume2]
    type = VolumePostprocessor
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu       NONZERO                superlu_dist'
  end_time = 100

  [TimeIntegrator]
    type = ExplicitSSPRungeKutta
    order = 2
  []
  l_tol = 1e-7

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = cfl_dt
  []
[]

[Outputs]
  exodus = true
[]

Impact

how urgent is this? I would suggest to try KT instead of HLLC. ping @joshuahansel

snschune commented 1 year ago

reported by @cbdutra

joshuahansel commented 1 year ago

@snschune It's not urgent. @cbdutra can use the other formulations for now.