Modélisation de l'évaporation d'un film de carburant




Evaporation of fuel layer




Geremia Giuliano alejandro

Mezguledi Ansar

Quentin Mouret



Legendre Dominique


March 2013






This site is a part of the work in the Industrial Engineering consulting  Energetic and Process.

The Industrial Engineering consulting BEI is a long project implemented within the Hydraulic Department and Fluid mechanics of ENSEEIHT by its various specialities Process and Energetic (E&P) and Computational fluid dynamics (MFN).

The aim of these BEI is to work in groups of 3 students on scientific projects proposed by industrial collaborators.

During our BEI, we will provide a study on Evaporation of fuel layer by elaborating a physical problem modeling. We will establish a theoretical analysis by treating all the physical phenomenon relevant to the problem such as diffusion and convection. The second part of the BEI is to implement a numerical simulation using OpenFoam for several configurations of this study.




The team that is working on this project comprises three students:

          Geremia Giuliano Alejandro: 3rd year at ENSEEIHT, Computational fluid dynamics (MFN)

          Mezguledi Ansar                   : 3rd year at ENSEEIHT,  Process and Energetic (E&P)

          Quentin Mouret                    : 3rd year at ENSEEIHT, Process and Energetic (E&P)




This BEI is supervised by Mr. Dominique Legendre a research and teaching professor at IMFT-ENSEEIHT.



Industrial environment and objectives

Industrial environment                          

& Objectives   


          1- Industrial environment

    The aeronautical sector is at present an ultra-competitive sector. With the arrival of new Russian and Chinese competitors, both big actors of this sector, Airbus and Boeing, have to step up their ingenuity and innovation to keep their leaderships and to deal with new markets.

    So they are developing at present new models of planes which are lighter, more aerodynamic and less consumer. They aim maintain a distance from their competitors and the new incomers.


le - Publié le

     It is in this rough economic environment that Airbus requested this study in collaboration with Altran, the IMFT (Institute of Fluid mechanics of Toulouse) as well as the Hydraulic Department and Fluid mechanics of the ENSEEIHT.


           2- Industrial partners

  • Airbus

    Airbus SAS is an aircraft manufacturing subsidiary of EADS, a European aerospace company. Based in Blagnac, France, a suburb of Toulouse, and with significant activity across Europe, the company produces approximately half of the world's jet airliners.



    Airbus employs around 63,000 people at sixteen sites in four European Union countries: France, Germany, the United Kingdom and Spain. Final assembly production is based at Toulouse, France; Hamburg, Germany; Seville, Spain; and, since 2009, Tianjin, China. Airbus has subsidiaries in the United States, Japan, China and India.
The company produces and markets the first commercially viable fly-by-wire airliner, the Airbus A320, and the world's largest airliner, the A380.


  • Altran Technologies


    Altran Technologies, SA is a European consulting firm founded in 1982 in France. Altran, self dubbed the European leader in high technology and innovation consultancy, operates primarily in technology and innovation consultancy, accounting for nearly half of its turnover. Administrative and information consultancy accounts for a third of its turnover with strategy and management consulting making up the rest.


Mission & Planning

Mission & Planing                                 


           1- Mission

    Under this project, it was asked to realize the study of a thin layer evaporation. This one is in line with the necessity of Airbus to reduce the weight of its planes.

Indeed, there are, at present, fuel leaks in the piping which allows to forward the fuel of the reservoir up to the engines of the turbomachinery of planes. This fuel is got back by reservoirs which are oversized.



                                                 Figure: Position of the studied geometry in red


Moreover, the flammability increases according to the concentration of the fuel in the air. It is for this that it will be necessary to limit the values of this concentration.

This project aims to study the evaporation of a thin layer of the liquid to be able to re-size reservoirs getting back fuels and this to decrease recirculations there.

 Both theoretical and numerical study on the evaporation of a layer of fuel in such reservoir will thus be performed in this BEI.


           2- Project's organization

​The following organization chart shows the principal activities that we have treated in this BEI. 


                                                          Figure: Project's organization chart



Theoretical analysis

Theoretical analysis                           


The first part of this project aims to achieve a theoretical analysis of the fuel layer evaporation in the container. The next sub-parts correspond to the steps that we have followed to determine the equations system solving the problem.

