Steady unidirectional flow on a 2-D domain

We want to examine the solution in comparison to a circular domain and a very wide channel.
Also, will show how to use “hydraulic diameter” to calculate pressure drop in a non-circular channel.

Copyright: Mark J. McCready
10/10/17

Three issues are considered in this notebook.  First, the numerical solution to steady, laminar flow in a rectangular domain is shown to be quite easy to get.  This obviates the motivation to learn the analytical separation of variables solution to the Poisson equation.

Second, flow in other shaped conduits is seen to be likewise easy to calculate numerically (if you know the tricks and lingo).

Finally, the utility of the friction factor -- Reynolds number relation is given.  That is, if you want to know the pressure drop in a pipe or channel, the standard way to get it is to use this relation, for a fluid of known physical properties and a design value of the flow rate, to get the pressure drop for a specified pipe size.  

Recall f≡((Δp D)/(2L ρ (v^2))) and Re≡ ((D v ρ)/μ)
For laminar flow in a circular pipe, we have shown that f = (16/Re).  

By analogy, for a non-circular pipe, we define an equivalent “diameter”, called the “hydraulic diameter” as

(Subscript[D, h])= ((4 cross-sectional area)/perimeter)

So for a rectangular channel that is h by w, the hydraulic diameter is 

(Subscript[D, h])= ((4 h w)/(2h + 2 w))

In the calculations below we show that a square channel behaves close to a circular pipe of the same size.  As the channel width is increased, the velocity is getting closer to the 1-D limit (which is not a surprise.)

We also see that flow in a triangular channel is as easy to get as for a square or rectangle.

Finally, the hydraulic diameter, does a pretty good job of correlating the “data” for all aspect ratios.  That is, you can say the f ≈ 16/Re for any aspect ratio channel if you use the hydraulic diameter.

Analytical solution of Poisson’s equation

Mathematica can do it

peqn = Laplacian[u[y,z], {y, z}] == gg

bc={u[y,0]==0,u[y,w]==0,u[0,z]==0,u[h,z]==0}

sol=DSolve[{peqn,bc},u[y,z],{y,z}, Assumptions->{h>0,w>0}]

Here is a more straightforward command.

DSolve[{D[u[y,z],{y,2}] + D[u[y,z],{z,2}] == gg, u[y,0]==0,u[y,w]==0,u[0,z]==0,u[h,z]==0},u[y,z],{y,z}]

However, if you try to put the axes in the middle there is a problem.

DSolve[{D[u[y,z],{y,2}] + D[u[y,z],{z,2}] == gg, u[y,-w/2]==0,u[y,w/2]==0,u[-h/2,z]==0,u[h/2,z]==0},u[y,z],{y,z}]

Let’s use the solution:

asol=u[y,z]/.sol[[1]]/.{∞->5}//Activate

The flow profile is as we expect.

