Numerical Solution of Plane Irrotational Flow

APPH 4200 Physics of Fluids
Columbia University

Solve Laplace's equation for the streamfunction, ψ[x,y], for two relatively simple examples. Since Mathematica can solve thousands of simultaneous equations quickly, the only challenge to this problem is defining the computational grid and the boundary conditions on the flow. Two dimensional solutions to elliptical PDE's can be easily and accurately solved using today's computers.

Define Laplace's Equation on a Grid

Off[General::spell1];

Define Laplace's Equation…

eqn[i_,j_] :=( ψ[i+1,j] - 2 ψ[i,j] + ψ[i-1,j])/δx^2 + (ψ[i,j+1] - 2 ψ[i,j] + ψ[i,j-1])/δy^2 == 0

eqn[i,j]

Define a Grid

height = 1.0 (* meter *);
length = 5.0 (* meter *);
wallHeight = 0.5 (* meter *);
wallLength = 0.5 (* meter *);

Grid points go from 1 to numX and from 1 to numY.

numX = 100;
numY = 80;
Print["Number of Grid Points = ", numX numY];

dx = length/(numX-1);
dy = height/(numY -1);

wallLength/dx

wallHeight/dy

flow = 5.0 (* m/s *);

grid = Table[If[length/2-wallLength/2≤ (ix-1) dx ≤ length/2 +wallLength/2 &&  (iy-1) dy≤ wallHeight || ix == 1 || ix == numX || iy == 1 || iy == numY, 1, 0], {ix, 1, numX},{iy, 1, numY}];

Dimensions[grid]

Count[Flatten[grid],x_/;x==1] (* number of known grid values *)

Count[Flatten[grid],x_/;x==0] (* number of unknown grid values *)

Length[Flatten[grid]]

ArrayPlot[Transpose[grid],Mesh->True]

wallBound[ix_] := If[length/2-wallLength/2≤ (ix-1) dx ≤ length/2 +wallLength/2, Floor[wallHeight/dy]+1, 1]

ListPlot[Table[wallBound[ix],{ix,Floor[(numX/4)],Floor[((3 numX)/4)]}],Frame->True,PlotStyle->PointSize[0.02`],Axes->False,FrameLabel->{"x","h","Wall Boundary Index",""}]

Define Boundary Conditions

Clear[ψ]; (* ψ is the variable used to store the streamfunction. *)

Top and bottom…

ψ[ix_,numY] = flow  height;
ψ[ix_,1] = 0;

Left and right sides…

ψ[1,iy_] = flow (iy-1) dy;
ψ[numX,iy_] = flow (iy-1) dy;

The wall…

Table[Table[ψ[ix,iy] = 0, {iy, 1, wallBound[ix]}], {ix, 1, numX}];

Define the Equations

eqns = Flatten[Table[Table[eqn[ix,iy],{iy, wallBound[ix]+1, numY-1}] ,{ix, 2, numX-1}]] /. δx -> dx /. δy -> dy;

eqns = Simplify[eqns];

Dimensions[eqns]

eqns[[{1,2,3}]]

unknowns = Flatten[Table[Table[ψ[ix,iy],{iy, wallBound[ix]+1, numY-1}] ,{ix, 2, numX-1}]];

unknowns[[{1,2,3}]]

Dimensions[unknowns]

Finding the Solution

Timing[sol = Flatten[Solve[eqns,unknowns]];]

ψSol = Interpolation[Flatten[Table[{(ix-1) dx, (iy-1) dy, ψ[ix,iy]}, {ix, 1, numX}, {iy, 1, numY}] /. sol,1]];

ContourPlot[ψSol[x,y],{x,0,length},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,FrameLabel->{"x","y"},PerformanceGoal -> "Speed", PlotLabel->"Streamfunction Around Wall"]

ContourPlot[ψSol[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,AspectRatio->Automatic,FrameLabel->{"x","y"},PerformanceGoal -> "Speed",PlotLabel->"Streamfunction Around Wall"]

Computing the Velocity

uX[xx_,yy_] = D[ψSol[x,y],y] /. x -> xx /. y -> yy;
uY[xx_,yy_] = -D[ψSol[x,y],x]/. x -> xx /. y -> yy;
u2[x_,y_] = uX[x,y]^2 + uY[x,y]^2;

ContourPlot[(1/2) u2[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,PerformanceGoal -> "Speed",AspectRatio->Automatic,FrameLabel->{"x","y"},PlotLabel->"U^2/2 Around Wall"]

Cp[x_,y_] = 1 - u2[x,y]/flow^2;

ContourPlot[Cp[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,PerformanceGoal -> "Speed",AspectRatio->Automatic,FrameLabel->{"x","y"},PlotLabel->"Normalized Pressure"]

Plot3D[Cp[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,PlotLabel->"Normalized Pressure"]

Show[%,ViewPoint->{1.2`,1.2`,1.2`}]

Plot[Cp[(length/2)-(height/2),y],{y,0,height},AxesLabel->{"y","Cp"},PlotLabel->"Pressure Along Leading Wall Edge"]

