# How can we improve transonic flow visualization?

With this code we can make 2D FEM simulation of transonic flow around airfoil NACA 0012 at Mach number of 0.925. It takes about 5 minutes on the XENIA-15 laptop of 32 GB memory with processor Intel Core i7 [email protected] GHz.

ClearAll[NACA0012];
NACA0012[{m_, p_, t_}, x_] := Module[{}, yc = 0;
yt = 5 t (0.17814 Sqrt[x] - 0.0756 x - 0.21096 x^2 +
0.17057999999999998 x^3 - 0.06090 x^4);
\[Theta] = 0;
{{x - yt Sin[\[Theta]],
yc + yt Cos[\[Theta]]}, {x + yt Sin[\[Theta]],
yc - yt Cos[\[Theta]]}}];

m = 0.0;
p = 0.0;
tk = 0.2;
pe = NACA0012[{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["NDSolveFEM"];
rt = RotationTransform[-4.04*Pi/180];(*angle of attack*)a =
Table[pe, {x, 0, 1,
0.01}];(*table of coordinates around aerofoil*)p0 = {.2, \
.02};(*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.01];
Show[mesh["Wireframe"], Frame -> True]
q = .4;
k = 80; Re0 = 1; U0 = 1; M0 = .925; 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 -> {"PDEDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 1}}}] //
Quiet;, {i, 1, k}]


First it been tested with Mathematica 11.3 and code been published on this forum and on https://community.wolfram.com/groups/-/m/t/1433064?p_p_auth=t9Zs9KjO

Then it been tested with vv. 12.0, 12.1, 12.2. My question not about code itself.
Suppose we have prepared picture like this one for Mach number distribution. Since it is transonic flow, we can visualize supersonic region with using DensityPlot[] and ColorFunction->Hue, we have

    DensityPlot[
Norm[{UX[i][x, y], VY[i][x, y]}]*M0/\[CapitalRho][i][x, y]^(q/2) /.
i -> k, {x, -.5, 1.5}, {y, -1, 1}, ColorFunction -> Hue,
Frame -> False, PlotRange -> All, PlotPoints -> 100,
Epilog -> {Line[coords[[5 ;; nn]]]}, PlotLegends -> Automatic,
MaxRecursion -> 3]


Now we can make frames and animation like this gif

Export["D:\\M11.3\\NACA0012M0925M6.gif",
Table[DensityPlot[
Norm[{UX[i][x, y], VY[i][x, y]}]*
M0/\[CapitalRho][i][x, y]^(q/2), {x, x1, x2}, {y, y1, y2},
ColorFunction -> Function[{x, y, z}, Hue[x^1.4 ]], Frame -> False,
PlotRange -> All, PlotPoints -> 100,
Epilog -> {Line[coords[[5 ;; nn]]]}, MaxRecursion -> 2], {i, 4, k}]]


It looks fine except airfoil contour and boundary layer flow. There are number of black points and gradient looks like piecewise function. How can we remove these points?

Update 1. We can remove black points with options MaxRecursion -> 5 as follows

frames = Table[
DensityPlot[
Norm[{UX[i][x, y], VY[i][x, y]}]*
M0/\[CapitalRho][i][x, y]^(q/2), {x, x1, x2}, {y, y1, y2},
ColorFunction -> Function[{x, y, z}, Hue[Log[Exp[x^1.4]] + 1]],
Frame -> False, PlotRange -> All, PlotPoints -> 100,
Epilog -> {Line[coords[[5 ;; nn]]]}, MaxRecursion -> 5,
ImageSize -> 300], {i, 4, k, 2}];


Animation looks better, but the boundary layer gradient not resolve well

Export["D:\\...\\NACA0012M0925FR5.2.gif", frames]


Update 2. After numerical experimenting with different color data we conclude, that we can improve visualization of transonic and supersonic flow with ColorFunction -> "StarryNightColors". Here shown some example for NACA2415 at Mach number of 1.25 and angle of attack of $$\pi/16$$

This picture can be improved if we use pressure distribution instead of Mach number, then we can recognize shock and rarefaction waves

• I think the end of your first code block {i, 41, k} should read {i, 1, k}. Commented Jan 21, 2021 at 14:20
• @TimLaska Thank you, it is corrected now. Commented Jan 21, 2021 at 15:47