  Example : Viscous Flow Script

TITLE 'Viscous flow in 2D channel, Re > 10'

SELECT errlim = 0.005

VARIABLES

u(0.1)
v(0.01)
p(1)

DEFINITIONS

Lx = 5 Ly = 1.5
p0 = 1
speed2 = u^2+v^2
speed = sqrt(speed2)
dens = 1
visc = 0.04
vxx = (p0/(2*visc*(2*Lx)))*(Ly-y)^2 { open-channel x-velocity }
rball=0.25
penalty = 100*visc/rball^2 { penalty term to enforce mass conservation }

INITIAL VALUES

u = 0.5*vxx
v = 0
p = p0*x/(2*Lx)

EQUATIONS

U : visc*div(grad(u)) - dx(p) = dens*(u*dx(u) + v*dy(u))
V : visc*div(grad(v)) - dy(p) = dens*(u*dx(v) + v*dy(v))
P : div(grad(p)) = 2*dens*(dx(U)*dy(V) - dy(U)*dx(V)) + penalty*(dx(u)+dy(v))

BOUNDARIES

REGION 1

start(-Lx,0)
line to (Lx/2-rball,0)

line to (Lx/2-rball,rball) to (Lx/2+rball,rball) to (Lx/2+rball,0)

line to (Lx,0)

load(u) = 0 value(v) = 0 value(p) = p0
line to (Lx,Ly)

value(u) = 0 value(v) = 0 load(p) = 0
line to (-Lx,Ly)

load(u) = 0 value(v) = 0 value(p) = 0
line to close

MONITORS

contour(speed)

PLOTS

grid(x,y)
contour(u)
contour(v)
contour(speed)
vector(u,v) as "flow"
contour(p) as "Pressure" painted
contour(dx(u)+dy(v)) as "Continuity Error"
elevation(u) from (-Lx,0) to (-Lx,Ly)
elevation(u) from (0,0) to (0,Ly)
elevation(u) from (Lx/2,0) to (Lx/2,Ly)
elevation(u) from (Lx,0) to (Lx,Ly)

END