The latest versions of Mathematica allows as to calculate convection even with $Ra=10^7$. Data below has been computed with the code in the attached file. It takes a time to finish computation at t=1. Code has been developed by Lion Sahara and debugged here in my answer. In the picture below shown velocity field animated with step dt=0.01. We see the turbulent like convection
.

Discussion and full code here:
Convection of air with higher Rayleigh number
https://community.wolfram.com/groups/-/m/t/3115061

Hi Alexander, Two quick questions:

1) Would it be possible to replace the cylinder by a square rod and still get a meaningful result? Or does this raise the turbulence to a level that either drives the CPU time to crazy values, or else becomes unstable?

2) Can Mathematica 12.1 also solve the unsteady compressible flow problem with temperatures in 3D (the problem with the cold cylinder and then computing several time steps, as you demonstrated further above)?

Cheers, Ron

Under the DFG Priority Research Program “Flow Simulation on High Performance Com- puters” solution methods for various flow problems have been developed

Someone above said "interesting to use old book equations and have mathematica show it" ... the above paper cited is not "old". Infact it does mostly to notate the speed of methods not show how it is solved (meaning, how to solve it).

My request is this. This is a great post for anyone to see. But "for the rest of us" it should be said at the beginning of the post if this is a heat conservation equations (or set of) or fluid flow conservation, if solved by fem instead system of pde matrix (if a choice), if this is heat flowing around a circular obstruction (what is the round thing, has it a temperature or it is a fluid obstruction) ... you know.

In other words for us beginners, please explain in english what we are physically representing / seeing and discuss if FEM is the only Mathematica solution, or if, say, monte carlo might be be easier and still show a plot (what do we get by using fem?). Explain setting up problem maybe (well, maybe a little complicated for that).

The cited .pdf file uses language which not everyone is going to know, like "laminar flow". on 2nd reading I see just above equations "fluid" is mentioned in just one sentence, easy to miss.

"12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)"

Hi Alexander, Yes, the streamflow looks like a mess in V12 as attached. Any idea why such a case and how would i be able to get the right flow?

Thank, appreciate it. GR

Attachments:

airfoil-v12.nb

Dear Alexander,

I couldn’t have done it without you.

Thanks for being such a star :-)

Great. Thank you so much, Alexander :-)

May I have the used command or notebook for my future study?

I really appreciate your help.

Hi Alexander,

I tried several times to plot the velocity distribution profile in the open channel using Piecewise command. How can I plot this profile for sections 1 and 2?

Thanks for your time.

Mathematica FEM method can be used to calculate unsteady viscous compressible flows around aerofoil in subsonic, transonic and supersonic modes. Here an example is given for the wing profile NACA2415 for the Mach number $M = 0.55$, Reynolds number $Re = 1000$, angle of attack $\alpha =\pi/16$ .

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];
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}];
p0 = {p, 
      tk/2};
x1 = -2; x2 = 3;
y1 = -2; y2 = 2;
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];
Show[mesh["Wireframe"], Frame -> True]

q = 1.4;
k = 40; Re0 = 1; U0 = 1; M0 = 0.55; 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}];

ContourPlot[
Norm[{UX[i][x, y], VY[i][x, y]}] /. i -> 32, {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 -> 32, {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 -> 32, {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]
StreamPlot[{UX[i][x, y], VY[i][x, y]} /. i -> 32, {x, 0.9, 
  1.05}, {y, -0.22, -0.12}, Epilog -> {Line[coords[[5 ;; nn]]]}, 
 AspectRatio -> Automatic, StreamPoints -> Fine]

Very nice example, thank you Alexander. I'd like to combine some features of the two models and simulate forced convection in air onto a cylinder in 2D (as in the above examples in this post), but with a temperature difference between cylinder and air. Optimally, there should be no walls (or the simulation domain should be quite large, if "no walls" boundary conditions are unstable). I'm interested in the air temperature distribution in the wake, i.e., after the air has exchanged heat with the cylinder, given a certain convective surface conduction value (unit W/m^2/K) on the cylinder wall. The Reynolds number should be rather high, but below 5 x 10^5, i.e., before the onset of full turbulence. Cheers, Ron

I consulted a textbook on boundary layer flow which states that the turbulence onset on the surface is at Re = 5 x 10^5. It may well be that turbulence in the wake starts earlier.

The cylinder axis is vertical to the convection direction, as in your first example. The cylinder has a diameter of, say, 60 cm, and is a few Kelvin cooler than the incoming air. I'd like to model the resulting air cooling pattern in space and time.

I'm not sure which approximation to choose, I'm not an expert on CFD -- that's why I'm asking in this forum.

Thanks for the 2D flow simulation around the cooled cylinder, Alexander.

Would you be able to send me the Mathematica code that you used?

Thanks, Ron

Using the nonlinear solver would be great indeed, Alexander! Cheers, Ron

After the release of Mathematica version 12, I tested a non-linear FEM on the convection problem and compared it with my method. Surprisingly, the coincidence is one to one. I congratulate the team Wolfram with a great achievement. The following code generates the same output as my code implementing the method of the false transient from the paper Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.

Pr0 = .72; Ra = 10^5; R = Ra*Pr0; t0 = 1/100; a = 0;
\[CapitalOmega] = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];

{UX, VY, \[CapitalRho], TK, TX} = 
  NDSolveValue[{{Inactive[
          Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
            u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
        u[x, y]*D[u[x, y], x] + v[x, y]*D[u[x, y], y] - 
        R*T[x, y]*Sin[a], 
       Inactive[
          Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
            v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
        u[x, y]*D[v[x, y], x] + v[x, y]*D[v[x, y], y] - 
        R*T[x, y]*Cos[a], D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 
       0} /. \[Mu] -> Pr0, {
     DirichletCondition[{p[x, y] == 0}, y == 1 && x == 1], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
      y == 0 || y == 1]}, {u[x, y]*D[T[x, y], x] + 
       v[x, y]*D[T[x, y], y] - (D[T[x, y], x, x] + 
         D[T[x, y], y, y]) == NeumannValue[0, y == 0 || y == 1], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 1}, 
      x == 0], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 0}, 
      x == 1]}}, {u, v, p, T, 
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)}, {x, y} \[Element] \[CapitalOmega], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}];

I attach two notebooks for comparison

Attachments:

Ra100000TestMMv12.nb

Ra100000TestMMv12v1.nb

Thank you, Vitaliy! I'm trying to pass two tests for the cylinder, which are described above, using nonlinear FEM. I will report the results here.

I express my gratitude to the Moderation Team for appreciating my work. I'll add an animation for the 2D2 test here and the Navier-Stokes for an unsteady flow with an explanation of the statement of the problem for tests 2D2 and 2D3 from the cited article of M. Schäfer and S. Turek

- Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!

We almost passed the test, therefore, the solution is stable. Explicit Euler is widely used in solving problems in hydrodynamics. The stability of the Euler method depends on the step of integration. In this case, the difference in solutions in v11.01 and v11.3 is due to integration over the surface of the cylinder when calculating the drag coe?cient and lift coe?cient. In our example, 2000 integrals over the surface are calculated. From these data, the maximum values of the coefficients and the Strouhal number are calculated. This is probably the main source of errors.