<< Click to Display Table of Contents >> viscous 

{ VISCOUS.PDE
This example shows the application of FlexPDE to problems in
viscous flow.
The NavierStokes equation for steady incompressible flow in two
cartesian dimensions is
dens*(dt(U) + U*dx(U) + V*dy(U)) = visc*del2(U)  dx(P) + dens*Fx
dens*(dt(V) + U*dx(V) + V*dy(V)) = visc*del2(V)  dy(P) + dens*Fy
together with the continuity equation
div[U,V] = 0
where
U and V are the X and Y components of the flow velocity
P is the fluid pressure
dens is the fluid density
visc is the fluid viscosity
Fx and Fy are the X and Y components of the body force.
In principle, the third equation enforces incompressible mass conservation,
but the equation contains no reference to the pressure, which is nominally
the remaining variable to be determined.
In practice, we choose to solve a "slightly compressible" system by defining a
hypothetical equation of state
p(dens) = p0 + L*(densdens0)
where p0 and dens0 represent a reference density and pressure, and L is a large
number representing a strong response of pressure to changes of density. L is
chosen large enough to enforce the nearincompressibility of the fluid, yet not
so large as to erase the other terms of the equation in the finite precision
of the computer arithmetic.
The compressible form of the continuity equation is
dt(dens) + div(dens*U) = 0
which, together with the equation of state yields
dt(p) = L*dens0*div(U)
In steady state, we can replace the dt(p) by div(grad(p))
[see Help  Tech Notes  Smoothing Operators in PDEs"],
resulting in the final pressure equation:
div(grad(p)) = M*div(U)
Here M has the dimensions of density/time or viscosity/distance^2.
The problem posed here shows flow in a 2D channel blocked by a bar of square
crosssection. The channel is mirrored on the bottom face, and only the upper
half is computed.
We have chosen a "convenient" value of M, one that gives good
accuracy in reasonable time. The user can alter this value to find one
which is satisfactory for his application.
We have included three elevation plots of Xvelocity, at the inlet, center
and outlet of the channel. The integrals presented on these plots show the
consistency of mass transport across the channel.
}
title 'Viscous flow in 2D channel, Re < 0.1'
variables u(0.1) v(0.01) p(1)
select ngrid = 40
definitions Lx = 5 Ly = 1.5 Gx = 0 Gy = 0 p0 = 1 speed2 = u^2+v^2 speed = sqrt(speed2) dens = 1 visc = 1 
vxx = (p0/(2*visc*Lx))*(Lyy)^2 { openchannel xvelocity }
rball=0.25
cut = 0.05 { value for bevel at the corners of the obstruction }
pfactor = staged(1,10,100,1000)
penalty = pfactor*visc/rball^2 { "equation of state" }
Re = globalmax(speed)*(Ly/2)/(visc/dens)
initial values
u = 0.5*vxx v = 0 p = p0*x/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)) = penalty*(dx(u)+dy(v))
boundaries
region 1
start(0,0)
{ unspecified boundary conditions default to LOAD=0 }
value(v)=0 line to (Lx/2rball,0)
value(u)=0 line to (Lx/2rball,rball) bevel(cut)
line to (Lx/2+rball,rball) bevel(cut)
line to (Lx/2+rball,0)
load(u)=0 line to (Lx,0)
value(p)=p0 line to (Lx,Ly)
value(u)=0 load(p)=0 line to (0,Ly)
load(u)=0 value(p)=0 line to close
monitors
contour(speed) report(Re)
plots
grid(x,y)
contour(u) report(Re)
contour(v) report(Re)
contour(speed) painted report(Re)
vector(u,v) as "flow" report(Re)
contour(p) as "Pressure" painted
contour(dx(u)+dy(v)) as "Continuity Error"
contour(p) zoom(Lx/2,0,1,1) as "Pressure"
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