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

steady-state solutions of Navier-Stokes equations with FlexPDE.

We have included four elevation plots of X-velocity, 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.

************************************************************** }

title 'Viscous flow in 2D channel, Re > 40'

select errlim = 0.005

variables

  u(0.1)

  v(0.01)

  p(1)

definitions

  Lx = 5       Ly = 1.5

  Gx = 0       Gy = 0

  p0 = 2

  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

  cut = 0.05        { bevel the corners of the obstruction }

  penalty = 100*visc/rball^2

  Re = globalmax(speed)*(Ly/2)/visc

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)) = penalty*(dx(u)+dy(v))

Boundaries

  region 1

     start(-Lx,0)

     load(u) = 0   value(v) = 0   load(p) = 0

       line to (Lx/2-rball,0)

     value(u)=0  value(v)=0  load(p)= 0

       line to (Lx/2-rball,rball) bevel(cut)

            to (Lx/2+rball,rball) bevel(cut)

            to (Lx/2+rball,0)

     load(u) = 0  value(v) = 0  load(p) = 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

  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)

end