Plot3D[(asol/.{h->.1,w->.3,gg->-10})//Evaluate,{y,0,.1},{z,0,.3},Exclusions->None,PlotRange->All]

asol/.{gg->-10,h->.1,w->.3,y->.05}

Plot[asol/.{gg->-10,h->.1,w->.3,z->.05},{y,0,.1}]

wss2/.w-> h/ϵ

Normal[%]

Limit[%,ϵ0->0]

N[%]

Get the average shear stress on the bottom wall

asol=u[y,z]/.sol[[1]]/.{∞->5}//Activate

tauxywall = 1/w Integrate[((D[asol,y])/.y->0),{z,0,w}]

%/.w-> h/ϵ

ExpandAll[%]

Normal[Series[%59,{ϵ,ϵ0,2}]]

Expand[%69]

Limit[%,ϵ0->0]

Series[ϵ Tanh[1/ϵ],{ϵ,0,2}]

Limit[ϵ Tanh[α/ϵ],ϵ->0, Assumptions->{α>0}]

Normal[Series[%,{ϵ,ϵ0,2}]]

Limit[%,ϵ0->0]

N[%]

Let’s look at the analytical solution for a rectangle

First trust that they have the correct value for wall shear stress on the main wall. (eq 3.4.28)

wss =-dpdx h/2 ( 1- 16 h/w Sum[(-1)^i Tanh[( 2 i +1) π w /2/h]/( 2 i+1)^3/ π^3,{i,0,nn}])

First, lets see how many terms it takes to get a good value for the wall shear stress.

wss0 = wss/.{nn->0,h-> .1, w-> 1, dpdx->10}

Table[wss/.{h-> .1, w-> 1, dpdx->10},{nn,0,30}]

Table[wss/.{h-> .5, w-> 1, dpdx->10},{nn,0,30}]

Table[wss/.{h-> 1., w-> 1, dpdx->10},{nn,0,30}]

It does not take many terms to converge the solution.

wss4 = wss/.{nn->4}

Now how does the value change with h/w?

Plot[{wss4/.{dpdx->10, h->1}},{w,.3,30},PlotRange->All,AxesLabel->{"w/h","τ-wall"}]

This does not mean that the pressure drop is going down.  The side walls are just carrying more stress!

Let’s do an expansion in h/w.

wss4/.w-> h/ϵ

ExpandAll[%]

Series[%,{ϵ,0,2}]

Normal[%]

Limit[%41,ϵ->0]

So we have to try some more

Series[ %,{ϵ,ϵ0,2}]

Normal[%]

The answer is:

Limit[%,ϵ0->0]

N[%44]

Analytical solution of Poisson’s equation

Mathematica can do it, but has issues!

peqn = Laplacian[u[y,z], {y, z}] == gg

bc={u[y,0]==0,u[y,h]==0,u[0,z]==0,u[h,z]==0}

sol=DSolve[{peqn,bc},u[y,z],{y,z}, Assumptions->{h>0}]

asol=u[y,z]/.sol[[1]]/.{∞->5}//Activate

Plot3D[(asol/.{h->.1,gg->-10})//Evaluate,{y,0,.1},{z,0,.1},Exclusions->None,PlotRange->All]

asol/.{gg->-10,h->.1,y->.05}

Plot[asol/.{gg->-10,h->.1,z->.05},{y,0,.1}]

Analytical solution of Poisson’s equation

Mathematica can do it, but has issues!

peqn = Laplacian[u[y,z], {y, z}] == gg

bc={u[y,0]==0,u[y,w]==0,u[0,z]==0,u[h,z]==0}

sol=DSolve[{peqn,bc},u[y,z],{y,z}, Assumptions->{h>0,w>0}]

asol=u[y,z]/.sol[[1]]/.{∞->5}//Activate

Plot3D[(asol/.{h->.1,w->.3,gg->-10})//Evaluate,{y,0,.1},{z,0,.3},Exclusions->None,PlotRange->All]

asol/.{gg->-10,h->.1,y->.05}

Plot[asol/.{gg->-10,h->.1,z->.05},{y,0,.1}]

Here is the standard solution from the text

vfield =-dpdx h^2/8/μ ( 1- 4 y^2/h^2) + dpdx h^2/8/μ Sum[( 32 (-1)^i Cosh[(2 i +1)π z/h] Cos[( 2 i +1) π y/h]/(( 2 i+1)^3 π^3 Cosh[( 2 i+1) π w /2/h])),{i,0,3}]

Integrate[vfield,{y, -h/2,h/2},{z,-w/2,w/2}]

vfield/.{dpdx->10,μ->.01,h->.1,z->0}

Plot[vfield/.{dpdx->-10,μ->1,h->.1,z->0,w->.1},{y,-.05,.05}]

oops

Show[%,%87]

Plot[vfield/.{dpdx->-10,μ->.01,h->1,z->.45,w->1},{y,-.5,.5}]

h/w = 1

Plot[vfield/.{dpdx->-10,μ->.01,h->1,y->0,w->1},{z,-.5,.5}]

h/w = 10

Plot[vfield/.{dpdx->-10,μ->.01,h->1,y->0,w->10},{z,-5,5}, PlotRange->All]

We can get whatever we need if the cross section is rectangular.

Numerical Solution for a Square channel

In these calculations, the density is set to 1.  The viscosity is 0.1.  The pressure drop varies.  So we could be thinking of a liquid 10 times as viscous as water with the same density as water.  (This is easily made with glycerin and water.)

ans1 =NDSolve[{  .5+ .1 (D[v[x,y],{x,2}]  + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[1,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,1},{y,0,1}]

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

How does this compare to the velocity calculated from flow in a wide channel? 
vmax = dpdx h^2/(8 μ)

.5 (1^2)/(8 .1)

So it is much slower.  (This is a critical observation.  If the channel is a square, you will make a really big mistake to use the 2-D result!

How about flow in a circular pipe
vmax = dpdx R^2/(4 μ)

.5 (1/2^2)/(4 .1)

This just about matches.  We are happy!

Square channel pressure drop flow rate relation

Now let’s get a data set for a square channel.

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[1,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,1},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,1},{y,0,1}],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a square channel.

ff1= Table[{(int[i]/1/1) 1 1/.1, 1 dpdxxx[i]/2/(int[i]/1/1)^2},{i,1,100,5}]

square =ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"aspect = 1"]

Let’s compare to f = 16/Re

sixteenslashRe =LogLogPlot[{16/re,6/re},{re,.1,100},PlotLegends->{"16/Re", "6/Re"}]

The agreement is pretty good.

Show[square, sixteenslashRe]

We see that it is close, but about 10% lower pressure drop than a circular pipe.  This might make sense as you could fit a circular of D = h inside the square channel and there would be some additional flow area.π d^2/4./d^2 = .79.  So the pressure drop should be lower.  However, this extra area is in the corners, with wall on two sides, so the effect is not as much as the 20% extra space would suggest.

π d^2/4./d^2

Now do a wide channel, aspect =10

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[10,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,10},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,10},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,10}, PlotRange->All]

