Numerical Tools for the Control of the Unsteady
Heating of an Airfoil

Françoise Masson, Francisco Chinesta, Adrien Leygue, Chady Ghnatios, Elías Cueto, Laurent Dala, Craig Law

► To cite this version:

Françoise Masson, Francisco Chinesta, Adrien Leygue, Chady Ghnatios, Elías Cueto, et al.. Numerical
Tools for the Control of the Unsteady Heating of an Airfoil. Journal of Mechanics Engineering and
Automation, 2013, 3 (6), pp.339-351. hal-01515084

HAL Id: hal-01515084

https://hal.science/hal-01515084

Submitted on 27 Apr 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Numerical Tools for the Control of the Unsteady Heating of an Airfoil

Françoise Masson¹,⁴, Francisco Chinesta¹,², Adrien Leygue¹, Chady Ghnatios¹, Elias Cueto³, Laurent Dala⁴ and
Craig Law⁴

1. EADS Corporate Foundation International Chair, GEM UMR CNRS - Ecole Centrale Nantes, F-44321 Nantes Cedex 3, France

2. Institut Universitaire de France, Paris 75 000, France

3. Aragon Institute of Engineering Research - I3A, Universidad de Zaragoza, Zaragoza E-50018, Spain

4. School of Mechanical, Industrial and Aeronautical Engineering, University of the Witwatersrand, Johannesburg 2000, South Africa

**Abstract:** This paper concerns the real time control of the boundary layer on an aircraft wing. This new approach consists in heating the surface in an unsteady regime using electrically resistant strips embedded in the wing skin. The control of the boundary layer's separation and transition point will provide a reduction in friction drag, and hence a reduction in fuel consumption. This new method consists in applying the required thermal power in the different strips in order to ensure the desired temperatures on the aircraft wing. We also have to determine the optimum size of these strips (length, width and distance between two strips). This implies finding the best mathematical model corresponding to the physics enabling us to facilitate the calculation for any type of material used for the wings. Secondly, the heating being unsteady, and, as during a flight the flow conditions or the ambient temperatures vary, the thermal power needed changes and must be chosen as fast as possible in order to ensure optimal operating conditions.

**Key words:** Model reduction, PGD (proper generalized decomposition), heating of an airfoil, boundary layers, laminar-turbulent transition and separation point, friction drag, unsteady heating.

## 1. Introduction

A new proposal enabling us to control the boundary layer flow over an airfoil is under way [1]. It would influence the boundary layer's laminar-turbulent transition and separation point, allowing an improvement in economic efficiency and safety of airplanes. This new approach proposed is associated with the unsteady surface heating regime using electrically resistant strips embedded in the wing skin. The control of the boundary layer's separation and transition point will provide a reduction in friction drag, and hence a reduction in fuel consumption. The implementation of strips in the wing skin could be

done at a low cost for the manufacturer without weakening the structural integrity of the wing. Another possible advantage of this method is associated with taking-off and landing regimes. This method could enlarge effectiveness of control surfaces and possibly reduce the aerodynamic noise produced by the control surfaces because it could influence the boundary layer separation point.

In order to ensure the desired temperatures on the aircraft wing at a given time, we must determine which is the required thermal power in each strip. This means we have to solve the heat equation for the plate, strip and air surrounding the system. But if we take into account all the equations due to fluid dynamics [2-3], it will be impossible to obtain the result fast enough to heat the airfoil at the temperature wanted.

<footnote>**Corresponding author:** Dr. Françoise Masson, Ph.D. student, research fields: mechanics of materials, processing technology, aerodynamics. E-mail: Francoise.Retat@ec-Nantes.fr.</footnote>
So we have to find a suitable simplified model.

Once we have solved the simplified 2D model, we can go one step further and solve a 3D model. But in this case, the optimal design of the strips has to be worked out: which would be the ideal length, width and distance between two strips in order to obtain the wanted temperatures.

## 1.1 Building-Up Parametric Solutions

Usual models in computational mechanics could be enriched by considering all the sources of variability (e.g., model parameters, initial or boundary conditions, geometrical parameters, etc.) as extra-coordinates. For example, in our case, we are interested in solving the heat equation but we do not know the dimension and position of the source term, as it has to be defined as an optimum between power consumption and evolution of the temperature in time. We have three possibilities: (1) we wait to know the chosen design before solving the heat equation (a conservative solution); (2) we solve the equation for many values of the length, width and distance between two strips and then the work is done (a sort of brute force approach); or (3) we solve the heat equation only once for any length, width and position of the strips.

