coupled_contaminant_initeq

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Fluids >

coupled_contaminant_initeq

Previous pageReturn to chapter overviewNext page

{  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