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: http://www.openfoam.com/features/
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.
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 next figure shows the domain at initial instant of simulation
Figure: the domain at initial instant of simulation
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.
Numerical simulations
1- Dynamic boundary layer
a- Velocity field
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:
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 Re_{x} 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.
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 T_{p} also uniform but different from the air temperature T_{0}.
a- Temperature field
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.
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
Source:http://help.solidworks.com/2013/french/SolidWorks/cworks/c_Convection_Heat_Coefficient.htm?format=P
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.
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.
Numerical simulations
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
source: http://photos.linternaute.com/photo/1457769/1553166946/2092/elements-d-avion-de-ligne-boeing-757/
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.
Geometry
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
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 |
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.
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.