INS pressure boundary condition with gravity

[Eilloo]

Hi all,

I have a question about properly setting outlet boundary conditions in the INS module when the 'INSADGravityForce' kernel is being used.

When gravity is enabled, there is naturally a pressure gradient throughout the fluid domain. Setting a constant pressure at the outlet results in back-flow, which is expected due to the pressure difference inside the fluid domain at the top and bottom of the outlet.

To avoid this problem, we can instead account for the hydrostatic pressure at the outlet by setting p = p_0 + rho*g*(h-h_ref), where p is the fixed pressure at the outlet, p_0 is the pressure you would specify if gravity wasn't included, rho is the fluid density, g is gravity, and (h-h_ref) is the distance to some reference height.

If we consider the simple case where h_ref = 0, p_0 = 0, and gravity acts in the negative y direction, this becomes p = rho*g*y. We can also initialise the pressure in the fluid domain to be rho*g*y.

However, when I do this, I see immediate pressure oscillations in space at the outlet which eventually prevent the simulation from converging. I am using PSPG stabilisation without which the oscillations are much worse.

I have also tried using a finer mesh, as well as setting the inlet velocity to zero, Neither of these changes the behaviour at the outlet.

Am I treating this boundary correctly, and is there anything else I can try to prevent these pressure oscillations?

Thanks!

(Below are two images, one showing the initialised domain, and the other showing the first time step. The inlet is a constant velocity inlet at the bottom left, and the pressure outlet is as described in the top right. I have also included the input file for reference.)

vel_x = 1
filename = './example_output'

[Mesh]
    [main]
        type = GeneratedMeshGenerator
        dim = 2
        xmin = 0
        xmax = 15
        ymin = 0
        ymax = 5
        nx = 150
        ny = 50
    []
    [block_to_remove_left]
        type = ParsedSubdomainMeshGenerator
        combinatorial_geometry = 'x<5&y>1'
        input = main
        block_name = 'left_to_delete'
        block_id = 99
    []
    [delete_left_block]
        type = BlockDeletionGenerator
        input = block_to_remove_left
        block = 'left_to_delete'
        new_boundary = 'left_cutWalls'
    []
    [block_to_remove_right]
        type = ParsedSubdomainMeshGenerator
        combinatorial_geometry = 'x>10&y<4'
        input = delete_left_block
        block_id = 98
        block_name = 'right_to_delete'
    []
    [delete_right_block]
        type = BlockDeletionGenerator
        input = block_to_remove_right
        block = 'right_to_delete'
        new_boundary = 'right_cutWalls'
    []
    [inletPipe_topSideset]
        type = ParsedGenerateSideset
        input = delete_right_block
        included_boundaries = left_cutWalls
        combinatorial_geometry = 'x<=5'
        new_sideset_name = 'inletPipe_top'
    []
    [leftWall]
        type = ParsedGenerateSideset
        input = inletPipe_topSideset
        included_boundaries = left_cutWalls
        combinatorial_geometry = 'y>=1'
        new_sideset_name = leftWall
    []
    [outletPipe_bottomSideset]
        type = ParsedGenerateSideset
        input = leftWall
        included_boundaries = right_cutWalls
        combinatorial_geometry = 'x>=10'
        new_sideset_name = 'outletPipe_bottom'
    []
    [rightWall]
        type = ParsedGenerateSideset
        input = outletPipe_bottomSideset
        included_boundaries = right_cutWalls
        combinatorial_geometry = 'y<=4'
        new_sideset_name = 'rightWall'
    []
    [inlet]
        type = ParsedGenerateSideset
        input = rightWall
        combinatorial_geometry = 'y<=1'
        included_boundaries = 'left'
        new_sideset_name = 'inlet'
    []
    [outlet]
        type = ParsedGenerateSideset
        input = inlet
        combinatorial_geometry = 'y>=4'
        included_boundaries = 'right'
        new_sideset_name = 'outlet'
    []
    [pressurePin_sideset]
        type = ExtraNodesetGenerator
        input = outlet
        coord = '0 0'
        new_boundary = 'pressurePinPoint'
    []
