Using ParametricNDSolveValue with parameter-dependent region to solve for a flow field around an airfoil via FEM

Question

I am looking for a way to have a parametric solution to a potential flow problem using FEM around an airfoil (as an initial bc for the inviscid bc problem), where the parameter to be varied has direct impact on the mesh (region). I provide my so far example, which is based on this post Drag and lift force acting on an airfoil

(*parameterize airfoil*)
ClearAll[NACA2415];
ClearAll[myLoop];
myLoop[n1_, n2_] := 
 Join[Table[{n, n + 1}, {n, n1, n2 - 1, 1}], {{n2, n1}}]
Needs["NDSolve`FEM`"];
NACA2415[{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]]}}];

(*generate coordinates of airfoil in the domain, depending on the angle of attack (AOA)*)
coords[AOA_] := Flatten[{
    m = 0.02;
    pp = 0.4;
    tk = 0.15;
    p0 = {pp, tk/2};(*point inside aerofoil*)
    
    pe = NACA2415[{m, pp, tk}, x];
    (*angle of attack*)
    alpha = -AOA Pi/32;
    rt = RotationTransform[alpha];
    a = Table[
      pe, {x, 0, 1, 0.01}];(*table of coordinates around aerofoil*)
  
      x1 = -1;
    x2 = 2;(*domain dimensions*)
    y1 = -1; 
    y2 = 1(*domain dimensions*);
    Join[{{x1, y1}, {x2, y1}, {x2, y2}, {x1, y2}}, rt@a[[All, 2]], 
     rt@Reverse[a[[All, 1]]]]}, 1];

and now here the challenge, I am facing:

(* generate a parametric solution, whereas the lhs/rhs AND the region may depend on real valued parameters *)
sol = ParametricNDSolveValue[{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], Element[{x, y}, mesh[AOA]], {AOA}];

Defining a parameter dependent operators/rhs/lhs is no problem, but when the mesh is parameter-dependent, it is not working. I have not found a way to somehow "delay" the evaluation of the region mesh[AOA]. A simple parameter sweep via Do or so is not a problem. Since here one could simply generate the mesh in advance. My ultimate goal, however, is to move this example further on to the viscous flow equations (also provided here Drag and lift force acting on an airfoil) and to calculate drags/lifts via varying AOAs and (!) also the airfoil geometry (as defined by the parameters in the first section) and also velocities. So sol will provide the velocity field, which I then, in combination with the mesh, will use to solve for the viscous flow. Note that changing the BC according to the AOA may also work, but as I wrote, I also want to illustrate how airfoil change will affect the Lilienthal polars. From there I am planing to move on, exploring Mathematica's POD capabilities.

I haven't used Mathematica since my physics studies 15 years or so ago. So I also appreciate any further comments on code and runtime optimization.

Answer 1

As pointed out by Alex Trounev, I have modified the code to run within 'Module[]'. This allows to run the code locally inside the module, as so to use locally defined parameters, which can be define as functional arguments. If I am right, it is close to defining a call-by-value function (as is usual in other programming languages). So something like this

f[AOA_]:=Module[{local variables}, {some computation/code depending on AOA}];

So I have treated the desired variables and fields that come out of the computation as local variables, which I also pass to the output of f[AOA]. Here I am not sure, if there is a better way and if this is even the way of using WL for this kind of problem. I am not sure for instance, how Mathematica is handling memory and code generation. To me, this way of computation especially when it comes to optimization (Gradient, Hessian), seems way over-headed (but I am not very familiar with Mathematica in this respect).

Here is the current code I am using, to sweep through the angle of attack. I noticed, that the sign of angle of attack might be messed up. (I still have to and will completely follow the code from Alex. For so far I have focussed on setting up this small test problem.):

(*parameterize airfoil*)
ClearAll[NACA2415];
ClearAll[myLoop];
myLoop[n1_, n2_] := 
 Join[Table[{n, n + 1}, {n, n1, n2 - 1, 1}], {{n2, n1}}]
Needs["NDSolve`FEM`"];
NACA2415[{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]]}}];

Now here the use of Module[] (which is also used for defining the airfoil, see Naca2415 fct.). As you can see, I define my parameter analysis function Sol[] to pass several results from the Module[] when calling Sol[]. However, post-processing is quite tedious on a multidimensional list. Here, you may know better ways to deal with it - I am thinking of keywords like in Python or other ways to call variables/lists/... by name.