We expect that the max V will be....

.1/8/.1

So it matches

We plot using the hydraulic diameter

ff10 = Table[{int[i]/10 (40/22.)/.1, 40/22 dpdxxx[i]/2/(int[i]/10)^2},{i,1,100,5}]

Here is a longer f - 16/Re data set.

ff16Re = Table[{j, 16/j},{j,.1,500,5}];

The wide channel comes in as higher pressure drop than the circular relation if the hydraulic diameter is used.

ListLogLogPlot[{ff1,ff10,ff16Re},PlotLegends->{"aspect = 1","aspect = 10","16/Re"},Joined->{False,False,True}]

So the pressure drop is about 20% high

Now do a still wider channel, aspect = 20

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[20,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,20},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,20},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,20}, PlotRange->All]

We expect that the max V will be....

.1/8/.1

So it matches

ff20 = Table[{int[i]/20 (80/42.)/.1, 80/42 dpdxxx[i]/2/(int[i]/20)^2},{i,1,100,5}]

So it looks like a wide channel will be high and an aspect ratio of 10 is enough to see the behavior.  (Of course for a wide channel, you could derive an exact relation, f = 6/Re, using the channel height as the length scale.)

PLotted with hydraulic diameter

ListLogLogPlot[{ff1,ff10,ff20, ff16Re},PlotLegends->{"aspect = 1","aspect = 10","aspect = 20","16/Re"},Joined->{False,False,False,True}]

Now do channel aspect = 3

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[3,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,3},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,3},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,3}, PlotRange->All]

ff3 = Table[{int[i]/3 (12/8.)/.1, 12/8 dpdxxx[i]/2/(int[i]/3)^2},{i,1,100,5}]

For aspect = 3, the result almost matches exactly to the circular pipe.

ListLogLogPlot[{ff1,ff3,ff10,ff20, ff16Re},PlotLegends->{"aspect = 1","aspect = 3","aspect = 10","aspect = 20","16/Re"},Joined->{False,False,False,False,True},AxesLabel->{"Re", "f"}]

Now do channel aspect = 2

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[2,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,2},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,2},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,2}, PlotRange->All]

So it matches

ff2 = Table[{int[i]/2 (8/6.)/.1, 8/6 dpdxxx[i]/2/(int[i]/2)^2},{i,1,100,5}]

ListLogLogPlot[{ff1,ff2, ff3,ff10,ff20, ff16Re},PlotLegends->{"aspect = 1","aspect = 2","aspect = 3","aspect = 10","aspect = 20","16/Re"},Joined->{False,False,False,False,False,True},AxesLabel->{"Re", "f"}]

You can see that for aspect ratio = 1, 2, 3... 5, the hydraulic diameter works well.

Let’s take a look at f = 6/Re for the wide channel.

