# How to add solutions from potential flow calculations (including circulation)?

Now we have the finite element method how do you combine potential functions appropriately? A classic problem is to calculate the lift on an aerofoil by adding circulation to a potential flow calculation. Here is a standard aerofoil in the NACA series with help from here.

ClearAll[myNACA];
myNACA[{m_, p_, t_}, x_] := Module[{},
yc = Piecewise[{{m/p^2 (2 p x - x^2),
0 <= x < p}, {m/(1 - p)^2 ((1 - 2 p) + 2 p x - x^2),
p <= x <= 1}}];
yt = 5 t (0.2969 Sqrt[x] - 0.1260 x - 0.3516 x^2 + 0.2843 x^3 -
0.1015 x^4);
θ =
ArcTan@Piecewise[{{(m*(2*p - 2*x))/p^2,
0 <= x < p}, {(m*(2*p - 2*x))/(1 - p)^2, p <= x <= 1}}];
{{x - yt Sin[θ],
yc + yt Cos[θ]}, {x + yt Sin[θ],
yc - yt Cos[θ]}}
];

m = 0.04;
p = 0.4;
tk = 0.15;
pe = myNACA[{m, p, tk}, x];
ParametricPlot[pe, {x, 0, 1}, ImageSize -> Large, Exclusions -> None] Now we make a mesh. There is a helper function and domain dimensions.

ClearAll[myLoop];
myLoop[n1_, n2_] :=
Join[Table[{n, n + 1}, {n, n1, n2 - 1, 1}], {{n2, n1}}]
Needs["NDSolveFEM"];
rt = RotationTransform[-π/16]; (* angle of attack *)
a = Table[
pe, {x, 0, 1, 0.01}]; (* table of coordinates around aerofoil *)
p0 = {p, tk/2}; (* point inside aerofoil *)
x1 = -2; x2 = 3; (* domain dimensions *)
y1 = -2; y2 = 2; (* domain dimensions *)

coords = Join[
{{x1, y1}, {x2, y1}, {x2, y2}, {x1, y2}},
rt@a[[All, 2]],
rt@Reverse[a[[All, 1]]]
];
nn = Length@coords;
bmesh = ToBoundaryMesh[
"Coordinates" -> coords,
"BoundaryElements" -> {LineElement[myLoop[1, 4]],
LineElement[myLoop[5, nn]]},
"RegionHoles" -> {rt@p0}
];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.005];
Show[mesh["Wireframe"], Frame -> True] Next we solve for the air flow assuming a potential flow around the aerofoil. There are close-up views of the streamlines.

ClearAll[x, y, ϕ];
sol = NDSolveValue[{
D[ϕ[x, y], x, x] + D[ϕ[x, y], y, y] ==
NeumannValue[1, x == x1 && y1 <= y <= y2] +
NeumannValue[-1, x == x2 && y1 <= y <= y2],
DirichletCondition[ϕ[x, y] == 0, x == 0 && y == 0]
},
ϕ, {x, y} ∈ mesh
];
ClearAll[vel];
vel[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}]]
StreamPlot[vel[x, y], {x, -0.5, 1.5}, {y, -0.5, 0.5},
PlotRange -> {{-0.5, 1.5}, {-0.5, 0.5}},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]
StreamPlot[vel[x, y], {x, 0.9, 1.05}, {y, -0.22, -0.12},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]  The airflow around the trailing edge is unrealistic -the flow will not go back like this with a separation point on the upper surface. What we need to add is circulation. So we need a second potential function that has just a circulation around the aerofoil. This is where I am going to try but am not convinced that this is the way to go.

Here is a solution where I put the circulation on the potential function on the boundaries of the region. I do this using a DirichletCondition on each of the edges.

ClearAll[x, y, ϕ];
sol1 = NDSolveValue[{
D[ϕ[x, y], x, x] + D[ϕ[x, y], y, y] == 0,
DirichletCondition[ϕ[x, y] == x - x1,
x1 <= x <= x2 && y == y1],
DirichletCondition[ϕ[x, y] == (x2 - x1) + y - y1,
x == x2 && y1 <= y <= y2],
DirichletCondition[ϕ[x,
y] == (x2 - x1) + (y2 - y1) - (x - x2),
x1 <= x <= x2 && y == y2],
DirichletCondition[ϕ[x,
y] == (x2 - x1) + (y2 - y1) + (x2 - x1) - (y - y2),
x == x1 && y1 <= y <= y2]
},
ϕ, {x, y} ∈ mesh
];
ClearAll[vel1];
vel1[x_, y_] := Evaluate[Grad[sol1[x, y], {x, y}]];
StreamPlot[vel1[x, y], {x, x1, x2}, {y, y1, y2},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]
StreamPlot[vel1[x, y], {x, -0.5, 1.5}, {y, -0.5, 0.5},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]
StreamPlot[vel1[x, y], {x, 0.9, 1.05}, {y, -0.22, -0.12},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]   Although there is circulation on the boundaries I am still getting some inflow and outflow. The large scale view shows this. How can I avoid this?

