8
$\begingroup$

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]

Mathematica graphics

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["NDSolve`FEM`"];
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]

Mathematica graphics

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]

Mathematica graphics Mathematica graphics

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]

Mathematica graphics Mathematica graphics Mathematica graphics

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

]

Mathematica graphics Mathematica graphics

If I look at the potential function from the second solution I almost have what I need; it is

Mathematica graphics

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

$\endgroup$
4
$\begingroup$

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["NDSolve`FEM`"];
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[0][x_, y_] := Cos[alpha];
VY[0][x_, y_] := Sin[alpha];
\[CapitalRho][0][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[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             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[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             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"}]

fig1

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]

fig2

$\endgroup$
  • $\begingroup$ Looks like excellent work. Thank you. Is your solver code in the link? $\endgroup$ – Hugh Nov 12 '18 at 9:49
  • $\begingroup$ For this task, all the code is here. For transonic and supersonic flows, the code is in the debugging stage. $\endgroup$ – Alex Trounev Nov 12 '18 at 9:59
  • $\begingroup$ 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. $\endgroup$ – Hugh Nov 12 '18 at 15:04
  • $\begingroup$ 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 $\endgroup$ – Alex Trounev Nov 12 '18 at 15:13
  • $\begingroup$ 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. $\endgroup$ – Hugh Nov 12 '18 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.