ff6Re = Table[{j, 6/j},{j,.1,500,5}];

Aspect = 10

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[10,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,10},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,10},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,10}, PlotRange->All]

ffch10 = Table[{int[i]/10  /.1,   dpdxxx[i]/2/(int[i]/10)^2},{i,1,100,5}]

ListLogLogPlot[{ffch10,ff6Re},PlotLegends->{"aspect = 10","6/Re"},Joined->{False,True}]

Aspect = 20

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[20,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,20},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,20},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,20}, PlotRange->All]

ffch20 = Table[{int[i]/20  /.1,   dpdxxx[i]/2/(int[i]/20)^2},{i,1,100,5}]

Comparison of aspect ratio 10 and 20 with the 1-D limit

You can see that there is a little difference between 10 and 20, perhaps a few %.

ListLogLogPlot[{ffch10,ffch20,ff6Re},PlotLegends->{"aspect = 10","aspect = 20","6/Re"},Joined->{False,False,True}]

Just to make sure the horse is dead.... let’s do aspect = 1

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[1,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,1},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,1},{y,0,1}],{i,1,100,5}]

Plot3D[Evaluate[v[x,y]/.ans[1]],{y,0.,1},{x,0.,1}, PlotRange->All]

ffch1 = Table[{int[i]  /.1,   dpdxxx[i]/2/(int[i])^2},{i,1,100,5}]

Comparison of aspect ratio 10 and 20 with the 1-D limit

You can see that there is a little difference between 10 and 20, perhaps a few %.

ListLogLogPlot[{ffch1, ffch10,ffch20,ff6Re},PlotLegends->{"aspect = 1","aspect = 10","aspect = 20","6/Re"},Joined->{False,False,False,True}]

So, the hydraulic diameter will do a lot better than assuming that a square channel is 1-D!

Now look at other shapes

We will need a few new Mathematica functions.  So let’s use these and check to make sure that our answers match above for a square and then that a circle matches pipe flow.

Square channel pressure drop flow rate relation

Now let’s get a data set for a square channel.  We do the same domain, {0<=x<=1,0<=y<=1}

Ω=Rectangle[{0,0},{1,1}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

Graphics[Ω]

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
ans[i]=NDSolve[{  dpdxxx[i]+ .1 (D[v[x,y],{x,2}] + D[v[x,y],{y,2}]) ==0, v[0,y] ==0, v[1,y] ==0,v[x,0] ==0,v[x,1] ==0},v, {x,0,1},{y,0,1}];
int[i] = NIntegrate[v[x,y]/.ans[i][[1]],{x,0,1},{y,0,1}],{i,1,100,5}]

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a square channel.

ff1= Table[{int[i] /.1, dpdxxx[i]/2/int[i]^2},{i,1,100,5}];

ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"aspect = 1"]

We see the same agreement as above.

ListLogLogPlot[{ff1,ff16Re},PlotLegends->{"aspect = 1","16/Re"},Joined->{False,True}]

So we have verified that the different commands give the same answers as above.  Let’s do another check.

Circular pipe pressure drop -- flow rate relation