Obviously the third alternative is the most exciting one. To compute this parametric solution, it suffices to introduce the design parameters as extra-coordinates, playing the same role as the standard space and time coordinates, even if there are no derivatives concerning these extra-coordinates. This procedure runs, very well, and can be extended for introducing many other extra-coordinates: the power of the source term, initial conditions ... (See Ref. [4] and the references therein for an exhaustive and recent review). It is easy to understand that after performing this type of calculations, a posteriori inverse identification or optimization can be easily handled [5-8].

The price to pay is the solution of a model involving many coordinates. If the model is defined in

a space involving $N$ coordinates, standard mesh based discretization techniques require $M^N$ degrees of freedom when $M$ nodes are involved in the discretization of each coordinate. In practical applications with $M \sim 10^3$ and $N \sim 10$ the complexity reaches the value of $M^N \sim 10^{30}$ beyond the present computational capabilities. Thus, standard discretization techniques fail to solve multidimensional models that suffer the so-called curse of dimensionality.

## 1.2 Routes for Circumventing the Curse of Dimensionality

The construction of parametric solutions seems an exciting route, but the main question needs an answer: How can we circumvent the curse of dimensionality?

Different techniques have been proposed for circumventing the curse of dimensionality, Monte Carlo simulations being the most widely used. Their main drawback is the statistical noise. Other possibilities lie in the use of sparse grids [9], within the deterministic framework, but they suffer also when the dimension of the space increases beyond a certain value (about 20).

Separated representations could be a valuable alternative. Separated representations proceeds by expressing a generic multidimensional function $u(x_1, \dots, x_N)$ in a separated form:

$$u(x_1, \dots, x_N) \approx \sum_{i=1}^{i=Q} F_i^1(x_1) \times \dots \times F_i^N(x_N) \quad (1)$$

In this expression, the coordinates $x_i$ denote any coordinate, scalar or vector, involving the physical space, the time or any other extra-coordinate (e.g., the conductivity in the example previously discussed).

Separated representations were present within the Hartree-Fock-based approaches were widely employed in quantum chemistry [10]. In the 80s, space-time separated representations were considered by P. Ladeveze within an original and powerful non-incremental-non-linear solver called LATIN method [11-12]. A natural generalization was proposed by Ammar and Chinesta [13-14] for solving
the highly multidimensional models encountered in
the kinetic theory description of complex fluids.
Parametric models were addressed in Ref. [15].

Thus, if *M* nodes are used to discretize each
coordinate, the total number of unknowns involved in
the solution is *Q* × *N* × *M* instead of the *M*<sup>*N*</sup>
degrees of freedom involved in mesh based
discretizations. We must recall that these functions are
not “a priori” known, they are computed on the fly by
introducing the approximation separated
representation into the model weak form and then
solving the resulting non-linear problem. The
interested reader can refer to Ref. [12] and the
references therein for a detailed description of the
numerical and algorithmic aspects. The construction
of such approximation is called Proper Generalized
Decomposition (PGD) because this decomposition is
not orthogonal but in many cases the number of terms
in the finite sum is very close to the optimal
decomposition obtained by applying the Proper
Orthogonal Decomposition (POD) (or the Singular
Value Decomposition (SVD)) on the model solution.

As can be noticed in the expression of the
approximation separated representation the
complexity scales linearly with the dimension of the
space in which the model is defined, instead of the
exponential growth characteristic of mesh based
discretization strategies. In general, for many models,
the number of terms Q in the finite sum is quite
reduced (few tens) and in all cases the approximation
converges towards the solution associated with a full
tensor product of approximation bases considered in
each space x_i. Thus, we can conclude about the
generality of the separated representation, but its
optimality depends on the solution features.

1.3 *Structure of the Paper*

In this paper, we focus on the numerical tools
linked to the new approach, enabling us to control the
boundary layer flow over an airfoil. First, in section
2.1 we will find a suitable simplified 2D equation,

which represents in a realistic way the heat transfer
equation for the flat plate, strip and boundary layer,
and that will be used in section 2.2 in order to identify
the advection coefficient to be considered for
modeling the thermal boundary conditions avoiding
the necessity of solving the thermal balance in the air
and also for justifying that the ambient temperature
can be considered in all the boundary conditions
instead of the real temperature of the air within the
boundary layer. In section 2.3, in order to determine
the real configuration and optimize the strips location
and geometry we considered a 3D modeling. The
design of the strips (length, width and distance
between two strips) must be taken into account and
introduced as extra-coordinates of the model allowing
efficient optimizations. As the number of
extra-coordinates becomes important, we will use the
PGD method to solve this problem. The last ingredient,
which will be looked at in section 3, is, as mentioned
earlier, the unsteady heating of the strips. We will use
the linearity and the superposition principle to solve
this problem.

