cavity_1k

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Fluids >

cavity_1k

Previous pageReturn to chapter overviewNext page

{ 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