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\).

Airfoil

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
 2 int  S = 99;// wing label
 3 // u infty
 4 real theta = 8*pi/180;// // 1 degree on incidence =>  lift
 5 real lift = theta*0.151952/0.0872665; //  lift approximation formula
 6 real  uinfty1= cos(theta), uinfty2= sin(theta);
 7 // Mesh
 8 func naca12 = 0.17735*sqrt(x) - 0.075597*x - 0.212836*(x^2) + 0.17363*(x^3) - 0.06254*(x^4);
 9 border C(t=0., 2.*pi){x=5.*cos(t); y=5.*sin(t);}
10 border Splus(t=0., 1.){x=t; y=naca12; label=S;}
11 border Sminus(t=1., 0.){x=t; y=-naca12; label=S;}
12 mesh Th = buildmesh(C(50) + Splus(70) + Sminus(70));
13 
14 // Fespace
15 fespace Xh(Th, P2);
16 Xh psi, w;
17 macro grad(u) [dx(u),dy(u)]// def of grad operator
18 // Solve
19 solve 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    ;
26 
27 plot(psi, wait=1);

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

Potential

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.

PotentialHeat

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

Note

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\);

Note

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