﻿ Sample Problems > Applications > Fluids > cavity_1k

# cavity_1k

Navigation:  Sample Problems > Applications > Fluids >

# cavity_1k   { CAVITY_1K.PDE

This problem computes the flow velocities in a square cavity driven by a velocity

on the top surface, with a Reynolds number of 1000.

The initial conditions are zero interior velocity and zero pressure.

We use a non-dimensional form of the Navier-Stokes equations, with a "penalty pressure".

The pressure equation simulates a "somewhat compressible" fluid.

Because a very hard incompressibility is subject to instabilities when the

initial condition is far from equilibrium, we stage the pressure penalty,

starting with a very soft fluid and progressing to more and more strict incompressibility.

As a visual assistance, we also display the stream function and vorticity when each stage completes.

This is a common test problem, see for example Hendriana and Bathe, Int J Numer Meth Engng 47, 317-340 (2000)

}

title 'Lid-driven Cavity Re=1000'

coordinates cartesian2

select

ngrid=20

regrid=off

busymon     ! display monitors during iteration

variables

u,v

p

psi

definitions

Lx=Ly=1

Re = 1000

penalty = staged(1,10,100,1000) ! perform the computation in stages with increasing pressure penalty

w = dx(v)-dy(u)     ! report vorticity

u0 = 1

h=1/10

equations

u:  div(grad(u)) - Re*dx(p)  = Re*(u*dx(u) + v*dy(u))

v:  div(grad(v)) - Re*dy(p)  = Re*(u*dx(v) + v*dy(v))

then

psi: div(grad(psi)) +w = 0   ! Report the Stream function

boundaries

region 1

start (0,1)

! left wall

value(u)=0 value(v)=0 ! no fluid slip on bottom and sides

value(psi)=0   !  stream function BC: no flow across walls or top

line to (0,0)

! bottom

value(p)=0 !fix the pressure at the bottom

value(psi)=0

line to(1,0)

! right wall

natural(p) = 0

natural(psi)=normal(-v,u)

line to (1,1)

! open top with imposed velocity

! since a hard shift from 0 to 1 would invite catastrophe, we ramp the driving velocity to zero at the corners.

value(u) = (1-x)*u0/h

mesh_spacing=h/5 ! more mesh density at the corners helps define the rapid velocity transitions

line to(1-h,1)

! full driving velocity

value(u)=u0

value(v)=0

mesh_spacing=h

line to(h,1)

! ramp down to zero

value(u) = x*u0/h

mesh_spacing=h/5

line to close

monitors

contour(u)   report(Re) report(penalty)

contour(v)   report(Re) report(penalty)

contour(p)   report(Re) report(penalty)

plots

grid(x,y)   report(Re) report(penalty)

contour(u)   report(Re) report(penalty)

contour(v)   report(Re) report(penalty)

contour(p)   report(Re) report(penalty)

vector(u,v) report(Re) report(penalty)

contour(psi) as "Streamlines" report(Re) report(penalty)

contour(w)   as "Vorticity"   report(Re) report(penalty)

end