<< Click to Display Table of Contents >> lowvisc 

{ LOWVISC.PDE
This example is a modification of the VISCOUS.PDE problem, in which the
viscosity has been lowered to produce a Reynold's number of approximately
40. This seems to be the practical upper limit or Reynolds number for
steadystate solutions of NavierStokes equations with FlexPDE.
As the input pressure is raised, the disturbance in velocities propagates farther
down the channel. The channel must be long enough that the velocities
have returned to the openchannel values, or the P=0 boundary condition
at the outlet will be invalid and the solution will not succeed.
The problem computes half of the domain, with a reflective boundary at the bottom.
We have included four elevation plots of Xvelocity, at the inlet, channel
center, obstruction center and outlet of the channel. The integrals presented
on these plots show the consistency of mass transport across the channel.
We have added a variable psi to compute the stream function for plotting stream lines.
}
title 'Viscous flow in 2D channel, Re > 40'
variables u(0.1) v(0.01) p(1) psi
select ngrid = 40
definitions Lx = 5 Ly = 1.5 p0 = 2 
speed2 = u^2+v^2
speed = sqrt(speed2)
dens = 1
visc = 0.04 vxx = (p0/(2*visc*(2*Lx)))*(Ly^2y^2) { openchannel xvelocity }
rball=0.4
cut = 0.1 { value for bevel at the corners of the obstruction }
penalty = 100*visc/rball^2
Re = globalmax(speed)*(Ly/2)/(visc/dens)
w = zcomp(curl(u,v)) ! vorticity is the source for streamline equation
initial values
u = 0.5*vxx v = 0 p = p0*(Lx+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)) = penalty*(dx(u)+dy(v))
then
psi: div(grad(psi)) + w = 0 ! solve streamline equation separately from velocities
boundaries
region 1
start(Lx,0)
load(u) = 0 value(v) = 0 load(p) = 0 value(psi) = 0
line to (Lx/2rball,0)
value(u) = 0 value(v) = 0 load(p) = 0
mesh_spacing = rball/10 ! dense mesh to resolve obstruction
line to (Lx/2rball,rball) bevel(cut)
to (Lx/2+rball,rball) bevel(cut)
to (Lx/2+rball,0)
mesh_spacing = 10*rball ! cancel dense mesh requirement
load(u) = 0 value(v) = 0 load(p) = 0
line to (Lx,0)
load(u) = 0 value(v) = 0 value(p) = p0 natural(psi) = 0
line to (Lx,Ly)
value(u) = 0 value(v) = 0 load(p) = 0 natural(psi) = normal(v,u)
line to (Lx,Ly)
load(u) = 0 value(v) = 0 value(p) = 0 natural(psi) = 0
line to close
monitors
contour(speed) report(Re)
contour(psi) as "Streamlines"
contour(max(psi,0.003)) zoom(Lx/23*rball,0, 3*rball,3*rball) as "Vortex Streamlines"
vector(u,v) as "flow" zoom(Lx/23*rball,0, 3*rball,3*rball) norm
plots
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"
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)
contour(psi) as "Streamlines"
contour(max(psi,0.003)) zoom(Lx/23*rball,0, 3*rball,3*rball) as "Vortex Streamlines"
vector(u,v) as "flow" zoom(Lx/23*rball,0, 3*rball,3*rball) norm
Transfer(u,v,p) ! write flow solution as initial values for Coupled_Contaminant.pde
end