Ω=Disk[{0,0},{.5,.5}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

Graphics[Ω]

Area:

area= NIntegrate[1,{x,y}∈Ω]

Area[Ω]

Circumference

δΩ = Circle[{0,0},.5]

ArcLength[δΩ]

circumference= NIntegrate[1,{x,y}∈δΩ]

For our circle, the hydraulic diameter

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for circle.

ffcircle= Table[{int1[i]/area (dh)/.1, dh dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ffcircle,AxesLabel->{"Re","friction factor"} ,PlotLegends->"circle"]

Let’s compare to f = 16/Re

Looks like an exact match!

ListLogLogPlot[{ffcircle,ff16Re},PlotLegends->{"circle ","16/Re"},Joined->{False,True}]

Ellipse

Now try an ellipse: pipe pressure drop -- flow rate relation

Ω=Disk[{0,0},{2,1}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

Graphics[Ω]

Area:

area= NIntegrate[1,{x,y}∈Ω]

Area[Ω]

Circumference

δΩ = Circle[{0,0},{2,1}]

Graphics[δΩ]

circumference= NIntegrate[1,{x,y}∈δΩ]

ArcLength[δΩ]

For our circle, the hydraulic diameter

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for circle.

ffellipse= Table[{int1[i]/area (dh)/.1, dh dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ffellipse,AxesLabel->{"Re","friction factor"} ,PlotLegends->"ellipse"]

Let’s compare to f = 16/Re

Close, but not exact!

ListLogLogPlot[{ffellipse,ff16Re},PlotLegends->{"ellipse ","16/Re"},Joined->{False,True}]

Another Ellipse

Now try an ellipse: pipe pressure drop -- flow rate relation

Ω=Disk[{0,0},{4,1}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

Graphics[Ω]

Area:

area= NIntegrate[1,{x,y}∈Ω]

Area[Ω]

Circumference

δΩ = Circle[{0,0},{4,1}]

Graphics[δΩ]

circumference= NIntegrate[1,{x,y}∈δΩ]

ArcLength[δΩ]

For our circle, the hydraulic diameter

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for circle.

ffellipse= Table[{int1[i]/area (dh)/.1, dh dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ffellipse,AxesLabel->{"Re","friction factor"} ,PlotLegends->"ellipse"]

Let’s compare to f = 16/Re

Close, but not exact!

ListLogLogPlot[{ffellipse,ff16Re},PlotLegends->{"ellipse ","16/Re"},Joined->{False,True}]

Now a triangle!

Triangular channel pressure drop flow rate relation

Now let’s get a data set for a triangular channel.

Ω=Polygon[{{1,0},{0,Sqrt[3]},{-1,0}}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

δΩ = Line[{{1,0},{0,Sqrt[3]},{-1,0},{1,0}}]

Graphics[Ω]

Graphics[δΩ]

duval=NDSolveValue[{dpdxx+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];

ContourPlot[duval[x,y],{x,y}∈Ω, AspectRatio->1]

Plot3D[duval[x,y],{x,y}∈Ω]

Area:

area= Area[Polygon[{{1,0},{0,Sqrt[3]},{-1,0}}]]

Circumference

circumference=ArcLength[δΩ]

For our triangle, the hydraulic diameter

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a triangular channel.

Table[{int[i]/3 (12/8.)/.1, 12/8 dpdxxx[i]/2/(int[i]/3)^2},{i,1,100,5}]

ff1= Table[{int1[i]/area  (dh) /.1, (dh)dpdxxx[i]/2/(int1[i]/area )^2},{i,1,100,5}];

ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"aspect = 1"]

Let’s compare to f = 16/Re

ff16Re = Table[{i, 16/i},{i,1,100,5}];

Still the agreement looks nice.

ListLogLogPlot[{ff1,ff16Re},PlotLegends->{"equilateral triangle ","16/Re"},Joined->{False,True}]

Triangular channel pressure drop flow rate relation

Now let’s get a data set for a different triangular channel.

Ω=Polygon[{{1,0},{0,3},{-1,0}}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

δΩ = Line[{{1,0},{0,3},{-1,0},{1,0}}]

Graphics[Ω]

For our triangle, the hydraulic diameter

area =Area[Ω]

circumference=ArcLength[δΩ]

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a isosoles channel.

Table[{int[i]/area (dh)/.1, dh dpdxxx[i]/2/(area)^2},{i,1,100,5}]

ff1= Table[{int1[i]/area (dh) /.1, (dh)dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"aspect = 1"]

Let’s compare to f = 16/Re

ff16Re = Table[{i, 16/i},{i,1,100,5}];

Still the agreement looks nice.

ListLogLogPlot[{ff1,ff16Re},PlotLegends->{"equilateral triangle ","16/Re"},Joined->{False,True}]

Polygon channel pressure drop flow rate relation

Now let’s get a data set for a different triangular channel.

Ω=Polygon[{{1,0},{0,2},{-1,3},{-4,1},{-1,0}}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

δΩ = Line[{{1,0},{0,2},{-1,3},{-4,1},{-1,0},{1,0}}]

Graphics[Ω]

Graphics[δΩ]

duval=NDSolveValue[{3+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];

ContourPlot[duval[x,y],{x,y}∈Ω, AspectRatio->1]

Plot3D[duval[x,y],{x,y}∈Ω]

For our triangle, the hydraulic diameter

area =Area[Ω]

circumference=ArcLength[δΩ]

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a isosoles channel.

Table[{int1[i]/area (dh)/.1, dh dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}]

ff1= Table[{int1[i]/area (dh) /.1, (dh)dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"irregular polygon"]

Let’s compare to f = 16/Re

ff16Re = Table[{i, 16/i},{i,.1,500,10}];

Still the agreement looks nice.

ListLogLogPlot[{ff1,ff16Re},PlotLegends->{"irregular polygon ","16/Re"},Joined->{False,True}]

non-convex Polygon channel pressure drop flow rate relation

Now let’s get a data set for a different triangular channel.

Ω=Polygon[{{1,0},{0,2},{-1,1.5},{-.4,1},{-1,0}}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

Graphics[Ω]

δΩ = Line[{{1,0},{0,2},{-1,1.5},{-.4,1},{-1,0},{1,0}}]

Graphics[δΩ]

duval=NDSolveValue[{3+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];

ContourPlot[duval[x,y],{x,y}∈Ω, AspectRatio->1]

Plot3D[duval[x,y],{x,y}∈Ω]

For our triangle, the hydraulic diameter

area =Area[Ω]

circumference=ArcLength[δΩ]

dh = 4  area/circumference

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Table[{dpdxxx[i],int1[i]},{i,1,100,5}]

The relation is linear as it should be

ListPlot[%, AxesLabel->{"dpdx","flowrate"}]

Now calculate the friction factor, Reynolds number relation for a isosoles channel.

Table[{int1[i]/area (dh)/.1, dh dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}]

ff1= Table[{int1[i]/area (dh) /.1, (dh)dpdxxx[i]/2/(int1[i]/area)^2},{i,1,100,5}];

ListLogLogPlot[ff1,AxesLabel->{"Re","friction factor"} ,PlotLegends->"irregular polygon"]

Let’s compare to f = 16/Re

ff16Re = Table[{i, 16/i},{i,.1,500,10}];

Still the agreement looks nice.

ListLogLogPlot[{ff1,ff16Re},PlotLegends->{"irregular polygon ","16/Re"},Joined->{False,True}]

Let’s look at the analytical solution for a rectangle

First trust that they have the correct value for wall shear stress on the main wall. (eq 3.4.28)

wss =-dpdx h/2 ( 1- 16 h/w Sum[(-1)^i Tanh[( 2 i +1) π w /2/h]/( 2 i+1)^3/ π^3,{i,0,nn}])

First, lets see how many terms it takes to get a good value for the wall shear stress.

wss0 = wss/.{nn->0,h-> .1, w-> 1, dpdx->10}

Table[wss/.{h-> .1, w-> 1, dpdx->10},{nn,0,10}]

Table[wss/.{h-> .1, w-> .3, dpdx->10},{nn,0,10}]

wss2 = wss/.{nn->2}

-(1/2) dpdx h (1-((16 h (((Tanh[((π w)/(2 h))])/(π^3))-((Tanh[((3 π w)/(2 h))])/(27 (π^3)))+((Tanh[((5 π w)/(2 h))])/(125 (π^3)))))/w))

Plot[wss2/.{dpdx->10, h->.1},{w,.3,30},PlotRange->All]

Plot[wss/.{dpdx->10, h->.1,nn->30},{w,.3,10}]

Show[%,%700]

Now how does the value change with h/w?

wss2/.w-> h/ϵ

ExpandAll[%]

Series[%,{ϵ,0,2}]

Normal[%]

Series[ %,{ϵ,ϵ0,2}]

Normal[%]

Limit[%,ϵ0->0]

%84/.h->1

%/.dpdx->-1

Limit[%,ϵ0->0,Assumptions->{ϵ>0}]

N[%]

N[%]

wss/.w-> h/ϵ

ExpandAll[%/.nn->6]

Series[ %,{ϵ,ϵ0,2}]

Normal[%]

Limit[%,ϵ0->0]

N[%]

vfield =-dpdx h^2/8/μ ( 1- 4 y^2/h^2) + dpdx h^2/8/μ Sum[( 32 (-1)^i Cosh[(2 i +1)π z/h] Cos[( 2 i +1) π y/h]/(( 2 i+1)^3 π^3 Cosh[( 2 i+1) π w /2/h])),{i,0,3}]

vfield

vfield/.{dpdx->10,μ->.01,h->.1,z->0}

Plot[vfield/.{dpdx->-10,μ->.01,h->.1,z->0,w->10},{y,-.05,.05}]

Plot[vfield/.{dpdx->-10,μ->.01,h->1,z->0,w->1},{y,-.5,.5}]

Plot[vfield/.{dpdx->-10,μ->.01,h->1,z->.3,w->1},{y,-.5,.5}]

Square channel pressure drop flow rate relation

Now let’s get a data set for a square channel.

Ω=Polygon[{{-.5,-.5},{.5,-.5},{.5,.5},{-.5,.5}}];
Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

δΩ = Line[{{-.5,-.5},{.5,-.5},{.5,.5},{-.5,.5},{-.5,-.5}}]

Graphics[Ω]

Graphics[δΩ]

nvfield=NDSolve[{10+.01 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];

Plot3D[Evaluate[u[x,y]/.nvfield[[1]]],{x,y}∈Ω]

Plot3D[duval[x,y],{x,y}∈Ω]

nvfield

Plot[Evaluate[u[0,y]/.nvfield[[1]]],{y,-.5,.5}]

Show[%,%787]

Plot[Evaluate[u[.3,y]/.nvfield[[1]]],{y,-.5,.5}]

Show[%,%788]

Plot[Evaluate[u[.48,y]/.nvfield[[1]]],{y,-.5,.5}]

Plot[vfield/.{dpdx->-10,μ->.01,h->1,z->.48,w->1},{y,-.5,.5}]

Show[%,%826]

vfield =-dpdx h^2/8/μ ( 1- 4 y^2/h^2) + dpdx h^2/8/μ Sum[( 32 (-1)^i Cosh[(2 i +1)π z/h] Cos[( 2 i +1) π y/h]/(( 2 i+1)^3 π^3 Cosh[( 2 i+1) π w /2/h])),{i,0,11}]

Plot[vfield/.{dpdx->-10,μ->.01,h->1,z->.48,w->1},{y,-.5,.5}]

Show[%,%826]

How about an annulus

pressure drop  -- flow rate relation

Now let’s get a data set for an annular region.

Ω = Annulus[{0,0},{1,1.5}]

Graphics[Ω]

Subscript[Γ,D]=DirichletCondition[u[x,y]==0,True];

δΩ1 = Circle[{0,0},{1,1}]

δΩ2 = Circle[{0,0},{1.5,1.5}]

Show[ Graphics[ δΩ2],Graphics[ δΩ1]]

dpdxx=1

duval=NDSolveValue[{dpdxx+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];

ContourPlot[duval[x,y],{x,y}∈Ω, AspectRatio->1]

Plot3D[duval[x,y],{x,y}∈Ω]

Area:

area= Area[Ω]

Circumference

circumference=ArcLength[δΩ1] + ArcLength[δΩ2]

For our triangle, the hydraulic diameter

dh = 4  area/circumference

4 (π r2^2 - π r1^2)/( 2 π r1 + 2 π r2)

FullSimplify[%]

dpdxx = .1;

Do[dpdxxx[i] = dpdxx*i;
duval1[i]=NDSolveValue[{dpdxxx[i]+.1 Inactive[Laplacian][u[x,y],{x,y}]==0,Subscript[Γ,D]},u,{x,y}∈Ω];
int1[i] = NIntegrate[duval1[i][x,y],{x,y}∈Ω],{i,1,100,5}]

Note that “int”[i] is the volumetric flow rate.

Now calculate the friction factor, Reynolds number relation for a annular channel.

ff1= Table[{int1[i]/area  (dh) /.1, (dh)dpdxxx[i]/2/(int1[i]/area )^2},{i,1,100,5}];

ff1= Table[{int1[i]/area  (.5) /.1, (.5)dpdxxx[i]/2/(int1[i]/area )^2},{i,1,100,5}];

Let’s compare to f = 16/Re

ff16Re = Table[{i, 16/i},{i,.1,100,5}];

ff6Re = Table[{i, 6/i},{i,.1,100,5}];

The agreement with the wide channel is much better

ListLogLogPlot[{ff1,ff16Re,ff6Re},PlotLegends->{"annulus ","16/Re","6/Re"},Joined->{False,True,True}]