Potential flow over cylinder with FEM: Pressure field is flipped

Question

Much in the fashion of a few problems already answered on Mathematica.SE (1, 2), I am trying to solve the Laplace equation in velocity potential $\phi$ to simulate irrotational, inviscid flow over a cylinder (in Cartesian coordinates).

I find that my pressure field is flipped. i.e., at the front ($\theta$=0) the stagnation pressure must be a maximum while it is zero at the top and bottm ($\theta=\pi/2$ and $\theta=3\pi/2$). However, my stagnation pressure is maximum at the top and minimum in the front and back.

The Bernoulli equation was used to calculate the pressure field as follows: $$\frac{P_\text{stag}}{\rho g} = \frac{P_1}{\rho g} + \frac{\vec{V}_1^2}{2 g}$$ $$\vec{V}_1 = \frac{\partial \phi}{\partial y} + \frac{\partial \phi}{\partial x}$$ and $$\vec{V}_1^2 = \left(\frac{\partial \phi}{\partial y}\right)^2 + \left(\frac{\partial \phi}{\partial x}\right)^2$$

My code

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
Needs["NDSolve`FEM`"]

Lx = 40; Ly = 40; lx = 2; ly = 2; cx = Lx/2; cy = Ly/2; r = 5;
Ω = 
  RegionDifference[Rectangle[{0, 0}, {Lx, Ly}], Disk[{cx, cy}, r]];
RegionPlot[Ω]


cellmeasure = 1; iorder = 2;
sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
    NeumannValue[1., x == 0] + NeumannValue[-1, x == Lx], 
   DirichletCondition[u[0, 0] == 0, x == Lx && y == 0]}, 
  u, {x, y} ∈ Ω, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> cellmeasure}, 
      "IntegrationOrder" -> iorder}}]
{cylmesh} = InterpolatingFunctionCoordinates[sol];
cylmesh["Wireframe"]

Velocity vectors seem fine

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]]
StreamPlot[f[x, y], {x, 0, Lx}, {y, 0, Ly}, AspectRatio -> Automatic, 
 Frame -> True, 
 Epilog -> {Thickness[0.005], Line[{{0, lx}, {0, ly}}]}, 
 StreamPoints -> 50]

Pressure contours are INCORRECT

ClearAll[p];
P∞ = 0;
(*U∞= Evaluate[D[sol[x,y],y] + D[sol[x,y],x]];*)
ρ = 1; (*Air at 25 degree C*)
(*p=P∞+0.5ρ Evaluate[D[sol[x,y],y]^2 + \
D[sol[x,y],x]^2];*)
p = P∞ + 0.5 ρ Evaluate[Norm[f[x, y], 2]];
pplot = ContourPlot[p, {x, y} ∈ Ω, 
  PlotLegends -> Automatic, Mesh -> True, 
  ColorFunction -> "Temperature", Contours -> 25]

Velocity field around the cylinder is INCORRECT (non zero value at $\theta=0$)

I can only surmise that my "syntax" or manner of usage of some function is incorrect (assuming that the Bernoulli equation has been utilized correctly). What went wrong?

Answer 1

How are you calculating the velocity around the cylinder I am using

 Plot[Norm[f[20 + 5 Cos[\[Theta]], 20 + 5 Sin[ \[Theta]]]], {\[Theta], 
  0, 2 \[Pi]}]

to give

Which I think is correct.

For the pressure we need the correct form for Bernoulli.

Where you take the values of pressure at infinity as 0 but ignore the velocity at infinity. I am also unsure about your use of Norm. A dot product would do. Putting this into your solution gives

ClearAll[p];
P\[Infinity] = 0;
(*U\[Infinity]=Evaluate[D[sol[x,y],y]+D[sol[x,y],x]];*)
\[Rho] = 1;(*Air at 25 degree C*)(*p=P\[Infinity]+0.5\[Rho] \
Evaluate[D[sol[x,y],y]^2+D[sol[x,y],x]^2];*)
p[x_, y_] := 
 P\[Infinity] + 0.5 \[Rho] 1 - 0.5 \[Rho] Evaluate[f[x, y].f[x, y]];
pplot = ContourPlot[p[x, y], {x, y} \[Element] \[CapitalOmega], 
  PlotLegends -> Automatic, Mesh -> True, PlotRange -> All, 
  ColorFunction -> "Temperature", 
  Contours -> Table[h, {h, -1.6, 0.6, 0.05}]]

So the pressure is approximately 0.5 at the stagnation point and -1.5 at the top and bottom. Hope that helps.