Problem modeling

Problem modeling                             

---

 

       1- Lists of symbols & subscripts

               Symbols:

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 ($J.kg^{-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.

Equations

Equations                                                

---

 

       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:

$$Re_{in}=\frac{\rho_{in}LU_{in}}{\mu}$$

• 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.

$$Sc=\frac{\frac{\mu}{\rho_{in}}}{D}$$

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 :

$$\bar{Nu}=\frac{\bar{h_c}*D_h}{k}=\frac{q_t*D_h}{\lambda(T_{in}-T_{S})}$$

 

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}$.

 

     N.B:

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.

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:

 

 

 

 

 

 

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 m²/s. Therefore, the Reynolds number is approximately 327.15. Thus, the test case is a laminar flow.
• The plate's temperature is equal to T_p=90°C while the air temperature is T₀=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.

 

 

 

 

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

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

 

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.