**2. Solving the Problem in the Steady State**

## 2.1 A Simple Model for Calculating the Airfoil's Surface Temperature

In what follows, and for simplicity's sake, we solve
the heat transfer equation in a two dimensional
domain, i.e., we consider a flat plate and an electrical
strip resistance on top, both having the same length, as
depicted in Fig. 1.

There is no heat source in the plate, and the source

<figure><img src="image_1.png" /><figcaption>Fig. 1 Heat transfer on a flat plate with an electrically resistant strip on top.</figcaption></figure>
term in the strip is P, thus the equations are respectively

$$ \begin{cases} K_p \Delta U^{II} = 0 \\ K_s \Delta U^{I} = P \end{cases} \qquad (2) $$

where $K_p$ and $K_s$ are respectively the plate's and the strip's conduction coefficient. We assume a local system of the coordinates attached to the bottom boundary in each domain (plate and strip).

The boundary conditions are

• For the plate,

$$ \left\{ \begin{aligned} \left. -K_p \frac{\partial U^{II}}{\partial x} \right|_{x=0} &= h \cdot (U^{II}(0, y) - T_{amb}) \\ \left. K_p \frac{\partial U^{II}}{\partial x} \right|_{x=H_x} &= h \cdot (U^{II}(H_x, y) - T_{amb}) \\ \left. -K_p \frac{\partial U^{II}}{\partial y} \right|_{y=0} &= h \cdot (U^{II}(x, 0) - T_b) \end{aligned} \right. \qquad (3) $$

• For the strip,

$$ \left\{ \begin{aligned} \left. -K_s \frac{\partial U^I}{\partial x} \right|_{x=0} &= h \cdot (U^I(0, y) - T_{amb}) \\ \left. K_s \frac{\partial U^I}{\partial x} \right|_{x=H_x} &= h \cdot (U^I(H_x, y) - T_{amb}) \\ \left. K_s \frac{\partial U^I}{\partial y} \right|_{y=H_y^I} &= h \cdot (U^I(x, H_y^I) - T_a) \end{aligned} \right. \qquad (4) $$

where $h$ is the convection coefficient, $T_{amb}$ represents the ambient temperature, and $T_a$ and $T_b$ respectively the temperatures of the upper and lower boundary layers. $H_x$ is the length of the plate and strip, $H_y^I$ and $H_y^{II}$ are respectively the height of the strip and of the plate.

Between the strip and the plate, we have equality in temperature and thermal fluxes:

$$ \left\{ \begin{aligned} \left. U^{II}(x, y^I = H_y^{II}) \right| &= U^I(x, y^I = 0) \\ \left. K_p \frac{\partial U^{II}}{\partial y} \right|_{y^{II}=H_y^{II}} &= K_s \frac{\partial U^I}{\partial y} \Big|_{y^I=0} \end{aligned} \right. \qquad (5) $$

Furthermore, inside the upper boundary layer, we consider the energy balance in a volume $\Delta x$ during a time $\Delta t$, i.e.,

$$ \frac{h}{K_s} \cdot (T_a - U_{surf}) \cdot \Delta t \cdot \Delta x + v \cdot \rho \cdot c \cdot T_a \cdot \delta \cdot \Delta t - v \cdot \rho \cdot c \cdot T_a^I \cdot \delta \cdot \Delta t = 0 \qquad (6) $$

where $v$ represents the free stream velocity, $\rho$ the

density, $c$ the specific heat capacity, $h$ the convection coefficient and $\delta$ the boundary layer thickness. $U_{surf}$ is the temperature at the surface of the strip, i.e., $U^I(x, H_y^I)$.

Thus, taking into account

$$ T'_{a} \approx T_{a} + \frac{dT_{a}}{dx} \cdot \Delta x \qquad (7) $$

we can deduce that the temperature in the upper boundary layer is governed by

$$ \gamma^I \frac{dT_a}{dx} - T_a = U_{surf} \qquad (8) $$

where $\gamma^I = \rho c v \delta \frac{K_s}{h}$.

Using a similar approach, we obtain the following equation for the lower boundary layer:

$$ \gamma^{II} \frac{dT_b}{dx} - T_b = U_{lower} \qquad (9) $$

where $\gamma^{II} = \rho c v \delta \frac{K_p}{h}$ and $U_{lower}$ is the temperature at the plate's bottom surface.

We have the following boundary conditions for both boundary layers:

$$ T_a(x=0) = T_b(x=0) = T_{amb} \qquad (10) $$

We are looking for the temperature in the whole structure (plate and strip) that we denote by $U(x, y)$, by solving Eqs. (2)-(5) and (8)-(10), however we must first find a simple procedure to estimate the value of the convection coefficient $h$.

## 2.2 A Simple Model for Calculating the Exchange Coefficient

Thus, we consider a second system composed of the

<figure><img src="image_2.png" /><figcaption>Fig. 2 System composed of the flat plate, the strip and the upper and lower boundary layers.</figcaption></figure>
flat plate, the strip and the upper and lower boundary layers, as shown in Fig. 2.

The temperature inside the plate and strip responds to the steady state heat equation as in the first system considered, i.e., Eq. (2). Inside each boundary layer, the temperature is given by Eq. (11).

$$v(x, y) \frac{\partial T}{\partial x} = k \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \quad (11)$$

where $k$ is the air conductivity and $v(x, y)$ is the velocity in the boundary layer, and can be calculated thanks to the Karman equation:

$$v(x, y) = V \left( \frac{2y}{\delta} - \frac{y^2}{\delta^2} \right) \quad (12)$$

where $V$ is the free stream velocity and $y$ is the vertical distance to the solid surface.

The boundary conditions between each sub-domain of this system are

* Between the upper boundary layer and the strip,

$$\left\{ \begin{array}{l} -K' \left. \frac{\partial U^I}{\partial y} \right|_{y=H_y^I} = -k^{T_a} \left. \frac{\partial T_a}{\partial y} \right|_{y=0} \\ \\ U^I(x, H_y^I) = T_a(x, 0) \end{array} \right. \quad (13)$$

* Between the strip and the plate,

$$\left\{ \begin{array}{l} -K^{II} \left. \frac{\partial U^{II}}{\partial y} \right|_{y=H_y^{II}} = -K^I \left. \frac{\partial U^I}{\partial y} \right|_{y=0} \\ \\ U^{II}(x, H_y^{II}) = U^I(x, 0) \end{array} \right. \quad (14)$$

* Between the plate and the lower boundary layer,

$$\left\{ \begin{array}{l} -K^{II} \left. \frac{\partial U^{II}}{\partial y} \right|_{y=0} = -k^{T_b} \left. \frac{\partial T_b}{\partial y} \right|_{y=H_y^{T_b}} \\ \\ U^{II}(x, 0) = T_b(x, H_y^{T_b}) \end{array} \right. \quad (15)$$

Note that in all the equations, the 'y' coordinate is again taken locally, i.e., it is dependent of each layer, and takes its origin ($y = 0$) on the lower surface of the corresponding layer.

If we consider adiabatic boundary conditions along $x = 0$ and $x = H_x$, a quite reasonable assumption because the reduced thickness of the system under consideration, the energy flow is thereby concentrated on the top and lower domain boundaries, and writes

$$\begin{align} & \int_{x=0}^{H_x} -K^I \left. \frac{\partial U^I}{\partial y} \right|_{y=H_y^I} dx + \int_{x=0}^{H_x} K^{II} \left. \frac{\partial U^{II}}{\partial y} \right|_{y=0} dx = \nonumber \\ & \int_{x=0}^{H_x} h \left( U^I(x, H_y^I) - T_{amb} \right) + \int_{x=0}^{H_x} h \left( U^{II}(x, 0) - T_{amb} \right) \tag{16} \end{align}$$

Eq. (16) enables us to estimate the convection coefficient $h$.

Once we have the value of the convection coefficient, we can come back to the simplified model described in the previous section:

* First we solve Eq. (2) with the given boundary conditions, taking an arbitrary $T_a(x)$ and $T_b(x)$ thus obtaining the temperature inside the plate and the strip.

* Then we solve Eqs. (8) and (9) using the temperatures $U^I(x, H_y^I)$ and $U^{II}(x, H_y^{II})$ on the upper and lower surface of the plate just calculated to obtain $T_a(x)$ and $T_b(x)$.

We repeat both steps until convergence.

Now it is possible to determine the value of power $P$ needed in order to achieve the target temperature on the surface of the airfoil simply by doing an inverse calculation. The numerical experiments reported later will prove that the energy balance in both boundary layers is unnecessary and that no significant error is introduced if we consider that everywhere in the boundary layer the air temperature is the ambient one.

## 2.3 Finding the Optimized Location and Dimensions for the Strips

In the previous section, we only considered two dimensions: the length and height of the plate. But eventually, we will have to determine what is the optimum width $W$ of each strip and distance $D$ between two consecutive strips. And this requires solving the problem in 3D, as in Fig. 3. However, thanks to the previous simplified analysis we can justify the use of simplified boundary conditions avoiding the consideration of energy balances in the surrounding air layers. Moreover the exchange coefficient to be considered in the boundary conditions was properly identified. In order to find the
<figure><img src="image_3.png" /><figcaption>Fig. 3 3D representation of the flat plate with the strips on top.</figcaption></figure>

ideal width and distance, W and D, we could attempt different trials until reaching an acceptable value of the cost functions describing the optimal choice. This would imply a 3D solution procedure for each tentative choice of the parameters.

Instead of the traditional procedure, we will define and solve a model considering all the design parameters as extra-coordinates. The resulting multidimensional model will only be solved once, and allows us to define a sort of abacus from which the optimization process can be carried out very fast even on light computing platforms. As we also do not know what the optimal length L of the airfoil will be, neither the power P given to heat the strips, we will add L and P as extra-coordinates. Thus, we arrive to 7 coordinates: x, y, z, W, D, L and P. Using a traditional mesh-based discretization method implies a too high number of degrees of freedom. This is why we use the PGD revisited and summarized in the first section of this work. The system we are considering is the flat plate (represented in white in Fig. 4) with an imbedded electrically resistant strip (in grey in Fig. 4).

This means, we have to solve the following heat transfer equation:

$$k\Delta U = P \qquad (17)$$

where the power source P only applies in the strip.

We have convective boundary conditions in the x and z directions, and symmetry conditions in the y-direction as there would be more than one strip to heat the airfoil, thus, we have the boundary conditions as Eq. (18):

<figure><img src="image_4.png" /><figcaption>Fig. 4 3D representation of the flat plate (in white) one electrically resistant strip (in grey).</figcaption></figure>

$$ \left\{ \begin{aligned} -k \left. \frac{\partial U}{\partial x} \right|_{x=0} &= h (T_{amb} - U_{x=0}) \\ k \left. \frac{\partial U}{\partial x} \right|_{x=L} &= h (T_{amb} - U_{x=L}) \\ \left. \frac{\partial U}{\partial y} \right|_{y=0} &= 0 \\ \left. \frac{\partial U}{\partial y} \right|_{y=D+W} &= 0 \\ -k \left. \frac{\partial U}{\partial z} \right|_{z=0} &= h (T_{amb} - U_{z=0}) \\ k \left. \frac{\partial U}{\partial z} \right|_{z=H} &= h (T_{amb} - U_{z=H}) \end{aligned} \right. \qquad (18) $$

By using the PGD method, we obtain a solution under the separated form:

$$U = \sum_{i=1}^{i=N} X_i(x) \cdot Y_i(y) \cdot Z_i(z) \cdot W_i(W) \cdot D_i(D) \cdot F_i(L) \cdot P_i(P)$$

For more details of such a separated representation construction, the interested reader can refer to Ref. [7]. As we now have the temperature of the airfoil for all width and distance between 2 consecutive strips and for all other extra-coordinates, we can determine the optimum width and distance between 2 strips (as the length of the airfoil is determined by other criteria) according to the power we want to apply on the strips.

### 3. Taking into Account the Unsteady Heating

For simplicity's sake, in the above calculations we only considered steady heating. But usually, the control of the boundary layer flow over the airfoil is achieved by applying unsteady heating in the electrical strips resistances. During the flight, we would need to know which power must be applied and
at which frequency, in order to obtain the desired temperature on the surface of the wing.

We are considering the flat plate/embedded strip system depicted in Fig. 5, which gives us the following parametric heat transfer equation to be solved:

$$ \rho c \frac{\partial u}{\partial t} = k \Delta u + P(x, y, z, t) \quad (19) $$

where the source term $P(x, y, z, t)$ only applies on the strip and has an arbitrary evolution in time.

The boundary conditions are as follows:

$$ \left\{ \begin{array}{l} -k \left. \frac{du}{dx} \right|_{x=0} = h (T_{amb} - u_{(x=0)}) \\ k \left. \frac{du}{dx} \right|_{x=L} = h (T_{amb} - u_{(x=L)}) \\ -k \left. \frac{du}{dx} \right|_{z=0} = h (T_{amb} - u_{(z=0)}) \\ k \left. \frac{du}{dx} \right|_{z=H} = h (T_{amb} - u_{(z=H)}) \\ -k \left. \frac{du}{dy} \right|_{y=0} = 0 \\ k \left. \frac{du}{dy} \right|_{y=W} = 0 \end{array} \right. \quad (20) $$

with the following initial condition:

$$ u(x, y, z, t = 0) = T_{amb} \quad (21) $$

As $P(x, y, z, t)$ can evolve arbitrarily in time but the heat equation being linear, we start by solving the simpler heat transfer equation associated with a unit power-step applied at the initial time:

$$ \rho c \frac{\partial u}{\partial t} = k \Delta u + p \quad (22) $$

<figure><img src="image_5.png" /><figcaption>Fig. 5 Flat plate with an embedded electrically resistant strip and corresponding boundary conditions.</figcaption></figure>

where $p = 1$ in the strip (as shown in Fig. 6), with the same boundary conditions described in Eq. (20), and initial condition described in Eq. (21). This can be done using different methods (finite differences, finite elements or even the non-incremental parametric PGD ...).

Once we have obtained $u(x, y, z, t; p = 1)$ corresponding to the solution of Eq. (22) with the boundary conditions (20), we can obtain $U(x, y, z, t; P)$ solution of the equation:

$$ \rho c \frac{\partial u}{\partial t} = k \Delta u + P \quad (23) $$

using the linearity, simply by multiplying $u(x, y, z, t; p = 1)$ by $P$, i.e.,

$$ U(x, y, z, t; P) = P \cdot u(x, y, z, t; p = 1) \quad (24) $$

The solution related to an arbitrary evolution of the source term can be now easily obtained by invoking the superposition principle.

If $P(t) = \sum_i \theta_i \cdot H_i(t)$, as seen in Fig. 7, by applying the superposition it results,

<figure><img src="image_6.png" /><figcaption>Fig. 6 Source term $p = 1$ for $t > 0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7 Source term $P$ varying in time.</figcaption></figure>
$$U(x, y, z, t) = \sum_{i} \theta_{i} \cdot u(x, y, z, t - t_{i}; p = 1) \quad (25)$$

# 4. Numerical Results

## 4.1 Steady Heating

### 4.1.1 Airfoil's Temperature in 2D

As explained in section 2.1, we first have to estimate the convection coefficient. This coefficient $h$ is assumed in first approximation independent of the power induced in the strip and the outside temperature. The only parameter that can influence this coefficient is the free-stream velocity $v$ as the conditions are then comparable to a forced convection.

The expected tendency is reflected in the numerical results, as $h$ changes only according to the free-stream velocity.

For $v = 20$ m/s, we obtain $h = 8.2$ W/m² K.

For $v = 100$ m/s, we obtain $h = 18.2$ W/m² K.

For a flat plate 15 cm long, 17.9 mm thick, and a strip of the same length and 0.1 mm thick, solving Eqs. (2), (6) and (7) with the boundary conditions described in Eq. (3), gives us the following temperatures (in the following figures, “I” represents the intensity administered in the electrically resistant strip).

Figs. 8-13 show the influence of the conduction coefficient $k$ and the intensity $I$ applied on the strip. Figs. 14-18 compared to the previous figures let us

<figure><img src="image_8.png" /><figcaption>Fig. 8 Plate's temperature for $v = 20$ m/s, $h = 8.2$ W/m²K, $I = 2$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 9 Boundary layer's temperature in contact with the plate, for $v = 20$ m/s, $h = 8.2$ W/m²K, $I = 2$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10 Plate's temperature for $v = 20$ m/s, $h = 8.2$ W/m²K, $I = 7$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11 Plate's temperature for $v = 20$ m/s, $I = 2$ A and $k = 100$ W/mK.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 12 Boundary layer's temperature in contact with the plate, for $v = 20$ m/s, $I = 2$ A and $k = 100$ W/mK.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13 Plate's temperature for $v = 20$ m/s, $I = 7$ A and $k = 100$ W/mK.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 14 Plate's temperature for $v = 100$ m/s, $h = 18.2$ W/m$^2$K, $I = 2$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 15 Boundary layer's temperature in contact with the plate, for $v = 100$ m/s, $h = 18.2$ W/m$^2$K, $I = 2$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 16 Plate's temperature for $v = 100$ m/s, $h = 18.2$ W/m$^2$K, $I = 7$ A and $k = 300$ W/mK.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 17 Plate's temperature for $v = 100$ m/s, $I = 2$ A and $k = 100$ W/mK.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 18 Boundary layer's temperature in contact with the plate, for $v = 100$ m/s, $I = 2$ A and $k = 100$ W/mK.</figcaption></figure>

appreciate the influence of the free stream velocity.

### 4.1.2 3D Plate Temperature—PGD Solution

In this example, we are going to solve the 3D heat transfer equation for a flat plate, i.e., Eq. (10) with the boundary conditions described in Eq. (11) for 7 coordinates: $x$, $y$, $z$ (space coordinates), the width of the strip $W$, the distance between two strips $D$, the length $L$ and the power in the strip $P$. Thus, the solution is found under the form:

$$U = \sum_{i=1}^{i=N} X_i(x) Y_i(y) Z_i(z) W_i(W) D_i(D) F_i(L) P_i(P) \quad (26)$$

The convection coefficient $h$ and the conduction coefficient $k$ are considered known. The results are shown for two values of each of these parameters: $h = 8.2$ W/m²·K and $h = 18.2$ W/m²·K; $k = 300$ W/mK and $k = 100$ W/mK.

The calculation time in order to obtain the solution for all values of the extra-coordinates is approximately 142s. This abacus contains the temperature of all the points on the plate and strip, for all length of the airfoil, width of the strips, distance between two strips and for any power applied. The meshing chosen was 100 values for $x$, 200 for $y$, 500 for $z$, 100 for $P$, 50 for $L$, and 20 for $W$ and $D$, that means $20 \times 10^{13}$ dof with a traditional meshing technique, or having to solve $2 \times 10^6$ 3D problems in order to obtain the same information as Eq.

<figure><img src="image_19.png" /><figcaption>Fig. 19 Plate's surface temperature for $h = 8.2$ W/m²·K, $k = 300$ W/mK, $W = 0.289$ m, $d = 0.289$ m, $P = 5 \times 10^6$ W/m³.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 20 Plate's surface temperature for $h = 8.2$ W/m²·K, $k = 300$ W/mK, $W = 0.1$ m, $d = 0.48$ m, $P = 5 \times 10^6$ W/m³.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 21 Plate's surface temperature for $h = 8.2$ W/m²·K, $k = 300$ W/mK, $W = 0.289$ m, $d = 0.289$ m, $P = 45.5 \times 10^6$ W/m³.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 22 Plate's surface temperature for h = 8.2 W/m²K, k = 100 W/mK, W = 0.289 m, d = 0.289 m, P = 5 × 10⁶ W/m³.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 23 Plate's surface temperature for h = 18.2 W/m²K, k = 300 W/mK, W = 0.289 m, d = 0.289 m, P = 5.10⁶ W/m³.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 24 Plate's surface temperature for h = 18.2 W/m²K, k = 100 W/mK, W = 0.289 m, d = 0.289 m, P = 5 × 10⁶ W/m³.</figcaption></figure>

(26) contains. The parametric solution particularization only takes 0.3 s.

From now on, in order to enable comparison, we only show the plate's surface temperature in Figs. 19-24. But the temperature of any part of the plate could be shown if wanted.

We can clearly notice in Figs. 19-24 the influence of the convection coefficient (which changes according to the free-stream velocity), the conduction coefficient, the applied power on the strip and the width of these strips.

## 4.2 Unsteady Heating

In the presented example, we consider a flat plate with the following dimensions:

$$x \in [0; 15], y \in [0; 20], z \in [0; 0.18]$$

The outside temperature is considered to be $T_{amb} = 20$ °C.

The heating strip is located at $x \in [3; 12], y \in [4; 8], z \in [0.16; 0.18]$.

We first calculate the solution of $\rho c \frac{\partial u}{\partial t} = k \Delta u + p$ where $p = 1$ as previously described.

Fig. 25 represents the temperature on the surface of the plate for $P = 1$. To obtain the solution for $P = 2$, as shown in Fig. 26, we make use of the linearity.

The temperature on the surface of the plate can easily be seen for one step, two or more, by using superposition.

Fig. 27 shows the temperature evolution at the surface of the strip ($z = 0.18$ cm), where $x = 7.5$ cm and $y = 6$ cm with one power step at $t_1 = 1$ and $\theta_1 = 6$.

Fig. 28 shows the temperature of the same point and same thermal step than in Fig. 27, but with a second step applied at $t_2 = 4$ with $\theta_2 = -6$. We can clearly see the temperature of the plate decreasing as a result of this second impulse. In Figs. 29-30, we can see the temperature at the same point under the influence of a pulsed wave, as we apply $P$ and $-P$ alternatively at every time step. In Fig. 29, the pulse frequency is $\Delta t = 3$, and $P = 6$. In Fig 30, $P$ is doubled, and $\Delta t$ is divided by 2.
<figure><img src="image_25.png" /><figcaption>Fig. 25 Temperature on the surface of the plate, at time $t = 50$, for $P = 1$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 26 Temperature on the surface of the plate, at time $t = 50$, for $P = 2$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 27 One impulse at time $t_1 = 0$ and $\theta_1 = 6$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 28 Two Impulses at $t_1 = 0$ where $\theta_1 = 6$, and at $t_1 = 25$ where $\theta_1 = -6$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 29 Impulse at $t_1 = 0$ and $\theta_1 = 6$ and at $t_2 = 25$ and $\theta_2 = -6$. Thereafter pulse wave every $\Delta t = 3$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 30 Impulse at $t_1 = 0$ and $\theta_1 = 12$ and at $t_2 = 25$ and $\theta_2 = -12$. Thereafter pulse wave every $\Delta t = 1.5$.</figcaption></figure>
## 5. Conclusions

A new approach to improve fuel consumption during an airplane flight consists in heating the wings. This requires the creation of a numerical abacus in order to know which thermal power must be given to obtain the required temperature as a function of the flight conditions.

In this work, we have seen that it is possible to obtain a simple mathematical model corresponding to the physics: indeed, the 2D models enabled us to define the convection coefficient *h*, and determine that it is not compulsory to take the boundary layer's temperature into account to obtain an accurate temperature of the surface of the plate. Instead, it suffices to consider the ambient temperature.

The 3D parametric model helps the design of the airfoil and strip by enabling the evaluation of all design configurations thanks to a faster calculation time. And finally, the 3D unsteady heating model (parametric or not) allows fast unsteady calculus for any given signal.

The ongoing work concerns

* We must still define the optimum and control criteria: As the final aim is to improve fuel consumption during a flight, we must determine which would be the compromise between the energy consumption and temperature, both depending on the intensity, length and power step frequency;

* A control system based on the numerical PGD abacus (containing all the different parameters) ensuring the right temperature on the surface of the strip, i.e., the inverse identification, must still be defined.

## References

[1] F. Masson, F. Chinesta, A. Leygue, E. Cueto, L. Dala, C. Law, Real-time control of the heating of an airfoil, Journal of Materials Science and Engineering A 2 (5) (2012) 478-487.

[2] Dr.H. Schlichting, Boundary-Layer Theory, McGraw-Hill, 1968.

[3] J.D. Anderson, Fundamentals of Aerodynamics,

McGraw-Hill, 1985.

[4] F. Chinesta, P. Ladeveze, E. Cueto, A short review in model order reduction based on Proper Generalized Decomposition, Archives of Computational Methods in Engineering 18 (2011) 395-404.

[5] F. Chinesta, A. Ammar, E. Cueto, Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models, Archives of Computational Methods in Engineering—State of the Art Reviews 17 (4) (2010) 327-350.

[6] Ch. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, P. Breitkopf, P. Villon, Methodological approach to efficient modelling and optimization of thermal processes taking place in a die: Application to pultrusion, Composites Part A: Applied Science and Manufacturing 42 (2011) 1169-1178.

[7] Ch. Ghnatios, F. Masson, A. Huerta, E. Cueto, A. Leygue, F. Chinesta, Proper Generalized Decomposition based dynamic data-driven control of thermal processes, Computer Methods in Applied Mechanics and Engineering 213-216 (2012) 29-41.

[8] D. Gonzalez, F. Masson, F. Poulhaon, A. Leygue, E. Cueto, F. Chinesta, Proper Generalized Decomposition based dynamic data-driven inverse identification, Mathematics and Computers in Simulation 82 (2012) 1677-1695.

[9] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 1-123.

[10] E. Cancès, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational quantum chemistry: A primer, Handbook of Numerical Analysis X (2003) 3-270.

[11] P. Ladeveze, Nonlinear Computational Structural Mechanics, Springer, NY, 1999.

[12] P. Ladeveze, J.-C. Passieux, D. Neron, The LATIN multiscale computational method and the Proper Generalized Decomposition, Computer Methods in Applied Mechanics and Engineering 199 (21-22) (2010) 1287-1296.

[13] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids, J. Non-Newtonian Fluid Mech. (2006).

[14] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II. Transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech. 144 (2007) 98-121.

[15] E. Pruliere, F. Chinesta, A. Ammar, On the deterministic solution of multidimensional parametric models by using the proper generalized decomposition, Mathematics and Computer Simulation 81 (2010) 791-810.