The next stage is to combine the results from the two calculations and adjust the circulation to get the separation at the trailing edge of the aerofoil. My adjusting parameter is α Here is the code and the resultant streamlines. The first one is not quite right but the second looks good.

    α = -0.3;
StreamPlot[
vel[x, y] + α vel1[x, y], {x, 0.9, 1.05}, {y, -0.22, -0.12},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]
α = -0.16;
StreamPlot[
vel[x, y] + α vel1[x, y], {x, 0.9, 1.05}, {y, -0.22, -0.12},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine

]  If I look at the potential function from the second solution I almost have what I need; it is It has a jump at x = x1, y = y1 which is what I need but the potential function is smooth thus allowing for inlets and outlets which I feel I should not have.

Thus two questions:

1. Is this a correct procedure?

2. How can I adjust the second potential function automatically?

Thanks

Here we need a realistic example to compare models of viscous and non-viscous flow. I use the same geometry, but I calculate the parameters of a viscous compressible flow with a Reynolds number = 1000 and a free-stream Mach number M = 0.25. As a solver, I use the standard finite element method and the method of integrating the Navier-Stokes equations for unsteady flows that I have developed, which I tested on several problems.

ClearAll[myNACA];
myNACA[{m_, p_, t_}, x_] :=
Module[{},
yc = Piecewise[{{m/p^2 (2 p x - x^2),
0 <= x < p}, {m/(1 - p)^2 ((1 - 2 p) + 2 p x - x^2),
p <= x <= 1}}];
yt = 5 t (0.2969 Sqrt[x] - 0.1260 x - 0.3516 x^2 + 0.2843 x^3 -
0.1015 x^4);
\[Theta] =
ArcTan@Piecewise[{{(m*(2*p - 2*x))/p^2,
0 <= x < p}, {(m*(2*p - 2*x))/(1 - p)^2, p <= x <= 1}}];
{{x - yt Sin[\[Theta]],
yc + yt Cos[\[Theta]]}, {x + yt Sin[\[Theta]],
yc - yt Cos[\[Theta]]}}];

m = 0.04;
p = 0.4;
tk = 0.15;
pe = myNACA[{m, p, tk}, x];

ClearAll[myLoop];
myLoop[n1_, n2_] :=
Join[Table[{n, n + 1}, {n, n1, n2 - 1, 1}], {{n2, n1}}]
Needs["NDSolveFEM"];
rt = RotationTransform[-\[Pi]/16];(*angle of attack*)a =
Table[pe, {x, 0, 1,
0.01}];(*table of coordinates around aerofoil*)p0 = {p,
tk/2};(*point inside aerofoil*)x1 = -2; x2 = 3;(*domain \
dimensions*)y1 = -2; y2 = 2;(*domain dimensions*)coords =
Join[{{x1, y1}, {x2, y1}, {x2, y2}, {x1, y2}}, rt@a[[All, 2]],
rt@Reverse[a[[All, 1]]]];
nn = Length@coords;
bmesh = ToBoundaryMesh["Coordinates" -> coords,
"BoundaryElements" -> {LineElement[myLoop[1, 4]],
LineElement[myLoop[5, nn]]}, "RegionHoles" -> {rt@p0}];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.001];
q = 1.4;
k = 20; Re0 = 1; U0 = 1; M0 = 0.25; Re1 = Re0/M0^2; Pin = 1;
t0 = 1/20; alpha = 0;
UX[x_, y_] := Cos[alpha];
VY[x_, y_] := Sin[alpha];
\[CapitalRho][x_, y_] := Pin;
yU = Interpolation[rt@a[[All, 1]], InterpolationOrder -> 2];
yL = Interpolation[rt@a[[All, 2]], InterpolationOrder -> 2];
Do[
{UX[i], VY[i], \[CapitalRho][i]} =
NDSolveValue[{{Inactive[
u[x, y], {x, y}]), {x, y}] +
Re1*(Abs[\[CapitalRho][i - 1][x, y]]^q)*
\!$$\*SuperscriptBox[\(\[Rho]$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]/\[CapitalRho][i - 1][x, y] +
Re0*UX[i - 1][x, y]*D[u[x, y], x] +
Re0*VY[i - 1][x, y]*D[u[x, y], y] +
Re0*(u[x, y] - UX[i - 1][x, y])/t0,
Inactive[
v[x, y], {x, y}]), {x, y}] +
Re1*(Abs[\[CapitalRho][i - 1][x, y]^q])*
\!$$\*SuperscriptBox[\(\[Rho]$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]/\[CapitalRho][i - 1][x, y] +
Re0*UX[i - 1][x, y]*D[v[x, y], x] +
Re0*VY[i - 1][x, y]*D[v[x, y], y] +
Re0*(v[x, y] - VY[i - 1][x, y])/t0,
D[\[CapitalRho][i - 1][x, y]*u[x, y], x] +
D[\[CapitalRho][i - 1][x, y]*v[x, y],
y] + (\[Rho][x, y] - \[CapitalRho][i - 1][x, y])/t0} == {0,
0, 0} /. \[Mu] -> 1/1000, {
DirichletCondition[{u[x, y] == U0*Cos[alpha],
v[x, y] == U0*Sin[alpha]}, x == x1 && y1 < y < y2],
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
y == yU[x] && 0 <= x <= Cos[Pi/16]],
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
y == yL[x] && 0 <= x <= Cos[Pi/16]],
DirichletCondition[\[Rho][x, y] == 1, x == x2]}}, {u,
v, \[Rho]}, {x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 1}}], {i, 1,
k}];