First, we will show the assumptions that we take to model the problem, then we will study some particular cases before considering the evaporation in the whole container.













   Preliminary studies on a particular case

This first study involves to work on evaporation in enclosure where a  pure macro-layer is lied. This enclosure is considered adiabatic except the bottom wall where its temperature is fixed. For the first part, the study is in 2D but the edge effect along the vertical wall is neglected.

The second study involves to work on a particular case with different boundary layer : thermal and mass boundary layer.


Problem modeling

Problem modeling                             


       1- Lists of symbols & subscripts


Here are the list of all symbols and subscripts used in the study.

  • $h$    liquid thickness
  • $L$    container height dimension
  • $l_c$    capillary length
  • $a$    depth
  • $b$    length
  • $S_{flux}=a b$    flux surface ($m^2$)
  • $\rho$    density of the gas mixture ($kg/m^3$)
  • $P$    pressure ($Pa$)
  • $u$    the projection of velocity vector along the x-axis ($m/s$)
  • $\mu$    the dynamic viscosity of the gas mixture  ($Pa.s$)
  • $\nu$    the cinematic viscosity of the gas mixture ($m^2/s$)
  • $v$    the projection of velocity vector along the y-axis ($m/s$)
  • $g$    standard gravity ($m/s^2$)
  • $C_p$    specific heat ($^{-1}.K^{-1}$)
  • $T$    temperature of the mixture ($K$)
  • $\lambda$    the thermal conductivity ($W/m.K$)
  • $D$    diffusion coefficient of the vapour in the air ($m^2/s$)
  • $\rho_S$    mass concentration of vapour at the interface ($kg/m^3$)
  • $\omega$    mass fraction of vapour in the air with $\omega= \frac{\rho_v}{\rho_v+\rho_a}$
  • $\dot{m}$    mass flow rate ($kg/s$)
  • $h_{lg}$    latent heat ($kJ/kg$)
  • $h_c$    thermal exchange coefficient
  • $h_m$    mass exchange coefficient
  • $D_h$    hydraulic dimension

                Subscript :

  • $v$    for the vapour
  • $S$    for the vapour at the interface
  • $in$    for the input condition
  • $w$    for the wall condition
  • $l$   for the liquid


        2- Assumptions related to the vapour- state :

We consider a layer of liquid in a ventilated box. The following figure represents the problem:



                                     Figure: Problem modeling & Simplifications


During this study, we simplify the problem by considering the following assumptions:

  • Bidirectional incompressible turbulent flow.
  • Liquid film thin in relation to the reservoir but sufficient thick for neglect the interfacial resistance.
  • Air charged in vapour is in thermodynamic balance : the phase change happens in saturation conditions.
  • Viscous dissipation and pressure work is neglected.
  • Soret effect ( temperature gradient dependence for the mass flow) and Dufour effect (mass gradient dependence for the heat flux) neglected.
  • Radiation phenomena neglected.
  • Boussinesq approximation considered.

In further, the conservation equations are presented below.




       1- Non-dimensional equations

For a greatest lisibility of the mechanisms, the non dimensional equation can be written as :

$$X=\frac{x}{L} , Y=\frac{y}{L} , U=\frac{u}{U_{in}} , V=\frac{v}{V_{in}}, D_{hydraulic}=L$$ ,

$$ \theta=\frac{T-T_{in}}{T_s-T_{in}} ,  C=\frac{\omega-\omega_{in}}{\omega_s-\omega_{in}}$$

The non dimensional equations become :

$$U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=\frac{1}{Re}\left(\frac{\partial^{2}U}{\partial X^{2}}+\frac{\partial^{2}U}{\partial Y^{2}}\right)    (1*)$$

$$U\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial Y}=-\frac{\partial P_{m}}{\partial Y}+\frac{1}{Re}\left(\frac{\partial^{2}V}{\partial X^{2}}\right)      (2*)$$

$$U\frac{\partial \theta}{\partial X}+V\frac{\partial \theta}{\partial Y}=\frac{1}{RePr}\left(\frac{\partial^{2}\theta}{\partial X^{2}}+\frac{\partial^{2}\theta}{\partial Y^{2}}\right)    (3*)$$

$$U\frac{\partial C}{\partial X}+V\frac{\partial C}{\partial Y}=\frac{1}{ReSc}\left(\frac{\partial^{2}C}{\partial X^{2}}+\frac{\partial^{2}C}{\partial Y^{2}}\right)      (4*)$$


        2- Non-dimensional numbers

