17
$\begingroup$

I want to plot the streamlines around Joukowski Airfoil using conformal mapping of a circle solution.

I do know that there are a lot of solutions to plot the airfoil itself (for example this), but I'm having difficulties plotting the streamlines around the airfoil.

Generaly, the streamlines around a circle in 2D are described by the contours of the imaginary part of:

$$F\left(z\right)=U\left(z+\frac{a^2}{z}\right)$$

Where $a$ is the radius of the circle and $U$ is the velocity at infinite distance from it. To create the Joukowski airfiol you take a circle with center $\left(a,b\right)$ with radius $R>a$ and transform it under the following conformal map:

$$w\left(z\right)=z+\frac{1}{z}$$

And then you get something like shown in here.

My problem is to plot the streamlines (the contours of $\textrm{Im}\left[F\left(z\right)\right]$) under the conformal transformation $w$. I just can't understand how do you get the streamlines to transform.

I would appreciate any help, Thanks :)

EDIT:

To be more clear, what I want to get is something like that: Streamlines around Joukowski Airfoil. The airfoil was created by transforming with the circle described in the title.

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
29
$\begingroup$

There is no simple method to display streamlines using the Zhukovsky function. I will show an example with a Zhukovsky profile at an angle of attack. 

Clear[z]; U = rho = 1; chord = 4; thk = 0.5; alpha = 
 Pi/9; y0 = 0.2; x0 = -thk/5.2; L = chord/4; a = 
 Sqrt[y0^2 + L^2]; gamma = 4 Pi a U Sin[alpha + ArcCos[L/a]];
w[z_, sign_] := 
  Module[{zeta = (z + sign Sqrt[z^2 - 4 L^2])/2}, 
   zeta = (zeta - x0 - I y0) Exp[-I alpha]/Sqrt[(1 - x0)^2 + y0^2]; 
   U (zeta + a^2/zeta) + I gamma Log[zeta]/(2 Pi)];
sign[z_] := 
 Sign[Re[z]] If[
   Abs[Re[z]] < chord/2 && 
    0 < Im[z] < 2 y0 (1 - (2 Re[z]/chord)^2), -1, 1]; 
w[z_] := w[z, sign[z]]; V[z_] = D[w[z, sig], z] /. sig -> sign[z];
bg1 = ContourPlot[Im[w[(x + I y)]], {x, -5, 5}, {y, -3, 3}, 
  AspectRatio -> Automatic, ColorFunction -> "Rainbow", 
  Contours -> Table[x^3 + 0.0208, {x, -2, 2, 0.05}], 
  ContourStyle -> White, PlotPoints -> 40, Frame -> False, 
  Exclusions -> {Log[x + I y] == 0}]

J = Show[bg1, 
  StreamPlot[{Re[V[x + I y]], -Im[V[x + I y]]}, {x, -5, 5}, {y, -3, 
    3}, AspectRatio -> Automatic, StreamStyle -> LightGray, 
   Frame -> False, StreamPoints -> Fine]]

Figure 1

We can add animation for particles carried away by the flow

    Clear[z]; U = rho = 1; chord = 4; thk = 0.5; alpha = 
 Pi/9; y0 = 0.2; x0 = -thk/5.2; L = chord/4; a = 
 Sqrt[y0^2 + L^2]; gamma = 4 Pi a U Sin[alpha + ArcCos[L/a]];
w[z_, sign_] := 
  Module[{zeta = (z + sign Sqrt[z^2 - 4 L^2])/2}, 
   zeta = (zeta - x0 - I y0) Exp[-I alpha]/Sqrt[(1 - x0)^2 + y0^2]; 
   U (zeta + a^2/zeta) + I gamma Log[zeta]/(2 Pi)];
sign[z_] := 
 Sign[Re[z]] If[
   Abs[Re[z]] < chord/2 && 
    0 < Im[z] < 2 y0 (1 - (2 Re[z]/chord)^2), -1, 1]; 
w1[z_] := w[z, sign[z]]; VX = 
 Evaluate[D[w[z, s], z] /. {z -> x + I y, s -> sign[x + I y]}];



