Oscillatory pressure driven flow

As a final problem for the semester, we want to get some insight into oscillatory flow in a tube.

The equation for pipe flow is

ρ((∂(Subscript[v, z]))/(∂t))=-((∂p)/(∂z)) + (μ/r)(∂/(∂r))(((r ∂(Subscript[v, z]))/(∂r)))

However, this is not as easy to solve with the canned solver so we will do the cross-section of flow between parallel plates.  

((∂(Subscript[v, x]))/(∂t))=-((∂p)/(∂x)) + ν(((∂^2)(Subscript[v, x]))/(∂(y^2)))

Let’s first just do the numerical solution showing the approach to steady state.

ν=.1;
dpdx=-.1;

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,5},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,5},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"}]

At t= 5, we are not at steady state.  Try a little longer...

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,10},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,10},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"}]

It is not much more than 10 however.  What if the Reynolds number were bigger?

ν=.01;
dpdx=-.1;

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,10},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,10},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"}]

Hum.  It takes longer!  (Is this consistent with boundary-layer thickening?

Ok, now look at an oscillatory flow.  We will just add an oscillatory pressure gradient.

ν=.1;
dpdx=-.1 + .05 Cos[ω t];
ω= 2 π/10;
L = 1;
Wormersly = Sqrt[ L^2 ω/ν]

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,20},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,20},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"},PlotPoints->40]

ν=.1;
dpdx=-.1 + .05 Cos[ω t];
ω= 2 π/2;
L = 1;
Wormersly = Sqrt[ L^2 ω/ν]

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,20},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,20},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"},PlotPoints->40]

Let’s look at the wall shear stress.

Plot[1/(.05)ν(v[t,0]-v[t,.05])/.ans1,{t,0,20}]

Plot[ν/(.01)(v[t,0]-v[t,.01])/.ans1,{t,0,20}]

Let’s look at just the oscillatory flow case

ν=10;
dpdx= .05 Cos[ω t];
ω= 2 π/10;
L = 1.;
Wormersly = Sqrt[ L^2 ω/ν]

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,10},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,10},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"},PlotPoints->40]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.125],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.25],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.375],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.5],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.625],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.75],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.875],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->6],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%%,%%%%%%%%]

Animate[Plot[Evaluate[(v[t,x]/.ans1)],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}],{t,1,10}]

Let’s try a faster oscilation. The Wormersly number will be higher!

ν=.01;
dpdx= .05 Cos[ω t];
ω= 2 π/5;
L = 1.;
Wormersly = Sqrt[ L^2 ω/ν]

ans1=NDSolve[{D[v[t,y],t]==ν  D[  D[v[t,y] ,y],y] - dpdx ,v[0,y]==0,v[t,1]==0,v[t,0]==0},v,{t,0,10},{y,0,1}]

Plot3D[Evaluate[v[t,x]/.ans1],{t,0,10},{x,0,1}, PlotRange->All,AxesLabel->{"t","y"},PlotPoints->40]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.125],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.25],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.375],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.5],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.625],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.75],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->5.875],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Plot[Evaluate[(v[t,x]/.ans1)/.t->6],{x,0,1}, PlotRange->All,AxesLabel->{"y","v"}]

Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%%,%%%%%%%%]

ans1

Animate[Plot[(v[t,x]/.ans1),{x,0,1},AxesLabel->{"y","v"}],{t,3,10}]