Some interesting numbers appear in this non-dimensional equation.

  • Reynolds number

The Reynolds number is characteristic of pure flow conditions, inertial effect compared with momentum diffusion:


  • Prandlt number

Another one is the Prandlt number. It appears in thermal equation and it compares the momentum diffusion with the thermal diffusion

$$Pr=\frac{\frac{\mu}{\rho_{in}}}{\frac{\lambda}{\rho C_p}}$$

  • Schimdt number

The Schimdt number, quiet close with the Prandt number but it appear in mass equation. It compares the movement quantity diffusion with the molecular diffusivity.


In addition the expression of the thermal flows through the liquid film has to be written clearly because the study of the heat transfers can explain the energy transfer between the liquid phase and the vapour phase.

Firstly, the sensible heat flux can be modeled as:

$$ q''_s=-\lambda_{liq} \left(\frac{\partial T}{\partial y}\right)_{y=h}          (5) $$

That is the energy transfer through the thickness of the liquid  resulting from temperature gradient. We can suppose in this case the convection is neglected.


Then the local energy balance on the interface is :

$$q''_{evap}=q''_{fluid}+q''_ {gas}$$

that is to say :

$$h_{lg}h_m(\rho_{in}-\rho_{S})=h_c(T_{in}-T_S)+\lambda_l \left(\frac{\partial T}{\partial y}\right)_{y=h}$$

  • Nusselt number & Sherwood number

After all that, the Nusselt number and the Sherwood number are interesting to be studied. These numbers compare the convective effect and the diffusion effect.

The local Nusselt number compares the thermal convection flux with the thermal diffusion flux and it is defined as :

$$ Nu=\left(\frac{\partial \theta}{\partial Y}\right)_{Y=h/L}$$

Then we can deduce the global Nusselt Number as :



To finish this first part, the Sherwood number compares mass convective flux with molecular diffusion.

$Sh=\frac{h_m*D_h}{D}$ with $h_m=\frac{\dot{m}}{\rho(\omega_{in}-\omega_{S})}$

The total flux of vapour mass is $\dot{m}=\rho_l v_l=\rho_S v_S=\omega_S\rho v_S-\rho D(\frac{\partial \omega}{\partial y})_{y=h}$.



This work was made according to the thesis "Transfert couplé de chaleur et de masse par convection mixte avec changement de phase dans un canal" fulfilled by Othmane OULAID.

Diffusion dominant position

Diffusion dominant position                


In this part, the context is the same as in the precedent part  but the air charged in vapour is fixed. Then the molecular diffusion is the only mechanism which transports the vapour into the air.

For this, we will simplify the developed equations in the previous section depending on this specific case.


         1- Mass balance equation 

$$ \frac{d}{dt}(hDl\rho)=-\rho_S v_S S_{flux}$$

In this equation the liquid volume variation is equal to the vapour flux out going from the surface.


          2- Interface mass balance equation 

