# ============================================================================== # MODEL PARAMETERS # ============================================================================== v_in = <> # superficial inlet velocity # - bed geometry bed_diam = 0.12065 # bed diameter bed_height = 1.524 # total bed height for experiment bed_vol = ${fparse pi/4*bed_diam^2*bed_height} z3 = 0.762 # elevation of lower pressure tap z5 = 1.2172 # elevation of upper pressure tap dz1 = ${z3} dz2 = ${fparse (z5 - z3)} dz3 = ${fparse (bed_height - z5)} # height difference for pressure drop over length calculation dH = ${fparse z5 - z3} # - bed parameters for D/dp = 19 pebble_diameter = 0.635e-2 eps = 0.3781 # bulk porosity to yield average value of 0.385 num_radial = 10 # dR = 1.9 pebble diams. # - bed parameters for D/dp = 9.5 # pebble_diameter = 1.27e-2 # eps = 0.397 # num_radial = 5 # dR = 1.9 pebble diams. # - bed parameters for D/dp = 6.33 # pebble_diameter = 1.905e-2 # eps = 0.416 # num_radial = 3 # dR = 2.11 pebble diams. # num_radial = 4 # dR = 1.58 pebble diams. mu = 0.00079735 # water @ 30 C rho_fluid = 995.65 T_in = ${fparse 30 + 273.15} # Computation parameters advected_interp_method='upwind' velocity_interp_method='rc' # ============================================================================== # GEOMETRY AND MESH # ============================================================================== [Mesh] [rz_mesh] type = CartesianMeshGenerator dim = 2 dx = '${fparse bed_diam/2}' dy = '${dz1} ${dz2} ${dz3}' ix = '${num_radial}' iy = '20 15 10' subdomain_id = '1 2 3' [] [tap3_ss] type = SideSetsBetweenSubdomainsGenerator input = 'rz_mesh' new_boundary = 'tap3' paired_block = '2' primary_block = '1' [] [tap5_ss] type = SideSetsBetweenSubdomainsGenerator input = 'tap3_ss' new_boundary = 'tap5' paired_block = '3' primary_block = '2' [] [] [Problem] coord_type = 'RZ' rz_coord_axis = 'Y' [] [GlobalParams] fp = water # name of fluid properties object rho = ${rho_fluid} # constant value for PINSFV model mu = 'mu' # material property name for viscosity porosity = porosity # name of auxvar for porosity pebble_diameter = ${pebble_diameter} velocity_interp_method = ${velocity_interp_method} advected_interp_method = ${advected_interp_method} u = vel_x v = vel_y pressure = pressure # material property name for linear drag term: Darcy_name = 'Darcy_coefficient' # material property name for quadratic drag term: Forchheimer_name = 'Forchheimer_coefficient' [] # ============================================================================== # VARIABLES AND KERNELS # ============================================================================== [Variables] [vel_x] type = PINSFVSuperficialVelocityVariable initial_condition = 1e-12 [] [vel_y] type = PINSFVSuperficialVelocityVariable initial_condition = ${v_in} [] [pressure] type = INSFVPressureVariable initial_condition = 1.0E+05 [] [] [FVKernels] # Mass Equation. [mass] type = PINSFVMassAdvection variable = pressure vel = 'superficial_velocity' [] # Momentum x component equation. [vel_x_time] type = PINSFVMomentumTimeDerivative variable = vel_x [] [vel_x_advection] type = PINSFVMomentumAdvection variable = vel_x advected_quantity = 'superficial_rho_u' vel = 'superficial_velocity' [] [vel_x_viscosity] type = PINSFVMomentumDiffusion variable = vel_x [] [vel_x_pressure] type = PINSFVMomentumPressure variable = vel_x p = pressure momentum_component = 'x' [] [vel_x_friction] type = PNSFVMomentumFriction variable = vel_x momentum_component = 'x' [] # Momentum y component equation. [vel_y_time] type = PINSFVMomentumTimeDerivative variable = vel_y [] [vel_y_advection] type = PINSFVMomentumAdvection variable = vel_y advected_quantity = 'superficial_rho_v' vel = 'superficial_velocity' [] [vel_y_viscosity] type = PINSFVMomentumDiffusion variable = vel_y [] [vel_y_pressure] type = PINSFVMomentumPressure variable = vel_y p = pressure momentum_component = 'y' [] [vel_y_friction] type = PNSFVMomentumFriction variable = vel_y momentum_component = 'y' [] [] # ============================================================================== # AUXVARIABLES AND AUXKERNELS # ============================================================================== [AuxVariables] [porosity] family = MONOMIAL order = CONSTANT fv = true # initial_condition = ${eps} [] [T_fluid] family = MONOMIAL order = CONSTANT fv = true initial_condition = ${T_in} [] [T_solid] family = MONOMIAL order = CONSTANT fv = true initial_condition = ${T_in} [] [] [AuxKernels] [eps] type = FunctionAux variable = porosity function = eps_fn execute_on = 'INITIAL' [] [] # ============================================================================== # INITIAL CONDITIONS AND FUNCTIONS # ============================================================================== [Functions] [mu_func] # provides large initial value to help start calculation type = PiecewiseLinear x = '1 3 5 10' y = '1e3 1e2 1e1 1' [] [eps_fn] type = ExponentialPorosityFunction profile = radial infinite_porosity = ${eps} pebble_diameter = ${pebble_diameter} model = duToit N = 6 wall_distance = wall_dist # name of wall distance UO [] [] # ============================================================================== # BOUNDARY CONDITIONS # ============================================================================== [FVBCs] [inlet_vel_x] type = INSFVInletVelocityBC variable = vel_x function = 1e-12 boundary = 'bottom' [] [inlet_vel_y] type = INSFVInletVelocityBC variable = vel_y function = ${v_in} boundary = 'bottom' [] [free_slip_wall_x] type = INSFVNaturalFreeSlipBC boundary = 'right' variable = vel_x [] [free_slip_wall_y] type = INSFVNaturalFreeSlipBC boundary = 'right' variable = vel_y [] [outlet_p] type = INSFVOutletPressureBC variable = pressure function = 1.0E+05 boundary = 'top' [] # Symmetry BCs for centerline [axis-vel_x] type = PINSFVSymmetryVelocityBC boundary = 'left' variable = vel_x momentum_component = x [] [axis-vel_y] type = PINSFVSymmetryVelocityBC boundary = 'left' variable = vel_y momentum_component = y [] [axis-p] type = INSFVSymmetryPressureBC boundary = 'left' variable = pressure [] [] # ============================================================================== # FLUID PROPERTIES, MATERIALS AND USER OBJECTS # ============================================================================== [FluidProperties] [water] # NIST values at 30 C and 1 atm. type = PHFunctionFluidProperties # allow_imperfect_jacobians = true rho = ${rho_fluid} k = 0.6155 mu = ${mu} cp = 4179.8 cv = 4117.2 [] [] [Materials] # FLUID [fluidprops] type = PronghornFluidProps mu_multiplier = mu_func [] [ins_fv] type = INSFVPrimitiveSuperficialVarMaterial superficial_vel_x = 'vel_x' superficial_vel_y = 'vel_y' T_fluid = 'T_fluid' T_solid = 'T_solid' rho_constant = ${rho_fluid} # constant value for INSFV formulation [] # closures in the pebble bed [drag] type = ErgunDragCoefficients # type = KTADragCoefficients # T_solid = T_solid [] [] [UserObjects] [wall_dist] type = WallDistanceCylindricalBed outer_radius = ${fparse bed_diam/2} top = ${bed_height} [] [] # ============================================================================== # EXECUTION PARAMETERS # ============================================================================== [Executioner] type = Transient solve_type = 'NEWTON' # petsc_options_iname = '-pc_type -ksp_gmres_restart' # petsc_options_value = 'lu 100' petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount' petsc_options_value = 'lu NONZERO 1e-10' # Iterations parameters l_max_its = 100 l_tol = 1e-7 nl_max_its = 25 nl_rel_tol = 5e-7 nl_abs_tol = 1e-6 # Automatic scaling automatic_scaling = true off_diagonals_in_auto_scaling = true # Problem time parameters dtmin = 0.1 dtmax = 2e4 end_time = 100 num_steps = 100 [TimeStepper] type = IterationAdaptiveDT dt = 0.1 cutback_factor = 0.5 growth_factor = 1.4 [] # Steady state detection. steady_state_detection = true steady_state_tolerance = 1e-8 steady_state_start_time = 10 [] # ============================================================================== # POSTPROCESSORS DEBUG AND OUTPUTS # ============================================================================== [Postprocessors] [p_3] type = SideAverageValue boundary = 'tap3' variable = pressure [] [p_5] type = SideAverageValue boundary = 'tap5' variable = pressure [] [delta_P] type = DifferencePostprocessor value1 = p_3 value2 = p_5 [] [dP_L] type = ParsedPostprocessor function = 'delta_P / dL' pp_names = 'delta_P' constant_names = 'dL' constant_expressions = ${dH} [] [f-mod] type = ParsedPostprocessor function = 'dP_L*dp^2*avg_eps^3/v/mu/(1-avg_eps)^2' pp_names = 'dP_L avg_eps' constant_names = 'dp v mu ' constant_expressions = '${pebble_diameter} ${v_in} ${mu}' [] [int_eps] type = FunctionElementIntegral function = eps_fn [] [avg_eps] type = ParsedPostprocessor function = 'int_eps / bed_vol' pp_names = 'int_eps' constant_names = 'bed_vol' constant_expressions = '${bed_vol}' [] [] [Outputs] file_base = 'test' exodus = true csv = true perf_graph = true []