Flow velocity in magnitude and direction and pressure on the lower and upper surface of the aerofoil

ContourPlot[
Norm[{UX[i][x, y], VY[i][x, y]}] /. i -> k, {x, x1, x2}, {y, y1, y2},
PlotLegends -> Automatic, Contours -> 20,
ColorFunction -> "TemperatureMap", FrameLabel -> {"x", "y"},
PlotLabel -> V, PlotRange -> All, PlotPoints -> 50,
MaxRecursion -> 2]
sp = StreamPlot[{UX[i][x, y], VY[i][x, y]} /. i -> k, {x, -0.5,
1.5}, {y, -.5, .5}, PlotLegends -> Automatic,
ColorFunction -> "TemperatureMap", FrameLabel -> {"x", "y"},
PlotLabel -> V, PlotRange -> All, MaxRecursion -> 2,
StreamStyle -> LightGray, StreamPoints -> Fine,
AspectRatio -> Automatic];
cp = ContourPlot[
Norm[{UX[i][x, y], VY[i][x, y]}] /. i -> k, {x, -0.5,
1.5}, {y, -.5, .5}, PlotLegends -> Automatic,
ColorFunction -> "BlueGreenYellow", FrameLabel -> {"x", "y"},
PlotLabel -> V, PlotRange -> All, MaxRecursion -> 2,
Contours -> 40, AspectRatio -> Automatic, PlotPoints -> 50];

Show[cp, sp]
Plot[Evaluate[
Table[\[CapitalRho][i][x, s*(yL[x] - 10^-3) + (1 - s)*yU[x]]^(q) /.
i -> k, {s, 0, 1, 1}]], {x, 0, Cos[Pi/16]}, PlotRange -> All,
PlotPoints -> 50, MaxRecursion -> 2 ,
PlotLegends -> {"Upper Surface", "Lower Surface"},
AxesLabel -> {"x", "P"}] Streamlines on two scales (as the author chose for a potential flow). In the second figure, we see the separation and reconnection of the boundary layer.

StreamPlot[{UX[i][x, y], VY[i][x, y]} /. i -> k, {x, -0.5,
1.5}, {y, -0.5, 0.5}, PlotRange -> {{-0.5, 1.5}, {-0.5, 0.5}},
Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic,
StreamPoints -> Fine]
StreamPlot[{UX[i][x, y], VY[i][x, y]} /. i -> k, {x, 0.9,
1.05}, {y, -0.22, -0.12}, Epilog -> {Line[coords[[5 ;; nn]]]},
AspectRatio -> Automatic, StreamPoints -> Fine] • Looks like excellent work. Thank you. Is your solver code in the link?
– Hugh
Nov 12, 2018 at 9:49
• For this task, all the code is here. For transonic and supersonic flows, the code is in the debugging stage. Nov 12, 2018 at 9:59
• Thanks, your code works with minor errors in calling an interpolation function outside range. Takes awhile. I will need to study to see how it works. Good work.
– Hugh
Nov 12, 2018 at 15:04
• Thank you, do not forget that a non-stationary problem is being solved, therefore, when a flow descends from the wing, vortices may occur, as in this problem community.wolfram.com/groups/-/m/t/1433064 Nov 12, 2018 at 15:13
• Your solution is good FE work. Industry has used the method that I have begun in my question. The vortices descending from the wing are modelled as a cut in the potential flow field. There should be a jump across this cut in the potential flow calculation. This was the bit I was after in my original question. You solution is going in another direction but a good one.
– Hugh
Nov 12, 2018 at 15:17