<< Click to Display Table of Contents >> coupled_contaminant_initeq |
{ COUPLED_CONTAMINANT_INITEQ.PDE
This example is a modification of the example COUPLED_CONTAMINANT.PDE where
LOWVISC.PDE does not need to be run first. Instead it uses the INITIAL EQUATIONS
section to solve for the flow velocities.
}
title 'Contaminant transport in 2D channel, Re > 40'
variables
u(0.1)
v(0.01)
p(1)
c(0.01)
definitions
Lx = 5 Ly = 1.5
p0 = 2
speed = sqrt(u^2+v^2)
dens = 1
visc0 = 0.04
visc = visc0*(1+c)
vxx = -(p0/(2*visc*(2*Lx)))*(Ly^2-y^2) { open-channel x-velocity }
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)/(visc0/dens)
{ program a contaminant pulse in space and time
use SWAGE to eliminate discontinuous changes }
swagepulse(f,a,b,rise) = swage(f-a,0,1,rise)*swage(f-b,1,0,rise)
C0=2
cinput = C0*swage(y-0.4,1,0,0.08 )*swagepulse(t,0.4,1,0.08)
Kc = 0.002 { contaminant diffusivity }
Initial Values
u = 0.5*vxx v = 0 p = p0*(Lx+x)/(2*Lx)
initial equations
{ solve these steady-state equations to establish initial conditions for the time-dependent run }
u: visc0*div(grad(u)) - dx(p) = dens*(u*dx(u) + v*dy(u))
v: visc0*div(grad(v)) - dy(p) = dens*(u*dx(v) + v*dy(v))
p: div(grad(p)) = penalty*(dx(u)+dy(v))
equations
u: visc*div(grad(u)) - dx(p) = dens*dt(u) + dens*(u*dx(u) + v*dy(u))
v: visc*div(grad(v)) - dy(p) = dens*dt(v) + dens*(u*dx(v) + v*dy(v))
p: div(grad(p)) = penalty*(dx(u)+dy(v))
c: dt(c) + u*dx(c) + v*dy(c) = div(Kc*grad(c))
boundaries
region 1
start(-Lx,0)
load(u) = 0 value(v) = 0 load(p) = 0 load(c)=0
line to (Lx/2-rball,0)
value(u) = 0 value(v) = 0 load(p) = 0
mesh_spacing=rball/10 ! dense mesh to resolve obstruction
line to (Lx/2-rball,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
{ Introduce a lump of contaminant: }
value(c) = cinput
mesh_spacing=Ly/20
line to (Lx,Ly)
mesh_spacing=100
value(u)=0 value(v)=0 load(p)= 0 load(c)=0
line to(-Lx,Ly)
load(u) = 0 value(v) = 0 value(p) = 0
line to close
time 0 to 10
monitors
for cycle = 10
contour(speed) report(Re)
contour(c) range(0,1) report(Re)
elevation(cinput) from (Lx,-Ly) to (Lx,Ly) range=(0,C0)
plots
for t=0 by 0.5 to endtime
contour(c) range(0,1) report(Re)
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"
history(integral(c))
history(u) at (0,0.8) (2,0.8) (3,0.8) (4,0.8) (Lx,0)
history(v) at (0,0.8) (2,0.8) (3,0.8) (4,0.8)
end