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
11/18/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.

Limit[%,ϵ0->0]

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

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}]

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 annulus, 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}]

How about an offset annulus

First the open pipe result:

Ω = ImplicitRegion[x^2 + y^2≤1,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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}∈Ω]

flowrateopen= NIntegrate[duval[x,y],{x,y}∈Ω]

Now the annulus

pressure drop  -- flow rate relation

Ω = ImplicitRegion[x^2 + y^2≤1&&!x^2 + y^2≤0.025,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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}∈Ω]

flowratecenterp025= NIntegrate[duval[x,y],{x,y}∈Ω]

ratiocenterp025= flowratecenterp025/flowrateopen

Ω = ImplicitRegion[x^2 + y^2≤1&&!(x-.05)^2 + y^2≤0.1,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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]

flowrate05p1= NIntegrate[duval[x,y],{x,y}∈Ω]

ratio05p1 = flowrate05p1/flowrateopen

Ω = ImplicitRegion[x^2 + y^2≤1&&!(x-.1)^2 + y^2≤0.1,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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]

flowratep1p1= NIntegrate[duval[x,y],{x,y}∈Ω]

ratiop1p1 = flowratep1p1/flowrateopen

Ω = ImplicitRegion[x^2 + y^2≤1&&!(x-.2)^2 + y^2≤0.1,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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]

flowratep2p1= NIntegrate[duval[x,y],{x,y}∈Ω]

ratiop2p1 = flowratep2p1/flowrateopen

try a bigger inside tube

Ω = ImplicitRegion[x^2 + y^2≤1&&!(x-0)^2 + y^2≤0.8,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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]

flowrate0p8= NIntegrate[duval[x,y],{x,y}∈Ω]

ratio0p8 = flowrate0p8/flowrateopen

Ω = ImplicitRegion[x^2 + y^2≤1&&!(x-.05)^2 + y^2≤0.8,{x,y}];

RegionPlot[Ω,AspectRatio->Automatic,ImageSize->300]

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

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]

flowratep05p8= NIntegrate[duval[x,y],{x,y}∈Ω]

ratiop05p8 = flowratep05p8/flowrateopen

Area:

area= Area[Ω]

Circumference

??Perimeter

Perimeter[Ω]

π (.4 + .6)

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

For our annulus, 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}]

laminar flow in an annulus analytical solution

DSolve[{0 == - dpdz + μ/r D[ r D[vz[r],r],r],vz[R] ==0, vz[α R] ==0}, vz[r],r]

vzr = vz[r]/.%[[1]]

vzr = FullSimplify[vz[r]/.%23[[1]], Assumptions->{dpdz ∈ Reals, R>0, α>0,μ>0}]

Plot[vzr/.{α->.5, R->1, μ->.01, dpdz->-.1},{r,.5,1}]

Plot[vzr/.{α->.9, R->1, μ->.01, dpdz->-.1},{r,.9,1}]

Plot[vzr/.{α->.2, R->1, μ->.01, dpdz->-.1},{r,.2,1}]

vzave = Integrate[vzr r,{r, α R,R}]/Integrate[ r,{r, α R,R}, Assumptions->{α>0, α<1,R>0,μ>0}]

vzave = ((dpdz (R^4) (-1+(α^2)) (1-(α^2)+(1+(α^2)) Log[α]))/(16 (((R^2)/2)-(((R^2) (α^2))/2)) μ Log[α]))

vzaverage==

Solve[vzaverage ==vzave,dpdz]

((8 vzaverage μ Log[α])/((R^2) (1-(α^2)+Log[α]+(α^2) Log[α])))   (R-α R)/ρ/vzaverage^2 /2

Limit[%,α->0]

dpdz dhh/2/ρ/V^2 == 16 μ/ρ/V/dhh

Solve[%,dpdz]

dpdz/.%[[1]]

%/.dhh-> 4 (π R^2- π (α R)^2)/(2 π R + 2 π α R)

((8 vzaverage μ Log[α])/((R^2) (1-(α^2)+Log[α]+(α^2) Log[α])))

%49/.vzaverage->V

Plot[(-(%48-%53)/%53)/.{R->1,V->1,μ->.01},{α,.0001,.99}, PlotRange->All]

LogLogPlot[(-(%48-%53)/%53)/.{R->1,V->1,μ->.01},{α,.0001,.5}]

Plot[%49/.{R->1,α->.4,μ->.01,vzaverage->V},{V,0,10}]

Show[%,%%]

Limit[ 4 (π R^2- π (α R)^2)/(2 π R + 2 π α R),α->0]

The exact limit compared to α->0

ratio = ((dpdz (R^4) (-1+(α^2)) (1-(α^2)+(1+(α^2)) Log[α]))/(16 (((R^2)/2)-(((R^2) (α^2))/2)) μ Log[α]))/Limit[((dpdz (R^4) (-1+(α^2)) (1-(α^2)+(1+(α^2)) Log[α]))/(16 (((R^2)/2)-(((R^2) (α^2))/2)) μ Log[α])),α->0]

LogLogPlot[%68/.{R->1},{α,.01,.1},PlotRange->All]