|
<< Click to Display Table of Contents >> cavity_1k |
  
|
|
<< Click to Display Table of Contents >> 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))
p: div(grad(p)) = penalty*(dx(u)+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