Sol[aoa_] := 
  Module[{xVel, yVel, pressure, fdrag, flift, coords, mesh, bmesh, 
    bmeshFoil},
   (*parameters here define the airfoil in the naca terminology m,tk,
   pp*)
   m = 0.02; pp = 0.4;
   tk = 0.15;
   pe = NACA2415[{m, pp, tk}, x];
   (*angle of attack*)
   alpha = -aoa Pi/32;
   rt = RotationTransform[alpha];
   a = Table[
     pe, {x, 0, 1, 0.01}];(*table of coordinates around aerofoil*)
   
   p0 = {pp, tk/2};(*point inside aerofoil*)
   x1 = -1;
   x2 = 2;(*domain dimensions*)
   y1 = -1; 
   y2 = 1(*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, AccuracyGoal -> 5, PrecisionGoal -> 5, 
     "MaxCellMeasure" -> 0.0005, "MaxBoundaryCellMeasure" -> 0.01];
   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 = Evaluate[Grad[sol[x, y], {x, y}]];
   (*Now we use potential flow as a boundary condition for viscous \
flow*)
   
   bcs := {DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, x == x1], 
     DirichletCondition[{u[x, y] == vel[[1]], v[x, y] == vel[[2]]}, 
      y == y1 || y == y2], 
     DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 0 <= x <= 1], 
     DirichletCondition[{p[x, y] == 1}, x == x2]};
   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}] + Derivative[1, 0][p][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}] + Derivative[0, 1][p][x, y], 
      Derivative[1, 0][u][x, y] + 
       Derivative[0, 1][v][x, y]} /. {\[Mu] -> 10^(-3), \[Rho] -> 1};
   pde = op == {0, 0, 0}; {xVel, yVel, pressure} = 
    NDSolveValue[{pde, bcs}, {u, v, p}, Element[{x, y}, mesh], 
     Method -> {"FiniteElement", 
       "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
   ydw = Interpolation[Take[coords[[5 ;; nn]], 101]]; 
   yup = Interpolation[Take[coords[[5 ;; nn]], -101]];
   force = 
    With[{umean = 1, Y2 = ydw'[x], 
      Y1 = yup'[x], \[Rho] = 1, \[Mu] = 10^-3, dux = D[xVel[x, y], x],
       duy = D[xVel[x, y], y], dvx = D[yVel[x, y], x], 
      dvy = D[yVel[x, y], y]}, 
     Function[X, Block[{x, y, nx, ny, fx, fy, p}, {x, y} = X;
       p = pressure[x, y];
       nx = 
        If[y > x Tan[alpha], -Y1/Sqrt[1 + Y1^2], Y2/Sqrt[1 + Y2^2]];
       ny = 
        If[y > x Tan[alpha], 1/Sqrt[1 + Y1^2], -1/Sqrt[1 + Y2^2]];
       fx = nx*p + \[Mu]*(-2*nx*dux - ny*(duy + dvx));
       fy = ny*p + \[Mu]*(-nx*(dvx + duy) - 2*ny*dvy);
       {fx, fy}]]];
   bmeshFoil = 
    ToBoundaryMesh["Coordinates" -> coords[[5 ;; nn]], 
     "BoundaryElements" -> {LineElement[
        Partition[Range[Length[coords[[5 ;; nn]]]], 2, 1, 1]]}];
   {fdrag, flift} = 
    NIntegrate[force[{x, y}], {x, y} \[Element] bmeshFoil, 
     AccuracyGoal -> 3, PrecisionGoal -> 3];
   {xVel, yVel, pressure, fdrag, flift, coords, mesh, bmesh, 
    bmeshFoil}];

To perform the parameter analysis wrt. angle of attack (AOA) I have used ParallelTable which gives me 20% speedup on my Macbook i5 (early 2015) with 2 Kernels and Mathematica 12. So here is the call of the parameter scan (takes about 10 minutes on my machine):

(*parameter scan AOA*)

data = ParallelTable [{- AOA Pi/32, Sol[AOA]}, {AOA, -1.2, 1.2, .1}]// AbsoluteTiming;

Then plotting part 1/2

(*plotting of results part I/II*)
GraphicsGrid[{{ListLinePlot[
    data[[2, All, 2, 4 ;; 5]], AxesOrigin -> 0, 
    AxesLabel -> {"drag", "lift"}], 
   ListLinePlot[{Transpose[{data[[2, All, 1]], data[[2, All, 2, 4]]}],
      Transpose[{data[[2, All, 1]], data[[2, All, 2, 5]]}]}, 
    AxesOrigin -> 0, PlotRange -> All, AxesLabel -> {"AOA", "Force"}, 
    PlotLegends -> {"drag", "lift"}]}}, ImageSize -> Large]

And plotting part 2/2

(*plotting and gif export of results part II/II*)

nn = Length[data[[2, 1, 2]][[6]]];
nSamples = Length[data[[2, All]]];
framesVel = 
  Table[Show[
    ContourPlot[
     Norm[{data[[2, i, 2]][[1]][x, y], data[[2, i, 2]][[2]][x, y]}], 
     Element[{x, y}, data[[2, i, 2]][[7]]], 
     PlotLegends -> BarLegend[{Automatic, {0, 2}}, 20], 
     ColorFunction -> ColorData[{"Rainbow", {0, 2}}], 
     ColorFunctionScaling -> False, PlotRange -> All, 
     AspectRatio -> Automatic, 
     Epilog -> {Line[data[[2, i, 2]][[6]][[5 ;; nn]]]}, 
     Contours -> 20], 
    StreamPlot[{data[[2, i, 2]][[1]][x, y], 
      data[[2, i, 2]][[2]][x, y]}, 
     Element[{x, y}, data[[2, i, 2]][[7]]], StreamStyle -> LightGray, 
     AspectRatio -> Automatic], Frame -> True, 
    FrameLabel -> {"velocity"}, ImageSize -> Medium], {i, 1, 
    nSamples}];
framesPres = 
  Table[ContourPlot[data[[2, i, 2]][[3]][x, y], 
    Element[{x, y}, data[[2, i, 2]][[7]]], 
    PlotLegends -> BarLegend[{Automatic, {0, 3}}, 20], 
    ColorFunction -> ColorData[{"Rainbow", {0, 3}}], 
    ColorFunctionScaling -> False, PlotLegends -> Automatic, 
    PlotRange -> All, AspectRatio -> Automatic, 
    Epilog -> {Line[data[[2, i, 2]][[6]][[5 ;; nn]]]}, Contours -> 20,
     ImageSize -> Medium, Frame -> True, 
    FrameLabel -> {"pressure"}], {i, 1, nSamples}];
framesDragLiftAOA = 
  Table[Show[{ListLinePlot[{Transpose[{data[[2, All, 1]], 
         data[[2, All, 2, 4]]}], 
       Transpose[{data[[2, All, 1]], data[[2, All, 2, 5]]}]}, 
      AxesOrigin -> 0, PlotRange -> All, 
      AxesLabel -> {"AOA", "Force"}, PlotLegends -> {"drag", "lift"}, 
      ImageSize -> Medium, Frame -> True], 
     ListLinePlot[{{data[[2, i, 1]], 
        data[[2, i, 2, 5]]}, {data[[2, i, 1]], data[[2, i, 2, 4]]}}, 
      PlotStyle -> Black, PlotMarkers -> {Automatic, 10}]}], {i, 1, 
    nSamples}];
framesPolar = 
  Table[Show[{ListLinePlot[data[[2, All, 2, 4 ;; 5]], AxesOrigin -> 0,
       PlotRange -> All, FrameLabel -> {"drag", "lift"}, 
      PlotStyle -> Black, ImageSize -> Medium, Frame -> True], 
     ListPlot[{{data[[2, i, 2, 4]], data[[2, i, 2, 5]]}}, 
      PlotStyle -> Black, PlotMarkers -> {Automatic, 10}]}], {i, 1, 
    nSamples}];
framesBinder = 
  Table[Grid[{{framesVel[[i]], 
      framesDragLiftAOA[[i]]}, {framesPres[[i]], 
      framesPolar[[i]]}}], {i, 1, nSamples}];
(*FFanim=ListAnimate[framesBinder]*)
Export["animationnaca_AOA.gif", \
framesBinder, "DisplayDurations" -> (0.25), 
 "AnimationRepetitions" -> \[Infinity]]

So to sum up: Using Module[] allows me to parameterize the flow solving process, including mesh generation. This is nice, since it allows to analyze several things, including

• mesh dependencies
• airfoil geometry
• Reynold's number
• boundary effects etc.

for the stationary 2D case. It also opens up (at least for me now), ways to explore machine learning and optimization capabilities of Mathematica in the context of CFD and related applications.

Open questions and work remain:

• efficient referencing of the results that come out of Sol[] (my defined function for parameter analysis of a flow around some airfoil)

• respectively a results from a complex computation based on Module[]

• error handling, e.g. how to deal with convergency problems

• monitoring of the solution process

• validation of the polar (e.g. with some NACA database)

• validate integration

• the polar looks a little bit awkward though...

• numeric efficiency of the entire approach (e.g. is compilation somehow possible and an option?)

• demonstrate geometry optimization on an acceptable time scale

• try out some POD and ML approaches.