Fluid Mechanics

{  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.

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 open-channel 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 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.

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^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)/(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/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 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/2-3*rball,0, 3*rball,3*rball) as "Vortex Streamlines"
  vector(u,v) as "flow"   zoom(Lx/2-3*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/2-3*rball,0, 3*rball,3*rball) as "Vortex Streamlines"
  vector(u,v) as "flow"   zoom(Lx/2-3*rball,0, 3*rball,3*rball) norm

  Transfer(u,v,p)   ! write flow solution as initial values for Coupled_Contaminant.pde

end    

 

---

{ ============= Comment ============

{  VISCOUS.PDE

This example shows the application of FlexPDE to problems in viscous flow.

The Navier-Stokes 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*(dens-dens0)
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 near-incompressibility 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 cross-section. 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 X-velocity, at the inlet, center and outlet of the channel.  The integrals presented on these plots show the consistency of mass transport across the channel.

 }

==================End Comment  ========== }