Boundary Conditions

The inlet-u parameter controls the inlet velocity. Its function value can be changed to match the experimental inlet velocities. To reduce computational resources The entire assembly was split in half in the direction of flow. Since the pipe and valve is symmetric along this axis, a symmetry boundary condition was applied along this surface to approximate simulating the entire assembly. Additionally, the outlet was set as a constant pressure outlet of , and all walls were set as no-slip boundary conditions. The boundary conditions are created by the Physics/NavierStokes/Flow/flow syntax.

Variables

All simulations were done in three dimensions. Because of the turbulent nature of the simulations, a relatively slow initial velocity in the x-direction was chosen. This prevented the simulations from diverging early on. The variables are created by the Physics/NavierStokes/Flow/flow syntax.

[Physics]
  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'

        density = '${rho}'
        dynamic_viscosity = '${mu}'

        # Initial conditions
        initial_velocity = '0.1 0 0'
        initial_pressure = 0

        # Boundary conditions
        inlet_boundaries = 'Inlet'
        momentum_inlet_types = 'fixed-velocity'
        momentum_inlet_functors = '2.051481762 0 0'

        wall_boundaries = 'Walls Valve Symmetry'
        momentum_wall_types = 'noslip noslip symmetry'

        outlet_boundaries = 'Outlet'
        momentum_outlet_types = 'fixed-pressure'
        pressure_functors = '0'

        mass_advection_interpolation = '${advected_interp_method}'
        momentum_advection_interpolation = '${advected_interp_method}'
        velocity_interpolation = '${velocity_interp_method}'
      []
    []
  []
[]

AuxVariables

To implement the mixing-length turbulence model, a mixing_length auxillary variable needs to be included.

The mixing length variable is created by the Physics/NavierStokes/Turbulence/mixing_length syntax.

AuxKernels

This kernel controls the behavior of the mixing-length model. Walls 2 and 4 correspond to the wall of the pipe and the valve, while the delta parameter controls turbulent behavior near these walls.

[Physics]
  [NavierStokes]
    [Turbulence]
      [mixing_length]
        turbulence_handling = 'mixing-length'
        von_karman_const = ${von_karman_const}
        mixing_length_delta = 0.005
      []
    []
  []
[]

Executioner

While these are steady-state simulations, a transient solver must be used in order for the viscosity rampdown to function as intended. This operates by having the viscosity ramp down over time while every other parameter remains unchanged. With each successful timestep, the viscosity approaches the desired value. The simulation terminates once it ramps down all the way to the desired level.

[Executioner]
  type = Transient

  # Time-stepping parameters
  start_time = 0.0
  end_time = 8

  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 10
    # Start small, or improve the initial condition
    dt = 0.0125
    timestep_limiting_postprocessor = 'dt_limit'
  []

  # Time integration scheme
  scheme = 'implicit-euler'

  # Solver Parameters
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -snes_max_it -l_max_it -max_nl_its'
  petsc_options_value = 'lu       NONZERO                25           50        20        '
  line_search = 'none'
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-10
  automatic_scaling = true
  scaling_group_variables = 'vel_x vel_y vel_z; pressure'
[]

Viscosity Rampdown

The viscosity rampdown is controlled but two input blocks, the Functions and Materials. The functions block determines how quickly the rampdown occurs while the materials block is what defines the value of mu as show in the listing below. For these simulations, mu was first ramped down to a value of 1e-3 PaS. Once it converged to that point, it was restarted and a new rampdown function was used to quickly bring the viscosity down to 5.4e-4 PaS.

[Functions]
  [rampdown_mu_func]
    type = ParsedFunction
    expression = 'if(t<= 5, mu*(100*exp(-3*t)+1), 0.000463*exp(-3*(t-5))+0.000537)'
    symbol_names = 'mu'
    symbol_values = ${mu}
  []
[]

[FunctorMaterials]
  [mu]
    type = ADGenericFunctorMaterial #defines mu artificially for numerical convergence
    prop_names = 'mu rho' #it converges to the real mu eventually.
    prop_values = 'rampdown_mu_func ${rho}'
  []
  [velocity_vector]
    type = ADGenericVectorFunctorMaterial
    prop_names = vel_vec
    prop_values = 'vel_x vel_y vel_z'
  []
[]