FreeFEM Documentation on GitHub

stars - forks

Irrotational Fan Blade Flow and Thermal effects

Summary : Here we will learn how to deal with a multi-physics system of PDEs on a complex geometry, with multiple meshes within one problem. We also learn how to manipulate the region indicator and see how smooth is the projection operator from one mesh to another.

Incompressible flow

Without viscosity and vorticity incompressible flows have a velocity given by:

\[\begin{split}\boldsymbol{u}=\left(\begin{matrix}{\partial \psi \over \partial x }\\ -\partial \psi \over \partial y \end{matrix}\right), \quad \mbox{ where }\psi\mbox{ is solution of }\quad \Delta \psi =0\end{split}\]

This equation expresses both incompressibility (\(\nabla\cdot u=0\)) and absence of vortex (\(\nabla\times u =0\)).

As the fluid slips along the walls, normal velocity is zero, which means that \(\psi\) satisfies:

\[\psi \hbox{ constant on the walls}.\]

One can also prescribe the normal velocity at an artificial boundary, and this translates into non constant Dirichlet data for \(\psi\).


Let us consider a wing profile \(S\) in a uniform flow. Infinity will be represented by a large circle \(C\) where the flow is assumed to be of uniform velocity; one way to model this problem is to write:

\[\Delta \psi =0 \hbox{ in } \Omega, \qquad \psi |_{S}=-l, \quad \psi|_{C}= {\boldsymbol{u}_\infty}.\boldsymbol{x}^\perp\]

where \(\partial\Omega=C\cup S\) and \(l\) is the lift force.

The NACA0012 Airfoil

An equation for the upper surface of a NACA0012 (this is a classical wing profile in aerodynamics) is:

\[y = 0.17735\sqrt{x}-0.075597x- 0.212836x^2+0.17363x^3-0.06254x^4.\]
 1// Parameters
 2int  S = 99;// wing label
 3// u infty
 4real theta = 8*pi/180;// // 1 degree on incidence =>  lift
 5real lift = theta*0.151952/0.0872665; //  lift approximation formula
 6real  uinfty1= cos(theta), uinfty2= sin(theta);
 7// Mesh
 8func naca12 = 0.17735*sqrt(x) - 0.075597*x - 0.212836*(x^2) + 0.17363*(x^3) - 0.06254*(x^4);
 9border C(t=0., 2.*pi){x=5.*cos(t); y=5.*sin(t);}
10border Splus(t=0., 1.){x=t; y=naca12; label=S;}
11border Sminus(t=1., 0.){x=t; y=-naca12; label=S;}
12mesh Th = buildmesh(C(50) + Splus(70) + Sminus(70));
14// Fespace
15fespace Xh(Th, P2);
16Xh psi, w;
17macro grad(u) [dx(u),dy(u)]// def of grad operator
18// Solve
19solve potential(psi, w)
20   = int2d(Th)(
21        grad(psi)'*grad(w) //  scalar product
22   )
23   + on(C, psi = [uinfty1,uinfty2]'*[y,-x])
24   + on(S, psi=-lift) // to get a correct value
25   ;
27plot(psi, wait=1);

A zoom of the streamlines are shown on Fig. 15.


Fig. 15 Zoom around the NACA0012 airfoil showing the streamlines (curve \(\psi=\) constant). To obtain such a plot use the interactive graphic command: “+” and p.


Fig. 16 Temperature distribution at time T=25 (now the maximum is at 90 instead of 120).

The NACA0012 Airfoil

Heat Convection around the airfoil

Now let us assume that the airfoil is hot and that air is there to cool it. Much like in the previous section the heat equation for the temperature \(v\) is

\[\partial_t v -\nabla\cdot(\kappa\nabla v) + u\cdot\nabla v =0,~~v(t=0)=v_0, ~~\frac{\partial v}{\partial\boldsymbol{n}}|_{C,u\cdot n>0}=0, v|_{C,u\cdot n<0}=0\]

But now the domain is outside AND inside \(S\) and \(\kappa\) takes a different value in air and in steel. Furthermore there is convection of heat by the flow, hence the term \(u\cdot\nabla v\) above.

Consider the following, to be plugged at the end of the previous program:

 1// Corrected by F. Hecht may 2021
 2// Parameters
 3real S = 99;
 5border C(t=0, 2*pi){x=3*cos(t); y=3*sin(t);} // Label 1,2
 6border Splus(t=0, 1){x=t; y=0.17735*sqrt(t) - 0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4); label=S;}
 7border Sminus(t=1, 0){x=t; y=-(0.17735*sqrt(t) - 0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4)); label=S;}
 8mesh Th = buildmesh(C(50) + Splus(70) + Sminus(70));
 9// Fespace
10fespace Vh(Th, P2);
11Vh psi, w;
13// Problem
14solve potential(psi, w)
15= int2d(Th)(dx(psi)*dx(w)+dy(psi)*dy(w))
16     + on(C, psi = y)
17     + on(S, psi=0);
19// Plot
20plot(psi, wait=1);
22/// Thermic
23// Parameters
24real dt = 0.005, nbT = 50;
26// Mesh
27border D(t=0, 2){x=1+t; y=0;}
28mesh Sh = buildmesh(C(25) + Splus(-90) + Sminus(-90) + D(200));
29int steel = Sh(0.5, 0).region, air = Sh(-1, 0).region;
30// Change label to put BC on In flow
31// Fespace
32fespace Wh(Sh, P1);
33Wh  vv;
35fespace W0(Sh, P0);
36W0 k = 0.01*(region == air) + 0.1*(region == steel);
37W0 u1 = dy(psi)*(region == air), u2 = -dx(psi)*(region == air);
38Wh v = 120*(region == steel), vold;
39// set the label to 10 on inflow boundary to inforce the temperature.
40Sh = change(Sh,flabel = (label == C &&  [u1,u2]'*N<0) ? 10 : label);
41int i;
42problem thermic(v, vv, init=i, solver=LU)
43= int2d(Sh)(
44              v*vv/dt + k*(dx(v)*dx(vv) + dy(v)*dy(vv))
45            + 10*(u1*dx(v) + u2*dy(v))*vv
46           )
47- int2d(Sh)(vold*vv/dt)
48+ on(10, v= 0);
51for(i = 0; i < nbT; i++) {
52  vold[]= v[];
53  thermic;
54  plot(v);


How steel and air are identified by the mesh parameter region which is defined when buildmesh is called and takes an integer value corresponding to each connected component of \(\Omega\);


We use the change function to put label 10 on inflow boundary, remark the trick to chanhe only label C flabel = (label == C && [u1,u2]'*N<0) ? 10 : label

How the convection terms are added without upwinding. Upwinding is necessary when the Pecley number \(|u|L/\kappa\) is large (here is a typical length scale), The factor 10 in front of the convection terms is a quick way of multiplying the velocity by 10 (else it is too slow to see something).

The solver is Gauss’ LU factorization and when init \(\neq 0\) the LU decomposition is reused so it is much faster after the first iteration.

Table of content