bg = ContourPlot[Im[w1[(x + I y)]], {x, -3, 3}, {y, -2, 2}, 
  AspectRatio -> Automatic, ColorFunction -> "Rainbow", 
  Contours -> Table[x^3 + 0.0208, {x, -1.75, 1.75, 0.05}], 
  ContourStyle -> White, Exclusions -> {Log[x + I y] == 0}, 
  ClippingStyle -> Red, Frame -> False]

pX = ParametricNDSolveValue[{X'[t] == 
    Re[VX /. {x -> X[t], y -> Y[t]}], 
   Y'[t] == -Im[VX /. {x -> X[t], y -> Y[t]}], X[0] == -3, 
   Y[0] == yp}, X, {t, 0, 15}, {yp}]; pY = 
 ParametricNDSolveValue[{X'[t] == Re[VX /. {x -> X[t], y -> Y[t]}], 
   Y'[t] == -Im[VX /. {x -> X[t], y -> Y[t]}], X[0] == -3, 
   Y[0] == yp}, Y, {t, 0, 15}, {yp}];

pt = Table[
    Show[{bg, 
      Graphics[
       Table[{LightGray, PointSize[.01], 
         Point[Table[
           Evaluate[{pX[yp][t], pY[yp][t]}], {yp, -4, 1, .25}]]}, {t, 
         t0, t0 + 1, .1}]]}], {t0, 0, 12, .25}]; // Quiet

ListAnimate[pt]

Figure 2

It can be compared with the potential flow around aerodynamic profile NACA9415 (close in parameters to the Zhukovsky profile). Flow without circulation (with circulation see here)

ClearAll[NACA9415];
NACA9415[{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.09;
p = 0.4;
tk = 0.15;
pe = NACA9415[{m, p, tk}, x];
ParametricPlot[pe, {x, 0, 1}, ImageSize -> Large, Exclusions -> None]
ClearAll[myLoop];
myLoop[n1_, n2_] := 
 Join[Table[{n, n + 1}, {n, n1, n2 - 1, 1}], {{n2, n1}}]
Needs["NDSolve`FEM`"];
rt = RotationTransform[-\[Pi]/9];(*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];

ClearAll[x, y, \[Phi]];
sol = NDSolveValue[{D[\[Phi][x, y], x, x] + D[\[Phi][x, y], y, y] == 
     NeumannValue[1, x == x1 && y1 <= y <= y2] + 
      NeumannValue[-1, x == x2 && y1 <= y <= y2], 
    DirichletCondition[\[Phi][x, y] == 0, 
     x == 0 && y == 0]}, \[Phi], {x, y} \[Element] mesh];
ClearAll[vel];
vel[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}]]



st = StreamPlot[vel[x, y], {x, -.5, 1.5}, {y, -.5, .5}, 
  Epilog -> {Line[coords[[5 ;; nn]]]}, AspectRatio -> Automatic, 
  StreamPoints -> Fine, StreamStyle -> LightGray];
dp = ContourPlot[sol[x, y], {x, -.5, 1.5}, {y, -.5, .5}, 
  ColorFunction -> "Rainbow", Epilog -> {Line[coords[[5 ;; nn]]]}, 
  AspectRatio -> Automatic, Frame -> False, Contours -> 20]

bac = Show[dp, st]

Figure 3

Add animation