Plot[Cp[x,wallHeight],{x,0,length},AxesLabel->{"x","Cp"},PlotLabel->"Pressure Along Wall Height"]

Flow Past a Cylinder

Clear[ψ,uX,uY,u2];

Print["Height of channel = ",height];

radius = 0.2 (* meter *);
centerX = length/2 (* meter *);
centerY = height/2  (* meter *);

Define Boundary Conditions

Top and bottom…

ψ[ix_,numY] = flow  height;
ψ[ix_,1] = 0;

Left and right sides…

ψ[1,iy_] = flow (iy-1) dy;
ψ[numX,iy_] = flow (iy-1) dy;

The cylinder

Table[If[((ix-1) dx  - length/2)^2 + ((iy -1)dy - centerY)^2 ≤ radius^2, ψ[ix,iy] = flow centerY], {iy, 1,numY}, {ix, 1, numX}];

Define the Equations

eqns = Flatten[Table[If[((ix-1) dx  - length/2)^2 + ((iy -1)dy - centerY)^2 ≤ radius^2,-1,eqn[ix,iy]], {iy, 2,numY-1}, {ix, 2, numX-1}]] /. δx -> dx /. δy -> dy;
eqns = DeleteCases[eqns, _?Negative];

Dimensions[eqns]

eqns = Simplify[eqns];

eqns[[{1,2,3}]]

unknowns = Flatten[Table[If[((ix-1) dx  - length/2)^2 + ((iy -1)dy - centerY)^2 ≤ radius^2,-1,ψ[ix,iy]], {iy, 2,numY-1}, {ix, 2, numX-1}]];
unknowns = DeleteCases[unknowns, _?Negative];

unknowns[[{1,2,3}]]

Dimensions[unknowns]

Solution

The solution is found quickly using  Mathematica's built-in Solve function.

Timing[sol = Flatten[Solve[eqns,unknowns]];]

ψSol = Interpolation[Flatten[Table[{(ix-1) dx, (iy-1) dy, ψ[ix,iy]}, {ix, 1, numX}, {iy, 1, numY}] /. sol,1]];

ContourPlot[ψSol[x,y],{x,0,length},{y,0,height},PlotPoints->80,ContourShading->False,Contours->24,FrameLabel->{"x","y"},PerformanceGoal -> "Speed",PlotLabel->"Streamfunction Around Cylinder"]

plt=ContourPlot[ψSol[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,PerformanceGoal -> "Speed",Contours->24,FrameLabel->{"x","y"},PlotLabel->"Streamfunction Around Cylinder"]

cylinder=Graphics[{Thickness[0.02`],Circle[{centerX,centerY},radius]}]

Show[plt,cylinder]

Computing the Velocity

uX[xx_,yy_] = D[ψSol[x,y],y] /. x -> xx /. y -> yy;
uY[xx_,yy_] = -D[ψSol[x,y],x]/. x -> xx /. y -> yy;
u2[x_,y_] = uX[x,y]^2 + uY[x,y]^2;

ContourPlot[(1/2) u2[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,PerformanceGoal -> "Speed",AspectRatio->Automatic,FrameLabel->{"x","y"},PlotLabel->"U^2/2 Around Cylinder"]

Cp[x_,y_] = 1 - u2[x,y]/flow^2;

ContourPlot[Cp[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,ContourShading->False,Contours->18,PerformanceGoal -> "Speed",AspectRatio->Automatic,FrameLabel->{"x","y"},PlotLabel->"Normalized Pressure"]

Plot3D[Cp[x,y],{x,(length/2)-(height/2),(length/2)+(height/2)},{y,0,height},PlotPoints->80,PlotLabel->"Normalized Pressure"]

Show[%,ViewPoint->{1.2,1.2,1.2}]

Show[%%,PlotRange->All]

Plot[Cp[centerX,y],{y,0,height},AxesLabel->{"y","Cp"},PlotLabel->"Pressure Profile through Cylinder"]

Plot[{Cp[centerX+1.01 radius Cos[θ],centerY+1.01 radius Sin[θ]],Cp[centerX+1.1 radius Cos[θ],centerY+1.1 radius Sin[θ]]},{θ,0,2 π},AxesLabel->{"θ","Cp"},Ticks->{Range[0,2 π,(π/2)],Automatic},PlotLabel->"Pressure Profile around Cylinder"]

Drag and Lift

dDrag = P dy and dLift = P dx. We use the normalized pressure (i.e. the coefficient of drag, (Subscript[C, p]) = P/(ρ(U^2)/2)).

{totalDrag, totalLift} = {NIntegrate[Cp[centerX+1.01 radius Cos[θ],centerY+1.01 radius Sin[θ]]radius Cos[θ],{θ,0,2 π}],NIntegrate[Cp[centerX+1.01 radius Cos[θ],centerY+1.01 radius Sin[θ]]radius (-Sin[θ]),{θ,0,2 π}]}

Summary

Two dimensional potential flow requires the solution of an elliptic partial differential equation. These problems can be solved very easily with today's computers.