[]

[Variables]
    [velocity]
        family = LAGRANGE_VEC
    []
    [pressure]
    []
[]

[Functions]
    #Simple ramp function at the velocity inlet to help convergence
    rampTime = 1
    [vel_ramp]
        type = ParsedFunction
        expression = 'if(t<${rampTime},(t/${rampTime})*${vel_x},${vel_x})'
    []
[]

[ICs]
    [velocity]
        type = VectorConstantIC
        x_value = 0
        y_value = 0
        variable = velocity
    []
    [pressure]
        type = FunctionIC
        variable = pressure
        function = '-9.81*y'    #rho = 1
    []
[]

[Kernels]
    [vel_time]
        type = INSADMomentumTimeDerivative
        variable = velocity
    []
      [fluid_mass]
        type = INSADMass
        variable = pressure
    []
    [momentum_advection]
        type = INSADMomentumAdvection
        variable = velocity
    []
    [momentum_pressure]
        type = INSADMomentumPressure
        variable = velocity
        pressure = pressure
    []
    [momentum_viscous]
        type = INSADMomentumViscous
        variable = velocity
    []
    [momentum_supg]
        type = INSADMomentumSUPG
        variable = velocity
        velocity = velocity
    []
    [mass_pspg]
        type = INSADMassPSPG
        variable = pressure
    []
    [gravity]
        type = INSADGravityForce
        variable = velocity
        gravity = '0 -9.81 0'
    []
[]

[BCs]
    [inlet_vel]
        type = VectorFunctionDirichletBC
        variable = velocity
        function_x = vel_ramp
        boundary = 'inlet'
    []
    [pressureGradOutlet]                #Accounts for hydrostatic pressure, given rho = 1
        type = ADFunctionDirichletBC
        variable = pressure
        function = '-9.81*y'
        boundary = 'outlet'
    []
    [free_slip_topBottom]
        type = ADVectorFunctionDirichletBC
        variable = velocity
        boundary = 'top bottom inletPipe_top outletPipe_bottom'
        function_y = 0
        set_y_comp = true
        set_x_comp = false
    []
    [free_slip_leftRight]
        type = ADVectorFunctionDirichletBC
        variable = velocity
        boundary = 'leftWall rightWall'
        function_x = 0
        set_y_comp = false
        set_x_comp = true
    []
[]

[Materials]
    [fluid]
        type = ADGenericConstantMaterial
        prop_names = 'rho mu T cp k D alpha'
        prop_values = '1 1.789e-5 288 1005 0.026 0.01 1'
    []
    [ins_mat]
        type = INSADStabilized3Eqn
        velocity = velocity
        pressure = pressure
        temperature = '288'
        alpha = 1.0
    []
[]

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

[Executioner]
    type = Transient
    end_time = 30
    solve_type = NEWTON
    nl_rel_tol = 1e-9
    nl_max_its = 15
    l_tol = 1e-10
    l_max_its = 50

    automatic_scaling = true
    line_search = default
    
    petsc_options_iname = '-pc_type -pc_factor_shift_type'
    petsc_options_value = 'lu       NONZERO'

    dtmax = 0.1
    [TimeStepper]
        type = SolutionTimeAdaptiveDT
        dt = 0.001
    []
[]

[Outputs]
    [Exodus]
        type = Exodus
        file_base = ${filename}
    []
[]

[GiudGiud]

@lindsayad

[lindsayad]

If you want to impose the pressure condition strongly with a Dirichlet condition, then you should not integrate the pressure by parts. So in a GlobalParams block, I would put integrate_p_by_parts = false. Because right now you are getting competing conditions with your Dirichlet BC and with the natural condition at the outlet for the momentum equation which is:

p*n - n * mu \nabla u = \vec{0}

where n is the normal vector. I won't guarantee that this will fix your problem, but it is the first thing I would try

[Eilloo]

Ah, makes sense - I believe that's exactly the problem, thanks!