pX = ParametricNDSolveValue[{X'[t] == vel[X[t], Y[t]][[1]], 
   Y'[t] == vel[X[t], Y[t]][[2]], X[0] == -1/2, Y[0] == y0}, 
  X, {t, 0, 5}, {y0}]; pY = 
 ParametricNDSolveValue[{X'[t] == vel[X[t], Y[t]][[1]], 
   Y'[t] == vel[X[t], Y[t]][[2]], X[0] == -1/2, Y[0] == y0}, 
  Y, {t, 0, 5}, {y0}];

pt = Table[
    Show[{dp, 
      Graphics[
       Table[{LightGray, PointSize[.01], 
         Point[Table[
           Evaluate[{pX[y0][t], 
             pY[y0][t]}], {y0, -.5, .6, .0505}]]}, {t, t0, 
         t0 + .5, .05}]]}], {t0, 0, 2.2, .1}]; // Quiet

ListAnimate[pt]

Figure 4

Finally, using the nonlinear FEM implemented in version 12, it is possible to calculate the viscous flow for profile NACA9415. Here we see a different picture, not similar to the potential flow.

ClearAll[NACA9415];
NACA9415[{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.09;
pk = 0.4;
tk = 0.15;
pe = NACA9415[{m, pk, tk}, x];
ParametricPlot[pe, {x, 0, 1}, ImageSize -> Large, Exclusions -> None]
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 = {pk, 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.0005]; yU = 
 Interpolation[rt@a[[All, 1]], InterpolationOrder -> 2];
yL = Interpolation[rt@a[[All, 2]], 
  InterpolationOrder -> 2]; mesh["Wireframe"]
op = {Inactive[
       Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
         u[x, y], {x, y}]), {x, 
       y}] + \[Rho] {{u[x, y], v[x, y]}}.Inactive[Grad][
        u[x, y], {x, y}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y], 
    Inactive[
       Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
         v[x, y], {x, y}]), {x, 
       y}] + \[Rho] {{u[x, y], v[x, y]}}.Inactive[Grad][
        v[x, y], {x, y}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y], 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]} /. {\[Mu] -> 1/1000, \[Rho] -> 1};
pde = op == {0, 0, 0};
bcs = {DirichletCondition[u[x, y] == 1, x == x1 || y == y1 || y == y2],
   DirichletCondition[v[x, y] == 0, x == x1 || y == y1 || y == y2], 
   DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
    0 <= x <= Cos[Pi/16]], 
   DirichletCondition[p[x, y] == 0., x == x2]};
{xVel, yVel, pressure} = 
  NDSolveValue[{pde, bcs}, {u, v, p}, Element[{x, y}, mesh], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
sp = StreamPlot[{xVel[x, y], yVel[x, y]}, {x, -1, 3}, {y, -1, 1}, 
  PlotRange -> All, AspectRatio -> Automatic, StreamPoints -> Fine, 
  StreamStyle -> LightGray, Epilog -> {Line[coords[[5 ;; nn]]]}]

dp1 = DensityPlot[
  Norm[{xVel[x, y], yVel[x, y]}], {x, -1, 3}, {y, -1, 1}, 
  PlotRange -> All, AspectRatio -> Automatic, PlotPoints -> 60, 
  Frame -> False, ColorFunction -> Hue, 
  Epilog -> {Gray, Line[coords[[5 ;; nn]]]}]

Show[dp1, sp]

Figire 6

Improve this answer
$\endgroup$
5
  • 10
    $\begingroup$ I have no idea if this is correct, but +1 for pretty. :-) $\endgroup$
    – Mr.Wizard
    Dec 15, 2019 at 13:39
  • $\begingroup$ @Mr.Wizard Thank you. There is nothing to compare, but the code for w was published somewhere. I added V and removed the rotation angle. $\endgroup$
    – Alex Trounev
    Dec 15, 2019 at 13:47
  • $\begingroup$ @AlexTrounev I am trying to reproduce (on V12) the animation but bg is not defined. Using bg1 results in the error "Coordinate {ParametricFunction[1, InternalBag[{Hold[ParametricNDSolveValue[{Derivative[1][$CellContextX][$CellContextt] == Re[$CellContextV], Derivative[1][$CellContextY][$CellContextt] == -Im[$CellContextV], $CellContextX[0] == -3, $CellContext`Y[0] == $Ce should be a pair of numbers, or a Scaled or Offset form." $\endgroup$
    – Rohit Namjoshi
    Dec 15, 2019 at 17:37
  • 1
    $\begingroup$ @RohitNamjoshi I'm sorry. I added code for the animation. $\endgroup$
    – Alex Trounev
    Dec 15, 2019 at 18:32
  • $\begingroup$ @AlexTrounev Beautiful. Thanks! $\endgroup$
    – Rohit Namjoshi
    Dec 15, 2019 at 20:18
4
$\begingroup$

There are a number of related examples out there that may be of use to you such as Potential Flow over a NACA Four-Digit Airfoil and some of the other demonstrations by Richard Fearn, along with William Shaw's Complex Analysis with Mathematica book, especially the Chapter 16 (Conformal mapping I: simple mappings and Mobius transforms) and Chapter 19 (Elementary applications to two-dimensional physics).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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