# Fluid Mechanics

This problem examines viscous fluid flow in a 2D channel. FlexPDE solves for the X- and Y- velocities of a fluid, with fixed pressures applied at the ends of the channel. The reynolds number in this problem is approximately 20.

The Navier-Stokes equation for steady incompressible flow in two cartesian dimensions is

rho*(dt(U) + U*dx(U) + V*dy(U)) = mu*div(grad(U)) - dx(P)

rho*(dt(V) + U*dx(V) + V*dy(V)) = mu*div(grad(V)) - dy(P)

together with the continuity equation,

dx(U) + dy(V) = 0.

Here U and V are the X- and Y- velocities, P is the pressure (introduced as a surrogate for the continuity equation), rho is density, and mu is viscosity.

## Example: Viscous Flow in a 2D channel

### Adaptively Refined Mesh

### Fluid Speed

### Pressure

{ 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 ========== }