At the surface there is :

  • a production of vapour $\rho_S v_S S_{flux}$,
  • a transport by convection and diffusion (Fick's law)

The convection velocity is $v_c=(1-\omega)v_a+\omega v_v$ with $v_a$ air velocity and $v_v$ vapour velocity. The convection velocity in the case of fixed air is $v_c=\omega v_v$

At the surface $v_{cS}=\omega_S v_S$.


The output vapour is considered on the interface as : $\rho_S v_{cS} S=\rho_S \omega_S v_S S$.

The balance equation corresponds to "production=diffusion+evaporation",  then :

$$\rho_S v_S S=-D_S \rho_S (\frac{\partial \omega}{\partial y})_{y=h} S+\rho_S\omega_S v_S S$$

We then develop this expression and we obtain:

$$v_S=\frac{D_S (\frac{\partial \omega}{\partial y})_{y=h}}{\omega_S-1}$$

Afterwards the Spalding parameter $b=\frac{\omega}{\omega_S-1}$ then $v_S=D_S(\frac{db}{dy})_{y=h}  (8)$

After calculations, it's simplified to:

$$j_S=\rho_S \frac {D_S}{L-h}ln(1+(b_L-b_S))$$

and becomes:

$$j_S=\rho_S \frac {D_S}{L-h}ln(1+\frac{(\omega_S-\omega_L)}{1-\omega_S})$$


         3- Conclusion

We can observe on this evaporation model that the flux depends on the boundary conditions and thermodynamic balance.

The model is available for a mono-constituent. For a particular case we can deduce the thermodynamic balance parameters $T_{Sat}=T_{wall}$ and $\rho_{Sat}$ with tables.

The comparison inside the logarithm from the last formula ( the sign of the difference between $b_L$ and $b_S$)  gives us the sign of the flux:

  • If the mass concentration is more important far from the interface than the mass concentration at the interface, then flux is negative : there is condensation.
  • Else if $b_S>b_L$ there is evaporation. Else $b_L=b_S$ the system is in total balance state.

Moreover, we can deduce a characteristic time of the problem:


To have an order of magnitude, here there's a numerical application in the water case : 

$T_w=300K, \rho_S=2.38*10^{-2}kg/m^3, L=50 cm, D_S=5*10^{-5} m^2/s$

we obtain $\tau=60$ hours.

Aeration dominant position

Aeration dominant position                


        Evaporation leaded by the heat science

If the air velocity is very important, we can simplify the problem, it may still be necessary to consider that the convection transport is very effective. Then, in that case, it is the conduction inside the fluid which could master the mechanism. In other words, the convection takes the vapour faster than the vapour is produced.

At the interface, we can made a thermal flux balance. As we know, the interface is in a balance state :

$$ 0=q_s+q_L$$

$$0=-\lambda_{liq} S_{flux} \frac{\partial T}{\partial y}\mid_{y=h}+\dot{m_{vap}}h_{fg}$$

with that we can write

$\dot{m}=\rho_l v_l S_{flux}=\rho_s v_s S_{flux} $

and $v_s$ is a boundary condition for the air at the interface liquid-vapour. This is brought by the mass conservation.

Now we can write $v_l=-\frac{d h}{d t}$

with $h$ the effective length of the liquid film.

We have $q_s=q_l$

afterwards : $$\frac{\lambda_l \Delta T}{h}=\rho_l \frac{d h}{d t} h_{fg}$$

We can deduce : $$h^2(t)={h_0}^2-\frac{2\lambda_l \Delta T}{\rho_l h_{fg}} t$$

with a characteristic time $$ \tau=\frac{h_0 \rho_l h{fg}}{2 \lambda_l \Delta T}$$


           Numerical Application :

for $\Delta=T_{wall}-T_{sat}=1 K, h_0=100\mu m, h_{fg}=2400 kJ/kg, \lambda_l=0,67 W/Km, \rho_l=1000 kg/m^3$.

We obtain:

$$\tau=17,9 s$$






Fuel evaporation in the container

Fuel evaporation                                       

in the container


Evaporation rate model

Evaporation rate model                      


  In the previous parts, the general equations and the global context have been presented, as well as two examples of application in particular situations :

  • pure diffusion 
  • pure aeration

In this part, the industrial case is studied. This work requires both a theoretical and a numerical study in order to identify the parameters which give the mass concentration profile throughout the container. The aerodynamic equations is intrinsically linked with thermal and mass coefficient. In further, the link between the thermal and the mass mechanism is developed. As referencing with the first part, the Sherwood and the Nusselt number appear in the equations in their dimension form $h_c$ and $h_m$.


       1- Mass conservation equation

$$ab\rho v= ab\rho_s v_s   (1)$$


     2- Spalding mass transport equation

$$\rho v \frac{\partial b}{\partial y}=\frac{\partial}{\partial y}\left(\rho D\frac{\partial b}{\partial y}\right)   (2)$$

with the spalding number $b=\frac{\omega}{\omega_S-1}$


         3- Heat transport equation 

$$\rho v \frac{\partial  C_p T}{\partial y}=\frac{\partial}{\partial y}\left(\lambda\frac{\partial   T}{\partial y}\right)   (3)$$

$$\lambda_s\left(\frac{d T}{d y}\right)_s=h_c(T_\infty-T_s)$$

\[h_c(T_\infty-T_s)=\rho_s v_s L_v\]


          4-  Boundary conditions

The boundary conditions are expressed as follows:

  • $b=b_{\infty}$ for $y=\delta_\omega$
  • $b=b_{s}$ for $y=h$
  • $v_s= D_s \left(\frac{d b}{dy}\right)_s$
  • $T=T_\infty$ for $y=\delta_T$
  • $T=T_s$ for $y=h$


Therefore, we can simplify the previous expressions some simplifications, the result is :

$W_v(x)=b\rho_sD_s ln(b_{\infty}-b_s+1)\frac{1}{\delta_{\omega}(x)-\delta}$

$$W_v(x)=\frac{b\lambda}{C_p} ln\left( 1+\frac{W_v(x)}{b}\frac{C_p}{h_c(x)}\right)\frac{1}{\delta_T(x)-\delta}$$


The result is close to the droplet result of J.Reveillon ( Rouen University). The second equation is entirely dependent on thermal conditions. That is why we can have an expression of the mass flux of vapour with the values of the thermal boundary layer $\delta_T$ and the thermal convection coefficient $h_c$. Then we can deduce that the mass flux influences the mass boundary layer.


Parametric study

Parametric study                                    


In this section, we will use the expressions of different parameters that we've found to compute these parameters in multiple flight's cases.


                  1- Physical quantities

We need the values of the diffusion coefficient in the air and the saturation vapour pressure to can compute our parameters. 

The following table summarizes the values of these physical quantities in different flight's cases:




                 2- Parametric study

Now we can calculate the non dimensional numbers involved in this problem. We summarize all the computed values in the next table for different flight's cases:







Numerical simulations

Numerical simulations                          


Numerical simulations consist of two parts. First, we will simulate the simple case of evaporation trough a flat-plate in order to estimate the thickness of the boundary layer.  Then, we will treat the real problem of the BEI by simulating the evaporation inside the reservoir.

In fact, the numerical analysis of the reservoir evaporation will be splitted in three parts.  First a hydraulic study will be done on the fuel recuparator using the solver simpleFoam from OpenFoam. Then OpenFoam will be used to create a solver which is able to simulate correctly the evaporation of a droplet and of a liquid film.
Finally, we will try to couple these two solvers.

The following paragraph represents a brief overview of OpenFoam.

         OpenFoam Software

OpenFoam (Open Source Field Operation And Manipulation) is a software written in C++ language. It's dedicated to resolve problems of Continuum mechanics and fluids mechanics. It was implemented on 1980 in Imperial college in London in order to improve a new simulation platform which will be stronger and more flexible than that used till now (Fortran).

So OpenFoam's use is more complex than others industrial software. In facts it's able to solve complex fluid flows using the "Standards solvers" it contains. Therefore, OpenFoam can solve CFD basic codes as well as compressible, incompressible or multiphase flows problems. OpenFoam solvers allow also solving problems of combustion, heat transfers and electromagnetism... More thann that, OpenFoam's user can easily establish a new solver by modifying an existing solver or by coding its own solver.

OpenFoam offers "Standards utilities" which are used to generate 3D meshes, to convert these meshes to OpenFoam files format in order to manipulate it using a lot of options. OpenFoam offers also a rich tools for post-processiong the velocities fields and the scalar fields...

Moreover OpenFoam has a large library "Standard library". It proposes multiples models that the user can choose as RAS and LES turbulence models ...

In this study, we will use some standard solvers,in particular "SimpleFoam" and "IcoFoam". We also create new solver "MyIcoFoam" that describes the evaporation problem.

To have more details about OpenFoam, here is the following link:





Boundary layer test case

Numerical simulations                                             


The difficulty in CFD is to have correct simulation results. When a new Navier Stokes code is used, it has to be tested on some academical cases. Usually, the turbulent codes (RANS) are tested on the axysymmetric impicting flow and the laminar cases are tested on the boundary layer cases.

As it has been described before, the main point here is to give reliable thermal boundary layer thickness in order to estimate the evaporation fluxes at the bottom of the geometry. The Navier Stokes solver adapted for such a goal is Icofoam: a transient solver for incompressible, laminar flow of Newtonian fluids. As there is no equation of temperature in this solver, it has been hanmade implemented.

Once the solver is chosen, it has to be validated on a test case: the boundary layer case.










Problem position

Numerical simulations                          


       1- Problem's geometry

In this part, we simulate the boundary layer over a flat-plate. For this, we consider a flat plate of length L. This flat-plate is situated in a uniform flow with a uniform velocity U.

                                               Figure: Flat-plate in a uniform flow



       Flat-plate's mesh

We use OpenFoam tools to mesh the domain. The mesh is particularly refined above the plate in order to describe better the evolution of the boundary layer over the flat-plate. The global domain is divided to (100x200) meshes.

The picture below represents the meshed domain viewed by Paraview.

                                                    Figure: Domain mesh


       3- "MyIcoFoamT" solver

OpenFoam doesn't offer any standard solver that can resolve both thermic and dynamic equations. That's why, we create our own solver "MyIcoFoamT" in order to simulate the thermic boundary layer. In fact, we modify the standard solver "IcoFoam" that solves the incompressible laminar Navier - Stokes equations using the PISO algorithm, by adding temperature equations.

The next screenshots are  representation of modifications that we have done in the source code:




                                                     Figure: Screenshot of source code


       4- Test case

We study the boundary layer over a flat-plate under the following circumstances:

  • The flat-plate is situated in a uniform flow of air with a uniform velocity U=0.025m/s. We consider that the air kinematic viscosity is equal to 0.55 10-6 m2/s. Therefore, the Reynolds number is approximately 327.15. Thus, the test case is a laminar flow.
  • The plate's temperature is equal to Tp=90°C while the air temperature is T0=60°C.

The next figure shows the domain at initial instant of simulation

                          Figure:  the domain at initial instant of simulation





Estimation of the boundary layer thickness

Numerical simulations                            


      Estimation of the boundary layer thickness​


Once we have created the geometry and then generated the refined mesh of the domain above the flat-plate, we can launch the simulation to visualize the boundary layer over the flat-plate. We will now expose the results obtained.

We will first study the dynamic boundary layer, then the thermic boundary layer. The aim of this project's part is to estimate the evolution of the boundary layer thickness and to compare it with theoretical profiles corresponding to Blasius solutions.





Dynamic boundary layer

Numerical simulations                         


      1- Dynamic boundary layer


                  a- Velocity field

In this part, we consider that the fluid and the plate have the same temperature. The physical problem is then reduced to a dynamic study.
The following figure shows the evolution of the dynamic boundary layer over the flat plate.             


                          Figure: The evolution of the dynamic boundary layer over the flat plate


We notice that the boundary layer is quickly established trough the area over the plate. The velocity values increase as we go away from the plate.

Moreover, it is noticeable that:

  • All fluid particles in touch with the wall are immobile relatively  to the plate.
  • The flow near the plate is slowed down.
  • The normal gradient of velocity near to the wall is important due to the viscosity effect.
  • The viscosity, in the boundary layer, is small but it has an important impact on the shear stress on the plate $$ \tau_p = \mu \frac{\partial U}{\partial y} |_{paroi} $$  which could have high values.


              Velocity profile

In the figure below we plot the velocity profile versus the altitude above the plate in x= 0.9 cm.


                      Figure:  the velocity profile versus the altitude above the plate in x= 0.9 cm


                   b- Boundary layer thickness

The most important in this study is to estimate the boundary layer thickness values over the flat plate. We define this thickness as following  $$ U[x,\delta(x)] = 0.99  U_0 $$ .


             Blasius model for a flat plate



                                           source: "Mécanique des fluides" Patrick Chassing

In the case of a flat plate disposed in a uniform parallel flow of a viscous fluid, the thickness of the boundary layer is announced by Blasius as:  

  $$ \delta (x) = 4.92 \frac{x}{\sqrt {Re_x}}$$

such as  Rex is the local Reynolds number $$ Re_x= \frac {x   U(x)}{\nu}$$


                             c- Comparison

We  compare the numerical boundary layer thickness with the Blasius solution.

The following figure shows the dynamic thickness profile versus Blasius solution:


                                      Figure: the dynamic thickness profiles


It can be seen that numerical curve of thickness is similar to Blasius solution. They both have the same order of magnitude ~10-3 m. 



Thermal boundary layer

Numerical simulations                          


       Thermal boundary layer

The thermal study is very important for the further evaporation in the container which represents the main work of this BEI.

We suppose that the incident air flow is at a uniform temperature, and that the surface is maintained in a temperature Tp also uniform but different from the air temperature T0.


        a- Temperature field

We consider now the temperature equations while simulating the problem via OpenFoam. For this, we add the equations in the source code of our own solver as we have already explained.
The following figure shows the evolution of the thermal boundary layer over the flat plate.            
                                 Figure:  the thermal boundary layer over the flat plate
  • By exploring the temperature field perpendicularly to the plate,  we shall observe a progressive variation of Tp to the air temperature T0; in fact this variation is at first fast then it becomes slower and slower as we penetrate into the incident air flow.
  • The region in which T varies in a significant way corresponds the the thermal boundary layer. It is noticeable that the thermal boundary layer gets thicker when we go away from the leading edge: it has the same altitude as the dynamic boundary layer.
  • At the origin, the temperature profile is uniform, then the influence of Tp appears gradually in the fluid. As a consequence, a temperature gradient emerged as well as a heat flux directed from the plate towards the fluid.

So, as well as the dynamic boundary layer is the expression of momentum diffusion, the thermal boundary layer results from the thermal diffusion in the fluid in movement.

              Temperature profile
In the figure below we plot the temperature profile versus the altitude above the plate at x= 0.9 cm.
                  Figure:  the temperature profile versus the altitude above the plate at x= 0.9 cm



      b- Boundary layer thickness

The aim of this part is to estimate the thermal boundary layer thickness values over the flat plate.

The region in which T varies in a significant way corresponds the the thermal boundary layer. However this definition is too vague. The problem was moreover the same with the dynamic boundary layer, and will be approached on the same spirit, namely a conventional definition of the thermal boundary layer thickness.

We define this thickness as following  $$\frac{T(x,\delta_T(x)) - T_p (x)}{T_0 - T_p (x)}= 0.99   $$ .


              Blasius model for a flat plate


The Blasius solution estimates the thermal thickness in such case as follows:

$$ \frac {\delta_T}{x} = 4.64 \sqrt {\frac{a}{U_0  x}}$$ where a is the thermal diffusivity.

       c- Comparison

We  compare now the numerical boundary layer thickness with the Blasius solution.

The following figure shows the thermal thickness profile versus Blasius solution:


                                       Figure: the thermal thickness profiles

As is the case of dynamic boundary layer, the thermal boundary layer thickness has the same profile as the Blasius solution.


           d- Conclusion:

The results correspond perfectly to Blasius solution. Therefore, the study of the boundary layer has enabled us to validate our model in order to use it in the evaporation study of the fuel container.



Project case : Evaporation in the container

Numerical simulations                          



Starting point

Numerical simulations                          


The studied geometry is an airplane part which lies between the wing and the turboreactor. This part is called fuel container.

Figure: Position of the studied geometry in red


It is crossed by the fuel pipes which supply the engine of the reactor. Sometimes, it occurs a leak of fuel in this part of the airplane. The fuel is then gathered in the container and evacuated by its bottom. The fuel may be evaporated there and the concentration of fuel may increases. An air flow is also kept in this structure to make sure that the concentration of fuel does not increase too much.

Even though the  security of the passenger in the airplane is guarantee, Airbus wants to understand the mechanism which drives the evaporation and want to know if CFD is adapted to a such problem.



Airbus, the industrial partner, has given the following geometry

Figure: Initial geometry given by the industrial partner


It's composed by two circular holes at the left hand side situated at the two extremities of the fuel container. Each hole measures 20 mm. To capture the flow inside the geometry, the flow at the inlet and at the outlet has to be correctly solved and involve a large number of cells in this regions.


To avoid meshing problems the inlet and the outlet of the flow have been redefined. Plus the baffles have been enlarged for the same reason avoid meshing problems.
Here below is shown the geometry used for the numerical simulation:


Figure: Simplified geometry


Project case : meshing on IcemCFD

Numerical simulations                        


          IcemCFD meshing


The meshing is one of the critical points of CFD. Numericians are aware that to obtain good simulation results, a clever meshing has to be done.
Enseeiht's teacher used to say that  meshing represents 90% of the time dedicated to the simulation.

The geometry is meshed with tetrahedon and cells has been refined near the wall to capture the wall shear stress. The mesh is composed by 306814 elements.
The following two figures show the mesh. The first gives an overview of the general mesh and the second one is a cut plane of the mesh. It illustrates the mesh refinement near to the wall. 


Figure: Overview of the mesh


Figure: Cut plane of the mesh with refinement near the wall


Then, when the mesh is created, the mesh has to be qualify. And to do so, three different meshes are going to be tested:

Mesh number Total number of elements
1 31 013
2 115 142
3 306 814


Project case : meshing sensitivity

Numerical simulations                        


          meshing sensitivity

The mesh has to be qualified. The results must be independent from the mesh, so the results can be valid and used for further analysis.

Three different meshes have be chosen to the "mesh convergence". Here follows the total number of elements for each one.

Mesh number Total number of elements
1 31 013
2 115 142
3 306 814


The most interessant regions are the ones where the fuel is more likely to evaporate: the region, separated  between the two baffles, which defined one part of the geometry's bottom.
The following picture shows the geometry, the velocity on a slice and the line along where the mesh sensitivity is done.


Figure: Velocity on the slice geometry and line over which the sensitivity will be done


The mesh sensitivity along this line brings up to the light that the mesh used is not the best one. The maximal velocity (Position = 0.1m)  is obviously not the same for the different meshes. Although, the velocity elsewhere along the line seems to be good for both mesh 2 and 3. Mesh 1 can not be used for an analysis because it gives strange results for the maximal velocity. On the other hand, mesh 2 and 3 have a parabolic shape which is what it's expected for a flow between two walls.  Even if these two meshes seem to be good they are not perfect.

Although for a matter of time and CPU capacity, it has been here chosen to go on with mesh number 3.














Figure: mesh sensitivity along the white line


This study underlines the fact that the mesh is not ideal. But we are not pretending to have perfect results. the project's goal is to set up a path which can be used for further analysis.


Project Case: Parametric study

Numerical simulations                          


Once the mesh sensitivity has been done, the parametric study can be processed. The following figure shows the three flight conditions used. For each flight condition, the industrial partner has given the Reynolds number and the density.


We note that the following results correspond to the "lift off" flight's conditions. The next figure shows a screenshot of the temperature field in the container:


                       Figure: The temperature field in the container


        1-  $ \delta_T $ and $ h_c$

The first step to estimate the evaporation rate  $W_v$ is to compute $ \delta_T $ and $ h_c$. For this, the simulations gave us  $ \delta_T $ and $ h_c$. The figure below explains how we have done to extract these values from the simulation results:


                   Figure:  numerical estimation of $ \delta_T $ and $ h_c$ 


         2- Evaporation rate $W_v$

We remind that we have found the expression of the evaporation rate $W_v$ on account density, the boundary layer thickness  $ \delta_T $ and $ h_c$

$$W_v(x)=\frac{b\lambda}{C_p} ln\left( 1+\frac{W_v(x)}{b}\frac{C_p}{h_c(x)}\right)\frac{1}{\delta_T(x)-\delta}$$


Then with this formula established previously, $W_v$ is processed.

To better understand the distribution of the evaporation model, here there is the following figure:

                                   Figure:  distribution of the evaporation model


Otherwise, the following table gave the values of the evaporation rate in different flight cases:


It's important to note that the most accurate result is the one of the lift off flight case because all the simulations were done in this case.

Otherwise, the mass vapor rate seems to be very high. What can be interesting to validate this approach is to compare it to available experimental data but it doesn't exist yet.



Conclusion & Perspectives

Conclusion & perspectives                      


This project has been established in three parts:

First, we have determined the equations describing the evaporation in the container, we have estimated all the non dimensional numbers that appears in this study (Reynolds number Re,  Prandtl number Pr ...). Then we have concluded the expression of the evaporation rate of the fluid in the container.

The second part, we have conceived our own solver "MyIcoFoamT" by developing the standard solver "IcoFoam" under OpenFoam in order to simulate the problem numerically. We have validated this solver in a simple case: The evolution of the boundary layer in a flat plate. In fact, we have compared the profiles given by our solver with the theoretical Blasius model. The results were satisfactory.

Then, we have simulated the principal case of this project which is the evaporation in a container. For this, we have first created the geometry and meshed it with IcemCFD, we have then imported it to OpenFoam to can resolve the evaporation problem with our solver "MyIcoFoamT".

As consequence, we have had good results. Actually, we have succeeded to simulate the evolution of temperature field in the container in only one  specific flight case which is the "Lift off" because  of lack of time. However, we have computed the boundary layer thickness, the convective coefficient and then the rate of evaporation in all flight cases.


So, it is recommended that in further studies, we can simulate the evaporation in the other flight cases. Moreover, what can be interesting to validate the model we have produced is to compare it to available experimental data.