Cfl Courant

Friedrichs Lewy

Espace Technologique

Immeuble Discovery

CFD Computational fluid dynamics

Keywords: Fluids mechanics, solidification, numerical simulation, porous media, volume averaging Bibliographie

I would also like to thank my research supervisors, Professor Dominique Gobin and Doctor Hervé Duval for their contributions in this research project, each one in his own field has given me the scientific guidance much needed for the progress of our work. Special thank to Dominique Gobin who invested his time, pushed my limits with his encouragement and supported me with positive energy.

I would also like to extend my thanks to the jury members; reviewers and invited professors.

First I would like to send my gratitude to Dr. Hervé Combeau and Dr. Frédéric Topin for accepting to review my thesis and permitting me to defend my work. Second, I would like to thank Professor Eric Arquis for accepting to join and president the jury and for reading my work. Finally, I would like to thank Mr. Yves Gaillard for accepting to participate as an invited member and for his help and valuable information for the treatment of the metal specimen.

My special thanks are extended to the Laboratory stuff and my office colleagues that changed every now and then, but even though we managed to create a special bound that I wish will last with time. So, special thanks to Fabian , Sepideh, Carmen and David that already left the laboratory while ago. And, best wishes for Moustafa for his research project.

A special thank to my family. Words cannot express how grateful I am to my mother and father for all the sacrifices that they made for me to achieve my dream. The distance separated us for so long but you were always with me, making me stronger. Special thank to my two sisters Eva and Sabine and their families. And last but not least, special thank to my other half Rami Abou-Eid who is also working on his thesis and managed to be there for me and helped me deal with my stress.

Introduction

The motivation of the present study is the elaboration of metallic foams. A metal foam is a porous structure whose solid matrix has a large fraction of void cells (see Figure 1.1). It has a variety of applications depending on whether the cells are connected or not; e.g. opencells foams are used in heat exchangers to increase heat transfer, while closed-cells foams are employed as impact-absorbing material such as in vehicle's crash box. The competitiveness in the metal industry leads to severe product specifications in order to reduce the cost of the manufacturing process. Thus, metal foams present a promising material since it keeps the high mechanical properties of the metal while reducing the weight up to 90%. Several patented manufacturing processes have been introduced during the last two decades [1,4,[START_REF] Miyoshi | ALPORAS Aluminum Foam: Production Process, Properties, and Applications[END_REF]. However, the metallic melts are foamed either by injecting gas into liquid or by the use of blowing agents.

As a result, the structure of the foam, its homogeneity and its effective properties cannot be rigorously controlled.

The CTIF [3] proposed the casting process CastFoam [2] for metal foams manufacturing which produces well defined homogeneous foams, making use in particular of Kelvin cells to elaborate perfectly regular foams (see Figure 1.2 ). It consists in injecting liquid metal into a porous mould consisting of a regular arrangement of sand or salt beads which are destroyed after solidification of the metallic structure (see Figure 1.3). However, this process needs to be optimized before being employed in the metal industry. Thus, the infiltration and solidification of liquid metal inside the porous mould must be studied. Metal foams may be made of steel, copper, nickel and aluminium. However, the latter is the most widely used in industrial processes because of its low density and its good mechanical properties. Moreover, its fusion temperature is much lower than the other materials, which makes the manufacturing process less expensive and easier to handle. Therefore, liquid aluminium is chosen in our numerical study as well as in the validation experiments.

Our objective is to develop a model able to simulate the infiltration and solidification process inside a porous medium, in order to study the influence of the operating parameters of this manufacturing technique and to monitor the production of homogeneous open cell metal foams.

The metal foams are characterized by their high permeability and porosity values; hence the size of the mould pore is small and can be modelled as a capillary tube. Our modelling approach consists of two major steps. First, the infiltration and solidification mechanisms are studied at the pore level, namely in a capillary tube cooled by its lateral wall. In a second step, a macroscopic model will be developed using the volume averaging method, to simulate the process at the The modelling of the infiltration and solidification inside a capillary mould includes the treatment of four phases: the solid and liquid phases of the metal, the mould and the air. The model is developed based on a one domain approach where the conservation equations of an equivalent domain are solved. The metal phase change is accounted for using the enthalpy-porosity method where penalty terms were introduced to the momentum conservation equations in order to obtain a zero velocity field in the solidified metal and the mould. The metal/air interface shape and position are determined using the VOF method, which introduces some numerical errors known as spurious currents. To minimize the effect of these currents the calculation time step has to fulfil a given criterion which slow down drastically the numerical calculations. Another faced difficulty was the calculation of the effective thermo-physical properties, where at the mould/metal interface the properties were miscalculated on some calculation cells along the mould length.

Developing a macroscopic model of infiltration and solidification of liquid metal inside a porous medium would imply four phases: the solid and liquid phases of the metal, the mould and the air. The derivation of this model using the volumetric averaging method consists in upscaling the local conservation equations, in this case the mass, energy and momentum conservation equations. However, the two interfaces evolves as function of time: the solid/liquid metal interface where the latent heat is exchange, and the metal/air interface where there is a pressure jump caused by the surface tension of the liquid metal. The development of such model is challenging on both mathematical and numerical levels. Therefore, in the scope of this work, the problem is simplified by considering that the mould is initially completely filled by liquid metal at rest, UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE thus the solidification is by pure diffusion.

First in chapter 2, the state of art of the infiltration and phase change at both local and macroscopic scale is presented. In chapter 3 the filling and cooling process at the local level is studied.

The macroscopic model of diffusive solidification of a molten metal inside a porous media considering a non local thermal equilibrium between the solid matrix and the metal is the object of chapter 4. Finally, chapter 5 contains the main conclusions drawn from this study. In the frame of the modelling approach of the complete process of infiltration of the porous matrix by a liquid metal and its subsequent solidification, it is first relevant to study the heat transfer mechanism at the local (pore) scale. The bibliographical review reveals that a huge corpus of publications exists on this topic. Mainly two viewpoints are concerned with such studies:

• on one hand, metallurgists are interested by determining the "fluidity" of metals in capillaries. The development of fluidity tests has been particularly active in the 60s and 70s,

• on the other hand, the community of heat transfer and fluid mechanics has for a long time been concerned by the coupling between solid-liquid phase change and convection in the fluid phase, and in particular by freezing of a fluid flow in ducts.

We will review the main studies in both domains and then try to locate our contribution in this field.

Metallurgy studies on the infiltration of liquid metals in casting

In metallurgy, the efficiency of a number of casting processes such as mould filling depend on the ability of liquid metals or alloys to fill the mould cavity before the metal completely freezes.

For this purpose the concept of "fluidity" has been defined: many experimental tests have been realized consisting in pouring liquid metal by gravity in pyrex tubes and measuring the length of the tube filled with metal when solidification has blocked the fluid flow. The fluidity is the distance in millimetres reached by the metal. This parameter is not a physical property of the metal in question, and obviously it depends on the metal properties (latent heat, viscosity, etc.), but it is also largely resulting from the operational conditions of the test (initial liquid temperature, mould nature, shape and temperature, etc.).

The extreme dispersion of the results for a given metal or alloy is thus not a surprise, and the observations have led in general to qualitative trends and conclusions, as shown in the reviews UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE by Campbell [START_REF] Campbell | Review of fluidity concepts in castings[END_REF] or Ravi et al. [START_REF] Ravi | Fluidity of aluminum alloys and composites: A review[END_REF] where a simple linearity relation was found between the metal superheat and its fluidity. Phenomenological observations have been performed (Niesse et al. [START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF] ; Pai and Jones [START_REF] Pai | Effect of experimental variables on casting fluidity and fluid life of liquid tin[END_REF]) and led to a characterization of solidification modes both for pure metals and binary alloys.

The metal fluidity varies with hydrodynamic and thermal conditions such as: the metal composition, the metallostatic pressure, the metal superheat and the mould preheat. In casting, this property is measured experimentally by fluidity tests. Two main methods are used in these tests, spiral test and vacuum test [START_REF] Sabatino | A review on the fluidity of Al based alloys[END_REF] (see Figure 2.1). In the first test, the metal is poured inside a spiral shape mould, and the fluidity corresponds to the flown distance. The second test consists in suction of the metal using a vacuum pump from a filled recipient.

Solidification modes

After determining the fluidity of the metal, the specimen recovered from the fluidity test can be treated and observed using an optical lens in order to determine the solidification mode.

Solidification of pure metals

The solidification of pure metals is characterized by a distinct smooth liquid-solid interface known as "plane front" [START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF][START_REF] Voller | Fixed grid techniques for phase change problems: A review[END_REF]. However the solidification mechanism depends on other factors such as the metal composition, or the metal initial superheat. If a pure metal is poured in a fluidity test at its fusion temperature, then the solidification starts instantly at the metal-mould wall. As the metal continues to flow the solid phase continues to form till the mould section is completely blocked by the formed solid and flow ceases. Since solidification started first at the mould entrance it will be completed first at that section (see Figure 2.2 ). On the other hand, if the metal is initially superheated, it will flow in the mould till it looses its superheat and then the solidification starts at the metal-mould wall, near the metal-gas interface. The mechanism of solidification is the same as before but the blockage section is moved downstream near the metal-gas interface region where the solidification first started (see Figure 2.3).

Solidification of alloys

Adding alloys elements to a pure metal changes its solidification mode and affects its fluidity. As stated before, the solidification of pure metal is characterized by a sharp solid-liquid interface.

However, when impurities are found in the metal or small quantities of elements are added to obtain an alloy, this smooth interface is broken and dendritic grains are formed. The cessation of the flow is similar to the case of pure metal with a difference that the section is not completely blocked when the flow is stopped. As a result, the fluidity of a dilute alloy is typically lower than the fluidity of the corresponding pure metal. Niesse et al. [START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF] state that "as little as 15 -20% solid in the region where flow is chocked can be sufficient to stop flow in alloys ". ing" occurs at the tube entrance even though the section is not completely solidified [START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF] When the solute quantity is more important, nucleation starts at the tip of the flow and fine equiaxed grains are formed changing the solidification mode of the alloy from dendritic to equiaxed dendritic. Beside the grains shape and size, another difference between the front plane/columnar solidification and equiaxed dendritic solidification is the stopping mechanism of the flow. While in the first case the cessation of the flow is due to the obsturation of the channel section, in the second case the cessation is due to a high friction at the tip of the flow as result of the formation of fine grains which is know as mushy zone. When the solid concentration C superheat [START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF] reaches a critical value the flow is stopped.

Figure 2.6 shows a schematic of the equiaxed dendritic solidification mechanism, while Figure 2.7 shows a longitudinal macrostructures of a rich alloy where it is easy to identify the mushy zone at the tip of the specimen. We can also notice the grains size that increase in the opposite direction of the flow.

Fluidity

As stated earlier, the flown distance in mm of a metal inside a mould is an important characteristic in casting, known as fluidity. This parameter is not a physical property of the metal but the results of experimental work where the thermal and hydrodynamic conditions as well as metal composition are changed. There has been a lot of experimental work conducted in order to measure various casting metal alloys fluidities mainly in the 60 s by M.C. Flemings et al., the effect of several parameters on the metal fluidity were tested especially for aluminium alloys which is convenient for us since we deal with the same metal. However, most of the studies were concerned by the effect of adding elements on the fluidity of alloys [START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF] or the effect of impurities on the fluidity of pure element [START_REF] Ravi | Fluidity of aluminum alloys and composites: A review[END_REF]. Hereby, we will try to review the results of the effect of the metal foam manufacturing process operating conditions on metal fluidity; such as; the metal superheat, the pressure difference (in other word: metal or pressure head), and the mould temperature.

L f = L ρ metal v (f c L f + c p ∆T ) h * (T -T amb ) 1 + B 2 (2.1)
where:

B = h * √ π∆X k s ρ s c p s v (2.2)
L is the mould channel length, v is the velocity of the metal flow, f c is the solid critical fraction, h * is the heat transfer coefficient at the metal-mould interface, ∆X the chocking range in m, k s ρ s and c p s are consecutively the mould thermal conductivity, density and heat capacity. Mould initial temperature

Flemings [START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF] summarized the effect of the mould characteristics on the metal fluidity as shown in Figure 2. [START_REF] Bauer | A general analytical approach toward the thermal conductivity of porous media[END_REF]. As expected, metal fluidity increases when the mould is less conductive (Figure 2.10-a), therefore in the metal foam manufacturing process Ω-sphere sand is chosen as the mould sand because they are characterized by their low conductivity 0.17 W/mK in order to increase metal fluidity and ensure complete mould filling (section 3.2.1). Another interesting observation is the influence of the mould initial temperature on the metal fluidity (Figure 2.10-c), where the latter increases when the mould is preheated which is also expected as the temperature gradient at the mould metal interface decreases, consequently the mould capacity to extract heat from the metal decreases as well. The effect of mould temperature on the metal fluidity using a highly conductive mould is less pronounced as shown in Figure 2. mould temperature the metal fluidity increases as well to a limit when the latter becomes nearly independent of the mould temperature this is because the highly conductive mould (here: cast iron mould) looses its preheat to the ambient. evolution of the system finally may lead to the blockage of the fluid flow in the tube. Given the large number of industrial applications, there is an impressive corpus of scientific publications on the problem of solidification of liquids flowing inside ducts, mainly published in the decades of the 70s and 80s. The bibliography is summarized in the reviews by Epstein and Cheung [START_REF] Epstein | Complex Freezing-Melting Interfaces in Fluid Flow[END_REF], by Fukusako and Yamada [START_REF] Fukusako | Solidification of Pure Liquids and Liquid Mixtures Inside Ducts and Over External Bodies[END_REF] and by Weigand et al. [START_REF] Weigand | Freezing in forced convection flows inside ducts : A review[END_REF]. The review by Weigand et al. [START_REF] Weigand | Freezing in forced convection flows inside ducts : A review[END_REF] is specifically dedicated to freezing processes with forced convection inside ducts, while the papers by Cheung and Epstein [START_REF] Cheung | Solidification and Melting in Fluid Flow[END_REF] review the studies concerning solidification or melting in external or internal flows, including by natural convection. The viewpoint presented by Fukusako and Yamada [START_REF] Fukusako | Solidification of Pure Liquids and Liquid Mixtures Inside Ducts and Over External Bodies[END_REF] is to separately analyse the process of solidification in ducts or around cooled bodies of flowing liquids for pure substances and for binary systems. If we refer to the topic of the present studynamely solidification of an internal forced flow in a ductthe quasi totality of the published works reported in these reviews are concerned with an established flow in a tube with an imposed inlet velocity profile, the temperature of the wall tube being suddenly decreased below the freezing temperature of the fluid (Figure 2.12).

The solid crust formed at the duct wall creates a constriction of the flow, increasing the pressure drop in the pipe until a steady state is reached or obstruction of the duct. Following the early work reported by Zerkle and Sunderland [START_REF] Zerkle | The Effect of Liquid Solidification in a Tube Upon Laminar-Flow Heat Transfer and Pressure Drop[END_REF], most investigations in this domain were interested with the calculation of the evolution of the pressure drop during the solidification process and in characterizing the influence of natural convection. The configuration under study generally consists of a steady, fully developed laminar flow in a rectangular or circular duct whose external wall is suddenly cooled on a given length by imposing a temperature lower than the freezing point of the fluid or extracting a constant heat flux. A steady state situation is generally obtained where the thickness of the frozen crust does not evolve in time and the fluid flows in the As stated before, most of the previous investigations studied the solidification of an isothermal fully developed inner flow, where the temperature at the tube wall is suddenly decreased to T wall < T f usion . However, in most real life applications the scenario is different. In metal casting, the mould is initially empty and the superheated molten fluid initially at rest in a reservoir is poured into it, the metal then solidify by loosing its energy to the cold mould. To our knowledge Cheung and Baker [START_REF] Cheung | Transient Freezing of Liquids in Tube Flow[END_REF] were the first to study the solidification of an infiltrating flow inside a mould. They found an empirical relation between the metal penetration length inside the tube and the flow average velocity based on experimental results. Cheung and Epstein [START_REF] Cheung | Solidification and Melting in Fluid Flow[END_REF] reviewed the work on the solidification of a penetrating flow inside a tube where the penetration length of the metal is calculated as function of (Re, P r, T wall -T f usion ). However in all cases the liquid metal is considered to be initially at its fusion temperature (no superheat) thus the solidification crust closes the section at the tube inlet, and the evolution of the liquid-solid interface radial position is calculated using the analytical solution of a simple cylindrical inward solidification by diffusion, neglecting the effect of the convection in a flowing flow and the axial conduction.

In a more recent work, Seeniraj and Hari [START_REF] Seeniraj | Transient freezing of liquids in forced flow inside convectively cooled tubes[END_REF] studied the solidification of a flowing liquid inside a tube by convective heat transfer with a surrounding coolant. They used the same assumptions and simplifications as the former works cited above (steady state, fully developed flow, axial conduction neglected ...) but instead of decreasing the tube wall temperature lower than the metal fusion temperature, they considered a heat transfer coefficient between the tube wall and a surrounding coolant. They also introduced an interesting solidification commencement condition, that is the solidification will starts only if the energy extracted from the metal is greater that UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE the one added by the fluid flow, this condition is written using the non-dimensional numbers:

Bi Ste > c p 2k N u (2.3)
where c p and k are the metal heat capacity and thermal conductivity consecutively. In the previous presented work, the determination of the unknown variables of the solidification of a fluid flow as the solid-liquid interface position or the penetration length, is done analytically by applying energy balance and neglecting the axial conduction and considering steady state scenario, or empirically by fitting data to experimental works. However, a more accurate solution will be to resolve numerically the mass, momentum and energy conservation equations taking into account all the coupled phenomena of solidification by diffusion and convection. Hibbert et al.

[40] solved numerically the conservation equations of the solidification of a fluid flow using the finite volume method, and they account for the phase change using the enthalpy-porosity method (section 3. 3.3.2). In a more recent numerical work, Mahdaoui et al.

[50] also used the enthalpy-porosity method and studied the effect of thermal, hydrodynamic and chemical parameters (melt temperature, wall temperature, inlet maximum velocity, initial concentration) on the solidification rate of a binary mixture in a laminar fully developed flow.

Literature on the solidification of a penetrating fluid showed some lacks in coupling the thermal and hydrodynamic conditions and normally some strong assumptions were taken (steady state, fully developed flow, neglecting axial conduction ...) and often the thermal and hydrodynamic conditions are different from those employed in metal casting. Actually, the presented reviews showed that the treated flows are fully developed where a velocity profile is imposed at the tube inlet which does not represent the filling process of an empty mould by gravity. Plus the metal flows in a tube where the temperature of the tube wall is suddenly decreased to a temperature lower than its fusion temperature except for a few studies where a heat transfer coefficient is considered between the metal and the ambient or external coolant, where in metallurgy the metal solidifies by loosing its energy to the mould. Our objective in the present study is to develop a model able to simulate both the infiltration and solidification process of a molten metal inside a mould in order to analyse the influence of the operating parameters of the metal foam manufacturing process via casting.

At the macro scale

Infiltration of a porous medium

There is an extensive literature on flow inside a porous medium since it is an important process of interest in many research domains such as soil science, chemical engineering. But to our knowledge there are few studies concerning the infiltration and solidification of a porous medium by a molten metal. However, this subject is now gaining more and more interest especially in material sciences since infiltration processes are used in the manufacturing process of polymer, metal and ceramic matrix composites. The infiltration process consists in replacing one fluid by another, generally gas is pushed out by a flowing liquid, which applies to our case. The liquid flows instantly or demands some applied pressure depending on the fluid nature, wetting or non wetting. In case of a wetting fluid, the driving force of the flow is the Laplace force at the fluid-gas interface. On the other hand, for a non wetting fluid, which is the case of most molten metals, the Laplace force works against the flow and the fluid needs an applied pressure in order to infiltrate the porous medium.

For a two-phase flow inside a porous medium it was found that the capillary number and the dynamic viscosity ratio affect the behaviour and the nature of the flow. Lenormand et al. [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF] characterized three interface displacement patterns: capillary fingering, viscous fingering, and stable displacement. For capillary fingering the liquid flow is characterised by a finger like grows in all direction even backwards, on the other hand for viscous fingering the liquid grows in a tree like shape, both type of fingering are characterised by a low saturation coefficient 0.4 -0.7.

For stable displacement the liquid advances in an almost flat front and the medium saturation is the closest to 1. Huang et al. [START_REF] Huang | Study of immiscible displacements in porous media using a color-gradient-based multiphase lattice Boltzmann method[END_REF] used a Lattice Boltzamnn method to study the displacement of two phase flow inside a porous medium, 17 cases were simulated were the capillary number and cinematic viscous ratio (M = µ nw µ w ) were varied in order to cover the viscous and capillary fingering scenarios. It was found that for a large capillary number the displacement is stable, while for a low capillary number the flow nature becomes capillary fingering. On the other hand, for a viscosity ratio lower or close to 1 the flow pattern is viscous fingering, and for an intermediate capillary number and a viscosity ratio larger than 1 the displacement is stable. In our case, the capillary number is low (Ca ∼ O(10 -4 )) and the dynamic viscosity ratio is larger then 1 (M ∼ 2). The presented cases do not cover all the possible scenarios therefore using the information above do not permit us to clearly identify the nature of the flow in the mould filling process. However the metal foam obtained from the casting process do not present major defect which means that the filling saturation is close to 1 and the flow pattern is stable displacement.

To take the capillary forces into consideration, the Laplace pressure drop at the infiltration front is given by Young-Laplace equation:

∆p L = σκ cosθ (2.4)
where σ is the liquid metal -gas surface tension, θ is the fluid-solid contact angle, and κ is the curvature of the porous medium. The latter is unknown, but several models exist based on the porous medium geometry that permit the estimation of this parameter. For instance Mortensen et al. [START_REF] Mortensen | On the Infiltration of Metal Matrix Composites[END_REF] calculated the curvature as the fibre surface area per unit volume S f that is

S f = 4V f d f (1 -V f ) (2.
S f = 6V f d part (1 -V f ) (2.6)
for a porous medium formed by a packed spherical particles of a diameter d part .

While Michaud et al. [START_REF] Michaud | Melt infiltration of metal matrix composites[END_REF] calculated the curvature at the liquid-gas interface for a fibre porous medium as:

κ = A v 1 -V f (2.7)
where A v is the specific surface area, and V f is the fibre volume fraction. And they used the slug flow assumption where they introduced the "saturation" parameter S defined as a step function, which is equal 1 in the infiltrated part of the porous medium and 0 in the part filled with gas.

The extended Darcy's law for a two-phase system using the above parameter is written as:

(1 -V f )Su = - K µ ∇p (2.8)
where K is the permeability of the porous medium defined as K = K r K s . The saturated permeability K s or intrinsic permeability [START_REF] Helmig | Efficient modelling of flow and transport in porous media using multi-physics and multi-scale approaches[END_REF] depends on the porous medium geometry and several simple models already exist as Kozeny and Kozeny-Carman models [START_REF] Carman | Fluid flow through granular beds[END_REF][START_REF] Kozeny | Über kapillare leitung der wassers im boden[END_REF], in other words the intrinsic permeability is the permeability of a medium saturated by one phase. While, the relative permeability K r represents the difference of a two or more fluids permeability inside a porous medium compared to a one phase flow, thus this property depends on the saturation S defined earlier, and is defined empirically in soil science as:

K r = S n (2.9)
where n is typically between 1 and 3 [START_REF] Spitz | A Practical Guide to Groundwater and Solute Transport Modeling[END_REF].

The Laplace pressure drop at the interface between two phases can also be directly given as function of the saturation and the entry pressure p d which is the minimum pressure required for a non wetting fluid to infiltrate a porous medium. Such as the Brook-Corey model [START_REF] Helmig | Efficient modelling of flow and transport in porous media using multi-physics and multi-scale approaches[END_REF]:

∆p L = p d S -1/λ (2.10)
Where λ is called pore-size distribution and has a value between 0.2 and 3, determined by fitting experimental data.

Using the slug flow assumption, an analytical solution that gives the position of the liquid-gas interface as a function of time for an isothermal two-phase flow is given by Michaud et al. [START_REF] Michaud | Melt infiltration of metal matrix composites[END_REF]: (2.11) indicates that the infiltration kinetics is parabolic, different from the case of capillary tube infiltration where the infiltrated length varies linearly as function of time.

L = 2K(∆p -∆p L ) µ(1 -V f ) √ t (2.

Solidification in a porous medium

Modelling phase change is a challenging task even in a one phase domain, the main difficulty lies in the presence of two unknown variables in the energy equation, the temperature distribution and the solid-liquid interface position which evolves in time. These two variables are coupled, thus analytical solutions can be found for some simple configurations as we will see in chapter 2. Dealing with phase change becomes more complex in the presence of other phases in a porous medium. Therefore, the existing studies tend to simplify the problem by taking several assumptions as considering the effective properties invariant in temperature, time and space and to neglect convection. But the most important simplification is to consider thermal equilibrium between phases so that one energy equation is solved for the whole domain.

Rubin et al. [START_REF] Rubin | Heat Transfer in Porous Meida with Phase Change*[END_REF] modelled the liquid/vapour phase change inside a porous medium. They proposed a continuum model using Darcy law, but they did not take into consideration the effect of surface tension and they neglected natural convection. As most of the phase change studies, local thermal equilibrium between the solid matrix and fluid was considered. For a steady state assumption and constant effective properties, an analytical solution of the temperature distribution is proposed for two boundary conditions: fixed temperature, and constant heat flux.

To our knowledge, Beckermann et al. [START_REF] Beckermann | Natural convection solid/liquid phase change in porous media[END_REF] were the first to develop a solid-liquid phase change continuum model and solve it numerically. They considered natural convection in the liquid phase using the Darcy-Forchheimer law and they determined the permeability of the porous media as function of the pore diameter and the medium porosity using Kozeny-Carman model. However, they assumed a local thermal equilibrium between the metallic mixture and the solid matrix, which led to one energy equation solved in the continuum domain. This assumption simplifies drastically the problem and does not take into consideration cases where the conductivity ratio between the phases is high. Thus, they used an alternative mixture model [START_REF] Veinberg | Permeability, Electrical Conductivity, Dielectric Constant and Thermal Conductivity of a Medium with Spherical and Ellipsoidal Inclusions[END_REF] to calculate the effective conductivity of the metallic mixture/ solid matrix. Comparison between the numerical results and experiments showed some differences which may be the result of the LTE assumption.

Zhang et al. [START_REF] Zhang | Analytical solution for solidification of close-celled metal foams[END_REF] proposed an analytical solution of the liquid-solid interface position as function of time by extending the Neuman's analytical solution of the Stefan problem for two semiinfinite walls. They used an existing analytical model based on thermal equilibrium between the phases [START_REF] Bauer | A general analytical approach toward the thermal conductivity of porous media[END_REF] and a shape factor parameter β s to determine the effective thermal conductivity of a porous media which was considered invariant in temperature and time. where τ = α th t/L 2 is the non-dimensional time.

Recently, advances in the numerical modelling of phase change in a porous media have been obtained considering non local thermal equilibrium (LTNE) between the solid matrix and the fluid. Harris et al. [START_REF] Harris | Phase-change phenomena in porous media -a non-local thermal equilibrium model[END_REF] studied conduction in a porous media in the presence of phase change,

where LTNE between the solid matrix and fluid is considered. For a semi-infinite medium, and quasi steady-state assumption an analytical model is presented based on lumped heat capacitance method, where the energy equation is linearised by calculating the mean heat capacity.

Three parameters define the model; the mean heat capacity, the heat resistance in the solid matrix and the interstitial heat transfer coefficient that is evaluated for a square prism pore shape [START_REF] Minkowycz | On departure from local thermal equilibrium in porous media due to a rapidly changing heat source: The Sparrow number[END_REF]. Liu et al. [START_REF] Liu | Numerical modeling for solid-liquid phase change phenomena in porous media: Shell-and-tube type latent heat thermal energy storage[END_REF] modelled melting of a PCM in a porous medium with the presence of natural convection using the enthalpy-porosity method [START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF]. To solve the momentum equation, the permeability and inertia coefficient were determined using Calmidi model [START_REF] Calmidi | Forced Convection in High Porosity Metal Foams[END_REF]. They considered a non local thermal equilibrium between the solid matrix and the fluid, and wrote two energy equations coupled by a heat transfer coefficient calculated using an existing correlation [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a horizontal cylinder[END_REF] for forced convection in a horizontal cylinder. Three structured model (tetrakaidecahedron) was used for the determination of the effective conductivity in each zone (PCM and solid matrix). Comparison between experimental results of PCM melt inside a porous medium and the numerical results showed that the non local thermal equilibrium model was more accurate compared to local thermal equilibrium (LTE) model. In the same spirit, Alomar et al. [START_REF] Alomar | Simulation of complete liquidvapour phase change process inside porous evaporator using local thermal non-equilibrium model[END_REF] simulate liquid-vapour phase change under steady flow condition and considered non local thermal equilibrium between the fluid and the solid matrix. But they used simple mixture model and calculate the effective conductivity as the arithmetic mean value of the phases conductivities, which in our opinion, fails to represent cases where the conductivity ratio between the phases is high. They also coupled the two energy equations by a heat convection coefficient estimated using a model, where h = f (k, d p , P r, Re). To solve the momentum equation, the permeability of the porous medium was calculated using Carmen-Kozeny model for a packed bed of spheres.

They compared the numerical simulations of the water phase change inside an evaporator using the LTE and LTNE models. They found that there is hardly a difference between the temperature profiles and this is due to the high convective heat transfer used to couple the two energy equations. Therefore, this comparison does not prove the utility of the LTNE model.

Literature on solidification in a porous medium showed some lack in dealing properly with the problem. Existing studies are normally based on a LTE between the fluid and the solid matrix assumption. Which means that the solid matrix wall is considered to be at the fusion temperature of the solidifying metal. This is a strong hypothesis, that is usually unrealistic especially in metal foam casting where the fusion temperature is higher than the burning temperature of the sand mould. equations using a correlated heat transfer coefficient, which we believe, is less accurate than deducing it from the averaging of the energy equations at the local scale. Therefore, our goal is to derive a macroscopic model of the solidification of a liquid metal inside a porous medium considering LTNE between the metal and the solid matrix. This model is developed using the volume averaging method [START_REF] Whitaker | The Method of Volume Averaging[END_REF] by applying the averaging operator on the conservation equations at the local scale and then resolving the derived closure problems to calculate the effective properties presented in the closed formed of the 2 energy equations at macro scale.

Conclusion

A new metal foam manufacturing process is proposed by the CTIF [3], it consists in filling a sand preform by a molten metal, normally aluminium, copper, nickel and steel. The thermal and hydrodynamic conditions of this process should be controlled based on the properties, size and shape of the metal foam piece. Therefore, the infiltration and solidification of a liquid inside a porous medium must be study. Our approach consists of studying the infiltration and solidification of a liquid metal on a local scale inside a cylindrical tube, then upscaling to the mould scale using the volume average method to obtain a macroscopic model of solidification inside a porous medium.

At the local scale, literature showed that infiltration and solidification of a liquid metal inside a cylindrical tube interested metallurgists and heat transfer community. These two communities dealt differently in tackling with the problem. In metallurgy, studies were mainly experimental where a measured property is defined, that is the metal flown distance in mm called fluidity. This parameter is the fruit of experimental tests where the thermal and hydrodynamic conditions are changed accordingly empirical relations are found. Also, metallurgists were concerned about the solidification modes, observations of the solidified specimens under optical lens showed the differences between the solidification modes of pure metals and alloys. In heat transfer studies, the solidification of a fully developed flow inside a tube is treated, and by applying energy balance between metal and the tube relation of the evolution of the solid/liquid interface as function of time is found. However, many non representative assumptions and simplifications were taken (steady state, negligible axial conduction, constant velocity profile ...), and the imposed boundary conditions are also different from the metal casting process. In literature, the case of a fully developed flow in a long tube where the temperature is suddenly decreased to a temperature below the fusion temperature is treated. However in metal casting, the fluid is poured into the mould by gravity which is more accurately translated by a fixed pressure difference between the tube inlet and outlet, and the metal freezes by loosing its energy to the mould that have a certain volume instead of a thin wall with a Dirichlet condition of fixed temperature. where they account for the liquid/solid phase change using the enthalpy-porosity method. However, to our knowledge the infiltration problem was never coupled to the phase change inside a tube, therefore the penetration length or "fluidity" is still an unknown. Thus, our objective is to develop a model that treats the infiltration and solidification of a liquid metal and to determine the evolution of the fluidity and the infiltration time as function of the operating conditions (pressure difference, metal superheat, mould initial temperature) based on the numerical results.

At the macro scale, liquid flow inside a porous medium was extensively studied where the porous medium permeability is modelled and well known relation between the flow velocity and pressure are found whether for low Reynolds number (Darcy flow) or high Re (Forchheimer flow). Recently, the flow of 2 immiscible fluids inside a porous medium has been studied, where the flow pattern is determined by two parameter, the capillary number and the kinematic viscosities ratio. New relation between the velocity and pressure gradient is found where the Darcy law is extended taking into account the relative porous medium permeability as function of the medium saturation fraction S ∈ [0, 1]. For an isothermal flow, analytical solution of the infiltration length as function of time is given by Michaud et al. [START_REF] Michaud | Melt infiltration of metal matrix composites[END_REF].

Existing studies of the solidification of a liquid inside a porous medium are based on a local thermal equilibrium (LTE) between the metallic mixture and the mould. This assumption is not valid when treating metal phase change problem where the fusion temperature exceed the mould maximum temperature. Analytical solution of the liquid-solid interface position as function of time is proposed by considering LTE and using a shape factor to transform the existing Neumann solution. Recent studies considered non local thermal equilibrium (LTNE) between the metal mixture and the mould, however the two energy equations were coupled using a correlated heat transfer coefficient which does not represent properly the heat exchange between the two domains (metallic mixture, mould). Therefore, our objective is to model the solidification of liquid metal inside a porous medium by considering LTE between the liquid and solid phases of the metal and LTNE between the metallic mixture and the mould. Infiltration and solidification at local scale

Introduction

The infiltration and solidification of a liquid flowing in a capillary tubes or, by extension, in porous structures is a problem encountered in a great number of industrial applications or natural processes. Many examples may be found in human activities. In this chapter we will focus on modelling the infiltration and solidification of a molten metal inside a capillary tube. If we refer to the bibliography concerning the infiltration of a capillary tube with heat transfer, we can identify two classes of studies. The first one consists in the classical problem of forced convection in cylindrical ducts, while the other class considers infiltration process, where the injected fluid pushes another fluid (namely gas) initially present in the duct.

When the heat transfer problems involves phase change of the injected liquid, the transient character of the mechanism is crucial, since the time evolution of the system finally leads to blocking the fluid flow in the tube. If we refer to the topic of the present study, solidification of an internal forced flow inside a duct, the quasi totality of the published works reported in the reviews presented in the bibliography (section 2. 1.2) are concerned with an established flow in a tube with an imposed inlet velocity profile, where the temperature of the wall tube being suddenly decreased below the freezing temperature of the fluid. The solid crust formed at the duct wall creates a constriction of the flow, increasing the pressure drop in the pipe until a steady state is reached or occlusion of the duct.

The purpose of the work presented in this chapter is to study the solidification dynamics of a liquid metal infiltrating a capillary tube in the perspective of developing a macroscopic model of infiltration and solidification in a porous matrix. Therefore, the problem is characterized by two specific features which are not present in the heat transfer literature on the topic: 1. the fluid flow is due to a given pressure gradient between the duct inlet and outlet, in order to model infiltration by gravity, 2. the duct is embedded in an initially cold mould and the latent heat generated by the solidification process is absorbed by the mould.

These conditions are specified in order to approach the conditions of casting met in the process of elaboration of metal foams. Our main concern will be the determination of the blockage time and the infiltration length (fluidity) under various operational conditions.

UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE
We consider a capillary tube whose internal radius R 0 is small compared to its length L. A liquid metal is injected at the duct entrance and a pressure difference is applied between the inlet and the outlet of the duct. In a first phase of the process, liquid metal infiltrates the tube and pushes the air out. At this stage, the liquid-gas interface is tracked as function of time. Then, the liquid metal looses its latent heat and solidification starts. The liquid metal progressively solidifies at the walls of the mould up to a certain point where the solid phase blocks up the capillary tube and stops the liquid flow. Then, solidification proceeds by heat diffusion in the liquid trapped upstream.

Casting tests

Fluidity tests

The most employed tests to measure metal fluidity (defined in section 2. 1.1 ) are the vacuum fluidity test and the spiral test. CTIF [3] chose the latter to carry out their tests with a spiral having a triangular section of 3 mm length each side (see Figure 3 For one initial superheat value (∆T = 190 K), four aluminium (99.5%) fluidity tests were done from which two in silica and two in Ω -sphere sand moulds. The thermo-physical properties of the two sands type employed as well as the metal properties are presented in Table 3.1 and Table 3.2 consecutively. The metallo-static height used to drive the flow inside the mould is 135 mm, this values is kept constant for each test and during the infiltration process. The pouring temperature is measured at the beginning of each test, this value could not be controlled rigorously, but its fluctuation from a test to another is estimated to be less than 15 degrees. 3.4. Since the pouring temperature was hard to control, we can notice that the results obtained for the "same" hydrodynamic and thermal conditions have a difference of less than 10% for the infiltration length and 30% for the infiltration time. Even though a proper validation of our model against the results of the fluidity tests cannot be done properly due to the incertitude of the experimental measurements and to the difference in the hydrodynamic conditions, we will see in section 3.7.2 that the tendency is the same.

Solidification modes

Before examining the samples recovered from the fluidity tests, a metallo-graphic preparation was made. The preparation consists of two steps: first the samples were coated and polished with an automatic polisher, then macro and micro etching were made on the samples to reveal the grain boundaries and micro-structure which leads to the identification of the solidification modes.

• Metallo-graphic Preparation

The surface of a metallo-graphic specimen is prepared by various methods, mechanical preparation is the most common one. The specimen are polished mechanically using abrasive papers and diamond pastes. Then the micro-structure can be revealed by etching.

Automatic polishing is done with a constant bearing force of 10 N and a 150 revolutions / minute plate speed. The supplier of all the utilised products is Struers.

-Step 1: Pre-polishing plane polishing using a new abrasive paper (MD-Piano 220) of 68µm particle size for 2 minutes to remove the layer of the resin and reveal the metal.

- polishing using a used abrasive paper (MD-Largo) of 9µm particle size with a diamond paste (DiaDuo) for 30 seconds.

-Step 3: Polishing polishing using a used wool abrasive paper (MD-Mol) of 3µm particle size with a diamond paste (DiaDuo 3µm) for 1 minute, followed by an identical polishing with a MD-Mol paper size of 1µm.

-Step 4: Final Polishing polishing using a new velour paper (MD-Nap) of particle size < 1µm with a colloidal silica paste (OP-S) with a particle size∼ 0.04µm for 1 minute.

• Etching technique

Two types of attacks are made: macro and micro attack. For micro attack, that reveals grain micro-structure within the grain, the Dix-Keller reactant was used. This reactant is suitable for pure aluminium whose composition is as follows: 190 ml of water, 5 ml of nitric acid, 10 ml of hydrochloric acid and 2 ml of hydrofluoric acid. This acid is very dangerous and strong, for that all safety precautions must be taken when handling it and an attack of a few seconds is sufficient to reveal the micro-structure. After the etching, observations under the optical microscope must be made.

For the macro attack, which aims to reveal the solidified grains borders, two solutions were used, one of which had no real effect on aluminium, it has a composition of 100 ml of water and 15 g sodium hydroxide. The second macro attack, that was provided by the CTIF is iron chloride. It reveals the shape of the grains to the naked eye, therefore simple observation under the binocular is enough.

Two samples recovered from the fluidity tests of pure aluminium (99.5%) having an initial superheat of 190 K (T = 1123 K) were analysed. The fluidity tests were done in 2 moulds, one in silica sand and the other in Ω -sphere sand. The observation of these samples under the binocular allowed us to identify the solidification modes.

The first sample corresponds to the casting test of pure aluminium poured into a silica mould, having an approximate length of 15 cm. The analysed segments are located at 1.5 cm, 4 cm, 10 cm and 12.5 cm from the entrance of the tube. Shrinkage porosities have prevented us to test the area between 4 and 10 cm, and have affected the solidification mode along the sample.

Thus, we note in Figure 3.3 (a) the presence of equiaxed grains in the shrinkage porosity zone.

Otherwise, we can identify three solidification zones; region A spread to at least 6.5 cm from the entrance, where we have one or two long columnar grains per diameter. Then region B

where the columnar grains point upstream as Flemings et al. noted [START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF] for solidification of a The second sample corresponds to the aluminium cast into a mould of Ω -sphere sand having an approximate length of 30 cm. Thus, we note that for the same hydrodynamic and thermal conditions, the fluidity of the metal is twice larger when using Ω -sphere sand, since the latter is characterise by its low conductivity which is translated by a low heat transfer between the isolator mould and the metal. The analysed segments are located at 21 cm, 23.5 cm, 26 cm and 28.5 cm from the entrance of the tube. We can identify three solidification zones in Figure 3.4. As the analysed segments were taken from the final part of the sample, region A presented earlier in the analysis of specimen 1 is not presented here. But, we can identify region B where columnar grains grow in the opposite direction of the heat flow (from the mould inner surface to the tube axis), and where equiaxed grains are formed at the tube axis. Followed by region C (that was not identified in specimen 1) where the equiaxed grains multiply, and finally the region D where equiaxed grains are finer.

Figures 3.3 & 3.4 show that there are no difference in the solidification modes of the pure aluminium (99.5%) depending on the sand mould, silica or Ω -sphere. The main difference is the metal fluidity that is much more important in Ω -sphere mould, therefore this mould is chosen for the casting of metal foams and its properties are used in our mathematical model. The comparison of the macro-stucture of the aluminium (99.5%) with the bibliography presented earlier (section 2. 1.1.1) shows that the solidification mechanisms corresponds to the solidification of an impure metal or an alloy with small quantity of solute (Figure 2.4). This is probably because the industrial pure aluminium 99.5% is not pure enough to obtain a planar front solidification, but for simplification reasons we will consider that the metal is pure and that the solid-liquid interface is sharp.

UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Mathematical model

Flow of two immiscible fluids

The problem of liquid metal hydrodynamics in the infiltration of a capillary tube initially filled with air involves three phases: solid, liquid and gas. The flow is driven by a pressure difference between the tube inlet and outlet, thus the flow velocity depends on the liquid-air-solid contact angle and the surface tension. Therefore we must estimate the effect of the Laplace force. An exact value of the contact angle θ between the liquid aluminium and the sand mould could not be found in the literature, and would require a separate experimental study. However, literature gives an approximation of this angle for other materials: for liquid aluminium on alumina at fusion temperature the contact angle at equilibrium is θ e 160° [START_REF] Mortensen | Infiltration of fibrous preforms by a pure metal: Part III. capillary phenomena[END_REF], for liquid aluminium at T = 1173 K on ceramic material (SiC) θ e 150° [START_REF] Hashim | The wettability of SiC particles by molten aluminium alloy[END_REF], and for liquid aluminium at fusion 148° [START_REF] Shen | Reactive Wetting of SiO 2 Substrates by Molten Al[END_REF]. We can conclude that in general for T f usion < T < 1173 K liquid aluminium is a non-wetting fluid, and that the applied pressure should overtake the Laplace pressure to allow for the fluid to flow inside the tube. The capillary number is given by:

Ca = µv σ ∼ O(10 -4 ) (3.1)
where µ is the dynamic viscosity, σ is the surface tension and v is the velocity. As the expression above shows the dimensionless capillary number represents the ratio between viscous and surface tension forces. For a small capillary number (≤ 10 -5 ) the capillary forces are important and for a high capillary number (≥ 10 -3 ) the capillary forces are negligible compared to the viscous forces. Therefore in our case the capillary forces cannot be neglected due to the high surface tension of most of liquid metals e.g. for liquid aluminium σ = 0.85 N/m [START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF], and we need to focus on the deformation and transport of the liquid-air interface.

In general, there are two approaches to identify the interface in multiphase flow: moving and fixed grid techniques. The latter itself contains two groups of methods, interface tracking method and interface capturing method. The interface tracking method is explicit and the interface is delimited by markers that are advected by the flow, this requires re-meshing techniques that slow down the numerical simulation and has its limitation when the interface has a complicated shape. On the other hand, the interface capturing method is an implicit method that employs the transport equation of an indicator which may be a volume fraction [START_REF] Hirt | Volume of fluid ( vof ) method for the dynamics of free boundaries *[END_REF][START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF], a level set function [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF][START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] or a phase field variable [START_REF] Jacqmin | Calculation of two-phase navier-stokes flows using phase-field modeling[END_REF], where the evolution of the interface is determined from the indicator solution as function of time.

Several methods are based on the interface capturing approach, the most common are the Level Set method (LS) [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] , and the Volume Of Fluid method (VOF) [START_REF] Hirt | Volume of fluid ( vof ) method for the dynamics of free boundaries *[END_REF]. In the LS method the interface is assigned by a distance function, this was proven to give a sharp interface solution.

Whereas, the VOF method employs the resolution of a transport equation of a scalar that can take two values (0 or 1), and the interface is located in the transition region. That is the reason why using VOF method can lead to smearing the interface, but this problem was solved by employing a "compression term" [START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF] that will be discussed later on. And as VOF method guarantees much better mass conservation than the other methods [START_REF] Scardovelli | Direct Numerical Simulation of Free-Surface and Interfacial Flow[END_REF], it has been chosen for our study.

Volume Of Fluid method (VOF)

The Volume Of Fluid method was first presented by Hirt and Nichols [START_REF] Hirt | Volume of fluid ( vof ) method for the dynamics of free boundaries *[END_REF]. Since then many researchers worked on developing and improving this method especially in dealing with interface smearing [START_REF] Chen | The development of a bubble rising in a viscous liquid[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Rudman | Volume-Tracking Methods for Interfacial Flow Calculations[END_REF] and it is now widely employed in commercial CFD codes. As mentioned before, the VOF method consists in solving an additional transport equation for a volume fraction simultaneously with the continuity and momentum equations. (usually a gas). The interface between the two fluids is determined by the transition region from 0 to 1. Often, an iso-surface having the value of 0.5 is used to determine the interface between the two fluids. This method is based on a one-domain approach where one set of conservation equations is solved for the whole domain, and the phase properties are calculated as average values as function of the indicators and are discontinued at the interfaces. Hereby, we are going to present the VOF method implemented in OpenFOAM [5] which uses an artificial compression term in the transport equation of the volume fraction that gives a sharp interface solution while ensuring a bounded solution [START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF].

The fluid flow is governed by a modified Navier-Stokes equation for an incompressible fluid that accounts for the pressure jump at the liquid-gas interface, to which a transport equation of the volume fraction α is added.

• Continuity equation ∂ρ ∂t + ∇.(ρ v) = 0 (3.2)
• Momentum conservation equation

∂ρ v ∂t + ∇.(ρ vv) = -∇p + ∇.(µ∇v) + 1 V S(t) σ κ n δ dS (3.3) • Transport equation of the volume fraction α ∂α ∂t + ∇.(vα) + ∇. [v r α(1 -α)] = 0 (3.4)
The last term in eq.( 3.3) is a surface force that acts only at the interface thus the use of the Dirac δ function. This surface force is due to the surface tension σ at the liquid-gas interface where it generates a pressure jump. As the interface capturing methods are implicit, the location/shape of the interface are unknown and thus the integration of the surface tension on the interface S(t) cannot be done. To account for this term in the momentum equation, the surface force is transformed into a continuous volumetric force using the Continuum Surface Force (CFS) model [START_REF] Berberovic | Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution[END_REF][START_REF] Brackbill | A Continuum Method for Modeling Surface Tension[END_REF] where this term is evaluated from the volume fraction distribution: Where n is the normal vector at the interface. The above definition does not take into consideration the contact angle at the contact line between the 2 fluids and the wall, this will be discussed later on in the mathematical model section 3.4. 2.

1 V S(t) σ κ n δ dS = σ κ ∇α (3.
As mentioned before, the convection of a step function produces a diffusive interface which is the consequence of the numerical diffusion due to the discretization scheme of the convective term. Special treatment must be taken in order to obtain a sharp interface solution. On the other hand, using higher-order schemes gives an unbounded solution of the volume fraction.

Many researchers use a compressive differencing schemes [START_REF] Rider | Reconstructing Volume Tracking[END_REF][START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF][START_REF] Ubbink | Numerical prediction of two fluid systems with sharp interfaces[END_REF], but openFOAM uses

Weller scheme [START_REF] Weller | Interface capturing methodology[END_REF][START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF], where an additional "compression term" is employed in the transport equation of the volume fraction that corresponds to the third term of eq.( 3.4). It can be noticed that this artificial term is active only within the interface region where the multiplication factor α(1 -α) = 0. The velocity used to transport this term is taken to be the relative velocity between the two fluids.

Phase change modelling

Similarly to the two-phase flow, the phase change of a material presents an interface between two phases where there is a variable discontinuity. In case of two phase flow, the discontinuity is translated by the pressure jump. But in phase change problem, metal latent heat is exchanged at the solid-liquid interface which corresponds to the heat transfer discontinuity, this quantity should be tracked as function of time. This bring us back to the two modelling methods: fixed grid and front tracking also known in this case as 1 and 2-domain methods simultaneously. In 2-domain method, the conservation equations are solved in each domain separately and the two system of equations are coupled by the Stefan-boundary condition at the liquid-solid interface [START_REF] Sarler | Stefan ' s work on solid-liquid phase changes[END_REF]. This method is usually employed when dealing with pure metal solidification where the solid and liquid phases are separated by a front plane interface as described earlier. On the other hand, in 1-domain method the solid and liquid form a one equivalent domain having effective thermo-physical properties, therefore the continuity equations are solved in this one domain and the phase change is accounted for by various methods from which the enthalpy method that will be presented afterwards. This method can be used for the treatment of all solidification modes but it is mostly used in case of the solidification of alloys because it can represent the mushy solidification zone where the 2-domain methods fails For some simple solidification cases, analytical solutions can be found in literature, to our knowledge the first similarity solution of two semi-infinite domain in contact was proposed by Neumann [START_REF] Carslaw | Conduction of Heat in Solids[END_REF]. It is important to review some of these cases, because analytical solutions give insight of what to expect to find when dealing with more complicated cases such as phase change in 2 phase fluid flow.

Two semi-infinite domains

Considering a semi-infinite liquid metal domain at a temperature higher than its fusion temperature in perfect contact with a semi-infinite mould. The solidification is by pure conduction, the liquid is considered to be motionless. The position of the solid-liquid interface and the temperature distribution in the 3 domains (liquid/solid metal and the mould) as function of time can be found by solving a transcendental equation [START_REF] Dantzig | Solidification[END_REF]. The solution shows that the solid grows proportionally to √ t, and that the solid-liquid interface velocity is proportional to 1/ √ t.

Axial or radial symmetric processes

Using the quasi-stationary approximation Alexiades & Solomon [START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF] proposed a solution for the temperature distribution in both liquid and solid phases and the solid-liquid interface position as function of time for the solidification of a liquid metal at its fusion temperature (no initial superheat) in a cylinder or sphere. Two boundary conditions are treated: imposed cold temperature below the metal fusion temperature (Dirrichlet condition), and imposed constant negative thermal flux (Neumann condition) at the cylinder/sphere wall. The solutions of these simple configurations are very useful as they will permit us later on (section 3.5.2) to validate our model for one simple case of solidification by pure conduction of liquid metal with no superheat in a cylinder.

Enthalpy method

Modelling the solidification process using the fixed grid techniques enable us to have the freedom of choosing between pure metal or alloy since it can deal with the solidification of both materials. And for the same reasons that drove us to choose the VOF method to track the liquidgas interface, we preferred the fixed grid technique over the front tracking method to monitor the solid-liquid interface. In the fixed grid method the Stefan problem is not presented explicitly, instead the energy released at the solid-liquid interface is accounted for implicitly in the energy conservation equation. This can be done using several techniques, such as: apparent heat capacity, latent heat source term, total enthalpy. However, these methods are usually used in case of phase change by pure diffusion where the energy equation is not coupled to the momentum equation. To model the convective-diffusion phase change we choose the enthalpy-porosity method presented by Voller et al. [START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF] where the zero velocity condition of the solidified metal is accounted for using a Darcy source term in the momentum equation by considering the solidliquid domain as a porous medium.

The enthalpy-porosity method consists in writing the energy conservation equation in term of the enthalpy:

∂(ρh) ∂t + ∇.(ρh v) = ∇.(k∇T ) (3.7)
The enthalpy h is equal to the sensible heat plus the latent heat released when the material change the phase from liquid to solid. By making the assumption of constant thermo-physical properties between the solid and liquid phases, the enthalpy is written as:

h = h s + L f (T ) = T T ref c p dθ + g l L f (3.8)
where g l represents the liquid fraction and is function of the temperature. The relationship between the liquid fraction and the temperature can be simple as a linear variation or step function in case of isothermal solidification of pure metal, but it can also be complex if one choose to model the equiaxed solidification [START_REF] Rappaz | Modelling of microstructure formation in solidification processes no[END_REF]. Re-writing the energy equation as function of one variable, the temperature gives:

∂(ρc p T ) ∂t + ∇.(ρc p T v) = ∇.(k∇T ) + S th (3.9)
where S th is the thermal source term that accounts for the latent heat in the energy equation and it is written as function of the liquid fraction:

S th = -ρ L f ∂g l ∂t + ∇.(g l v) (3.10)
However, the solidification of pure metals occurs at the fusion temperature. Thus, the liquid volume fraction behaves as a Heaviside step function, having a Dirac delta function as divergence. As a result, the divergence term in eq. ( 3.10) has a non-null value only at the solid/liquid interface. However, the velocity at the solidification front is null therefore for solidification of pure metals the second term of the RHS of eq. (3.10) can be neglected. This is not the case when dealing with the solidification of alloys [START_REF] Voller | A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[END_REF].

As mentioned earlier in the enthalpy-porosity method a penalization term is added to the mo-

Infiltration and solidification of liquid metal

The heat transfer and fluid flow process of a liquid metal inside a capillary tube consists of three problems: invasion of the tube by the liquid metal and air-metal interface displacement, liquid metal flow and phase change at the mould-metal interface and heat conduction in the mould.

To simplify the problem we will assume that the fluid flow is laminar and incompressible, both in the liquid metal and in the air initially present in the tube. This hypothesis is justified for a molten metal in the frame of our study where the Reynolds number is less than 2300 and the metal density is considered to be invariant as function of temperature. We also considered that the natural convection inside the tube is negligible compared to forced convection. The non dimensional Richardson number (Ri) permits the validation of this assumption, it represents the ratio between Grashof (Gr) number and Reynolds (Re) number powered by 2. For liquid aluminium Ri number is calculated as follows:

Ri = Gr Re 2 = gβ e (T hot -T ref )D 3 /ν 2 (vD/ν) 2 = gβ e (T hot -T ref )D v 2 (3.12)
where β e is the thermal expansion coefficient, Ri number presents the importance of natural convection compared to the forced convection. solid metal, air and mould) are considered as one effective medium throughout the domain, having weighted physical properties as functions of the phase fractions. Let us call φ the duct indicator (φ = 1 inside the tube, 0 in the mould), and α the local volume fraction of metal in the tube (α = 1 in the metal, 0 in the air), finally we also introduce the liquid metal fraction g l that is function of the temperature (g l = 1 in the liquid metal and 0 in the solid phase).

Consequently the different sub-domains are defined as follows:

φ = 1, α = 1, g l = 0 solid 1 liquid
The thermophysical properties of the effective medium are obtained by assuming a linear variation between the physical properties of each phase at the different interfaces, except for the effective conductivity calculated using Patankar [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] formulation due to strong conductivity contrasts between the phases. Therefore, the effective properties are weighted as follows:

ρ = α ρ metal + (1 -α)ρ g (3.14) µ = α µ metal + (1 -α)µ g (3.15) ρc p = φ α ρ metal (g l c p l + (1 -g l ) c p s ) + (1 -α)ρ g c p g + (1 -φ) ρ m c p m (3.16) k = φ α g l k l + (1 -g l )k s + (1 -α) k g -1 + (1 -φ) [k m ] (3.17)
Where the subscripts l, s, m, g and metal stand for liquid metal, solid metal, mould , gas and metal consecutively.

Note that here we assume that the density difference between the liquid and solid phases is The fluid flow is governed by the modified Navier-Stokes equation for an incompressible fluid.

Therefore, the continuity and the momentum equations are given by:

1. continuity equation ∂ρ ∂t + ∇.(ρ v) = 0 (3.18)
2. momentum equation

∂(ρ v) ∂t + ∇.(ρ v v) = -∇p + ∇.(µ∇v) + σ κ ∇α -Av -Bv (3.19)
where κ = -∇ ∇α |∇α| is the curvature of the metal-gas interface

The momentum conservation equation (3.19) contains additional source terms allowing for a general form of the equation all over the domain. The third term in the RHS accounts for the pressure jump at the liquid-gas interface arising from the surface tension between the two fluids.

The interface position (meniscus shape) is solved using the volume-of-fluid (VOF) method [START_REF] Hirt | Volume of fluid ( vof ) method for the dynamics of free boundaries *[END_REF],

now widely used in many CFD codes for solving interface problems on fixed grids (section 3.3.2). It consists in solving a transport equation of the metal volumetric fraction α:

∂α ∂t + ∇.(α v) + ∇. [v r α(1 -α)] = 0 (3.20)
where v r = v l -v g is the relative velocity of the two fluids modelled as the maximum velocity magnitude at the liquid-air interface. Since the surface curvature required for the determination of the surface tension force is a function of the volumetric fraction, this fraction distribution must be accurately determined, especially in the case of high physical properties ratios (density, thermal conductivity, ...), where small errors in the evaluation of α lead to significant errors in calculating the effective physical properties. Thus, eq. (3.20) contains an additional convective term known as the "compression" term active only within the interface region where α(1-α) = 0 ; its role is to get a sharper liquid-gas interface [START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF].

The last two terms in the RHS of eq. (3.19) are penalty terms. The first one depends on the solid/liquid interface position and imposes a zero velocity field in the solidified metal, it acts only in the metal region (φ = 1, α = 1) and is function of the liquid metal fraction g l . A

Carmen-Kozeny formulation is chosen:

A = φ α (1 -g l ) 2 g l 3 + b (3.21)
where b is a small constant to avoid singularity when dividing by zero for g l = 0. When the calculated temperature T is higher than the melting temperature, the liquid metal fraction is

38 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE
equal to 1 and the penalty coefficient vanishes, thus the momentum equation recovers its normal form. Alternatively, when the calculated temperature is lower than the melting temperature, solidification occurs and the liquid metal fraction is equal to zero. Thus, the penalty coefficient has a large value and results in a zero velocity field.

The second additional source term in the momentum equation ensures a zero velocity field in the mould region. Thus, the penalty coefficient is function of the mould indicator φ:

B = (1 -φ)P (3. 22
)
where P is a large constant. In the air/metal domains φ = 1 and the penalty coefficient turns to zero. However, in the mould region where φ = 0, the penalty coefficient has a large value ensuring a zero velocity solution.

The phase change model of the solidification of the pure metal uses a procedure adapted from the enthalpy-porosity method [START_REF] Brent | The Enthalpy-Porosity Technique For Modeling Convection-Diffusion Phase Change : Application To The Melting Of A Pure Metal[END_REF]. The advantage of this method is that it respects the one domain approach. The general energy conservation equation:

∂H ∂t + ∇.(H v) = ∇.(k∇T ) (3.23)
is transformed by deriving the relationship between the enthalpy and the temperature with respect to time and space in eq. (3.23), which leads to a formulation where temperature is the depending variable:

ρc p ∂T ∂t + ρc p ∇.(T v) = ∇.(k∇T ) -φ α ρ metal L f ∂g l ∂t (3.24)
One may notice that the phase change term involving the latent heat (second term in the RHS) applies only in the metal (φ = α = 1), and that another variable g l is still present in the energy equation. But, the latter is function of the temperature, for a pure metal a simple linear function relates the liquid volume fraction to the temperature [START_REF] Rösler | Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments[END_REF]:

g l =          0 T < T s T -T s T s ≤ T ≤ T l 1 T > T l (3.25)
where = T l -T s , the temperature difference between the metal liquidus and solidus. For a pure metal T l = T s but in order to avoid discontinuity is taken to have a small value so that the liquid fraction will approach a Heaviside step function. This temperature difference must be small enough to obtain a sharp interface between the solid and liquid phases but large enough so that the temperature belongs to the temperature interval T ∈ ]T s , T l [ where the latent heat is exchanged. Thus after running numerical tests the value of = 0.

∂g l ∂t =          0 T < T s 1 ∂T ∂t T s ≤ T ≤ T l 0 T > T l (3.26)
The temporal derivative of g l has a non zero value only in the transition zone from liquid to solid, where the latent heat is exchanged at the liquid-solid interface (2nd term in the RHS of eq. (3.24)). Introducing a Dirac delta function δ:

δ =        0 T < T s 1 T s ≤ T ≤ T l 0 T > T l (3.27)
The temporal derivative of g l can be written as:

∂g l ∂t = δ 1 ∂T ∂t (3.28)
Finally, eq.( 3.24) can now be written as follows:

ρc p ∂T ∂t + ρc p ∇.(T v) = ∇.(k∇T ) -φ α ρ metal δ L f 1 ∂T ∂t (3.29)

Numerical procedure

The set of equations resulting from the mathematical model presented above (section 3.

Spurious currents

The literature shows that implicit tracking methods present some limits to which numerical parasite currents appear and can destroy the interface [START_REF] Lafaurie | Modelling Merging and Fragmentation in Multiphase Flows with SURFER*[END_REF][START_REF] Scardovelli | Direct Numerical Simulation of Free-Surface and Interfacial Flow[END_REF][START_REF] Worner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF]. These currents are due to vortices that appear in the interface region despite the absence of any external forces. These currents are common in many numerical methods that deal with interfaces including the CFS model. There is not enough theoretical analysis on their origin, that is why they are also called "parasite" currents. To reduce the effect of these vortices, we should use a very small time step which increases the computational time [START_REF] Saldi | Marangoni driven free surface flows in liquid weld pools[END_REF], Galusinski et al. [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF] defined a time step criterion ∆t ≤ max 10µ∆x/σ, ρ∆x 3 /σ that should be respected in order to avoid the presence of these spurious currents. Even though it is difficult to give a real estimation of the magnitude of the spurious currents because they fluctuate in time, numerical experiments done by Bruno et al. [START_REF] Lafaurie | Modelling Merging and Fragmentation in Multiphase Flows with SURFER*[END_REF] showed that the maximum Reynold number is proportional to the Laplace number, and that the corresponding maximum velocity of these vortices is given by:

u max = Re La σ µ = λ σ µ (3.30)
where, λ = 10 -2 , another study reported a different value of this proportionality factor (λ = 10 -5 ) but no reference is mentioned (Tryggvason, unpublished lecture notes) [START_REF] Scardovelli | Direct Numerical Simulation of Free-Surface and Interfacial Flow[END_REF].

The spurious currents are more important when the capillarity forces control the flow, which means that the surface tension is large compared to the fluid viscosity. A dimensionless number that illustrates the ratio between the two forces is the capillary number presented earlier. Numerical simulations become difficult when this number is on the order of 10 6 or higher [START_REF] Scardovelli | Direct Numerical Simulation of Free-Surface and Interfacial Flow[END_REF],

unfortunately using the aluminium properties and the capillary tube diameter gives a Laplace number O(10 6 ). Therefore to avoid spurious currents the CFL criterium for infiltration problem is set to be less than 0.1 which gives a small timestep that varies between 10 -05 s and 10 -07 s and increases the calculation time.

Contact angle

The problem with the one domain formulation is that liquid-air front is not a physical interface and boundary conditions such as; pressure drop, contact angle etc ... cannot be specified explicitly as in 2 domain approaches. Instead the contact angle at the mould-liquid-air contact line is implicitly present in the calculation of the liquid-air curvature. The contact angle is a physical property that depends on the fluids and solid in question, and thus should be imposed to the system. Therefore some additional manipulation to the code should be done in order to include the right contact angle in the curvature calculation.

Let us recall the definition of the normal vector at the liquid-air interface used in the VOF method:

n = ∇α |∇α| (3.32)
and, the curvature of the free surface is determined by

κ = -∇. n (3.33)
As we can see using the definitions above do not permit us to apply the correct angle, thus instead of using eq.

Ω = π -θ (3.34)
The rotation of a vector about the y axis is obtained by applying an inner product between the vector and the corresponding rotation matrix

n =     cos(π -θ) 0 sin(π -θ) 0 1 0 -sin(π -θ) 0 cos(π -θ)     . S f (3.35) =     cos(π -θ) 0 sin(π -θ) 0 1 0 -sin(π -θ) 0 cos(π -θ)     .     0 0 1     =     sin θ 0 -cos θ    
This correction of the calculation of the normal vector n at the liquid-air interface permits the use of the real contact angle θ which is a specific property of the fluid and solid in contact. In our case the value of 150°is used.

Numerical Validation

The numerical validation of the code is done by comparing the 2D numerical solutions with analytical solutions for simple test cases available in the literature. In this section, the mould is not taken into consideration. The capillary tube is modelled as a cylindrical tube of 1.2 mm radius, and 120 mm length (unless mentioned differently) to be in accordance with the pore size of a typical metal foam mould and the tube diameter used in the fluidity tests. In order to ensure that the solution is mesh-independent, a mesh refinement study is done and for 15 cells/radius and ∆x = 0.01 mm the numerical solutions were in good agreement with the analytical ones with a relative error less than 1% in the range of parameters for this study. The simulation of the infiltration of liquid metal inside a capillary tube presents a difficulty which is the interface tracking done using the VOF method (presented earlier). A tube initially filled with air is filled by liquid aluminium that pushes the air out. The flow is driven by an applied pressure difference between the tube inlet and outlet of 100P a. The maximum velocity at the tube axis as function of the axial distance is plotted in Figure 3.8, we can identify two maximum velocity values that correspond to the maximum velocity in the liquid and in the air.

However the liquid-air interface maximum velocity is equal to the flow average velocity, that is why Figure 3.8 shows a minimum value that corresponds to the interface location. In order to understand more the two-phase flow dynamics, the radial velocity profiles in the two region: liquid and air, as well as at the interface are plotted in figs. for pressure driven flow the state is established at every cross-section, but it is not steady as function of time contrary to a flow driven by a constant velocity boundary condition. In fact, the characteristic diffusion time over a distance L is of order of [START_REF] Guyon | Hydrodynamique physique[END_REF]:

τ = L 2 ν (3.37)
where, ν is the cinematic viscosity and L the characteristic length which is the tube radius R for an internal flow in a tube. Therefore the diffusion time is around 3.5 s for liquid aluminium and 0.1 s for the air which explains why the flow is fully developed and steady in the air region and not in the liquid aluminium at τ air < t = 0.18 s << τ Al .

A special attention is given to the liquid-air interface velocity. In fact, after a certain time the interface shape remains the same, and displaces at the fluid average velocity; as stated by Dussan [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF] "the speed of the contact line Q/πa 2 ". In order to understand the motion of the two fluids and the interface dynamics we refer to [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF] where it is demonstrated that the displacing fluid undergoes a rolling motion and that the displaced fluid has a similar but more complex behaviour. To visualize this, the relative velocity vectors in the interface frame of reference are shown in Figure 3.12.

Internal forced convection

After validating the hydrodynamic behaviour of a one-phase flow, the thermal aspect must be verified. Therefore, internal forced convection is simulated by imposing a uniform velocity at the tube inlet. There are two possibilities of thermal conditions at the tube wall: fixed temperature or constant thermal flux, accordingly to which the thermal entrance length and the correlated dimensionless number of Nusselt are given as follows [START_REF] Bejan | Convection Heat Transfer[END_REF] :

• Fixed temperature at the tube wall Using the aluminium properties, the thermal entrance length for both cases (L th < 2 mm) is negligible in comparison with the mechanical one. The validation of the thermal behaviour is done by fixing a constant flux at the wall that corresponds to the energy extracted from the metal in order to lower its temperature from T in to T out :

Q1 = ṁ c p (T in -T out ) = ρ S a v c p (T in -T out ) (3.40)
where, S a = πR 2 the tube section.

This energy should be lost to the tube wall:

Q2 = S c q = S c k ∇T (3.41)
where, S c = 2πRL the exchanging tube surface. Then equalling the two quantities Q1 = Q2 , gives:

q = ρ R v c p (T in -T out ) 2L = k ∂T ∂r ( 3.42)
Now in order to calculate the numerical Nusselt number, we have to recall its definition:

N u L = h L k (3.43)
where, h the convective heat transfer coefficient. At the tube wall, the heat exchange between the fluid and the wall is given by:

q = h(T w -T f ) (3.44)
where the subscripts w and f refer to the wall and the fluid consecutively. Therefore, the convective heat transfer coefficient can be calculated by:

h = k ∂T ∂r T w -T f (3.45)
Where T w the temperature obtained numerically at the tube wall (r = R ), and T f the weighted average temperature calculated from the radial temperature profile T (r) in a cross-section beyond the thermal and mechanical entrance lengths (x = 100 mm)(see fig. 3.13).

T f = R 0 ρ 2π r v(r) T (r) dr ṁ = R 0 ρ 2π r v(r) T (r) dr ρ π R 2 v av (3.46)
For a Poiseuille flow, the radial velocity profile is given by :

v(r) = ∆pD 2 16µL 1 - r R 2 = V max 1 - r R 2 (3.
T f = 4 R 2 R 0 r 1 - r R 2 T (r) dr (3 .48)
The numerical result gives a Nusselt number of 4.3 to be compared to the usual correlation N u = 4.363, giving a relative error less than 1.5%.

Phase change problem

An analytical solution for inward solidification in a cylinder for a liquid initially at the fusion temperature may be found in the literature [START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF]. The numerical solution is compared to the one-dimensional pure conduction analytical solution: a negative uniform heat flux (Q = -1000 kW/m 2 ) is imposed at the outer boundary of a cylindrical tube filled with liquid aluminium at its melting temperature.The analytical solution gives the radial position of the solidliquid interface as function of time:

r(t) = R 0 2 + 2R 0 ρ L f Qt (3.49)
and the solidification time:

t sol = ρ L f R 0 2|Q| (3.50)
The solid-liquid interface position obtained numerically and analytically is shown in Figure

t sol = ρ( L f + c p ∆T )R 0 2|Q| (3.51)
and the time delay caused by the additional sensible heat that should be lost before solidification begins is:

t delay = ρc p ∆T R 0 2|Q| (3.52)
Therefore the liquid-solid interface position as function of time for t ≥ t delay is written as follows: 3.16 shows that in case of superheat the solid-liquid interface evolution as function of time has a similar allure as in solidification of a metal at its fusion temperature but with a time offset that is larger as function of the superheat. This result is expected since the metal has to loose its sensible heat in addition of the latent one before solidifying.

r(t) = R 0 2 + 2R 0 ρ L f Q(t -

Preliminary Simulations

The performance and accuracy of the numerical code has been tested by solving various subproblems: the Poiseuille flow problem, the dynamic problem of infiltration of a cylindrical capillary tube and the problem of inward solidification by purely diffusive heat transfer. The classical solutions are obtained with a good level of accuracy as shown earlier.

In the presentation of the simulation results, we will progressively tackle the problem of infiltration with solidification. In a first step, we will consider a one-phase flow of liquid metal in a cylindrical tube driven by a pressure difference, and cooled by the tube wall. We will then use the developed code to study the infiltration problem and finally we will present the results concerning the dynamics of solidification. The difficulty of a parametric study in this problem is that the operating conditions are not appropriate to the definition and use of the classical dimensionless numbers:

• the Reynolds number, Re: essentially based on a reference velocity, here difficult to define because the fluid is accelerated from rest under a driving pressure difference • the Stefan number, Ste: should be built on the temperature difference which drives the solidification process, between some imposed wall temperature and the melting temperature. Here the "cold" temperature is the internal surface mould temperature which evolves from an initial to a final value during the process.

As a consequence we could just provide indicative values of these parameters, but since they evolve with time, they cannot be used for the analysis of the results and the conditions of the foregoing simulations will be given in dimensional terms.

Solidification of a liquid metal flowing in a tube

The first problem to be simulated using the numerical code is the solidification dynamics of a fluid flowing in a duct without infiltration (no liquid-air interface) and directly cooled by its external wall (no external mould): the configuration consists of a cylindrical tube (radius 1.2 mm and 50 mm long) initially filled with liquid aluminium at a temperature of 975 K (initial superheat of 42 K). In the absence of the liquid/air interface, there is no Laplace force UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE exerted therefore the pressure difference needed to drive the flow is much lower than when dealing with infiltration problem. Thus a pressure difference of 50 Pa is imposed between the tube inlet and outlet, the flow is fully developed at t = 0s and has an average velocity of 0.18 m/s when a constant heat flux of (-1000 kW/m 2 ) is imposed at the external wall. There is no mould around the tube and solidification is provoked by extracting a constant heat flux.

The imposed heat flux is uniform and the effect of gravity is neglected, thus an axi-symmetrical solidification is expected.

The progression of the solid phase is described in Figure 3.17, where the liquid aluminium phase is in red and solid aluminium in orange. The tube is initially (t = 0s) completely filled with liquid aluminium. At early times ( t = 0.07s) aluminium solidification starts at the tube wall at a distance downstream of the entry, the solidified aluminium crust rapidly grows inward in the radial direction (b., c.). In this process the diameter of the liquid flow decreases and the fluid velocity decreases as well (Figure 3.18 -b) up to a point (here t = 0.645s) where the solidified aluminium blocks the flow at a distance from the inlet. Later on, solidification proceeds by pure diffusion towards the entry (f., g.) until complete solidification is achieved (t = 0.8s). In The evolution of the interface position as function of time in case of solidification in a flowing liquid has the same tendency as in pure diffusive solidification. For the same metal superheat of 42 K the main difference is the solidification time which is slightly lower in case of solidification of diffusive solidification which is expected since for a flowing liquid fresh "hot" liquid is constantly flowing inside the tube until solidification blocks the tube section which is translated by an additional heat source that should be extracted in order to solidify the liquid metal. This effect is more clear for higher velocities, as Figure 3.19 shows for the same thermal conditions when the fluid flows faster the solidification of the latter is slower. As presented earlier, in the

Infiltration of a capillary tube

The second problem of interest is the problem of the infiltration of a capillary tube, inserted in a cylindrical mould, initially filled with air by a liquid metal, here aluminium. In a first step, we will just consider the fluid flow problem, before the beginning of solidification at the wall so that the heat transfer is not taken into account and the tube section is constant.

Compared to most studies reported in literature where a constant velocity or parabolic velocity profile are imposed at the inlet of the duct, one specific feature of the present problem is that the fluid motion is driven by a pressure difference which must overcome the Laplace force at the meniscus and the inertia of the initially motionless fluid. As a consequence, the fluid is accelerated and its average velocity increases with time, and since the infiltration time is less than the diffusion time (∼ 3.5 s, section 3.5. 1.2), the flow is never permanent in the frame of our study.

Let us consider the infiltration of a cylindrical tube of radius R by a fluid of density ρ under a pressure difference ∆p. Let l(t) be the length of the infiltrate liquid in the tube. Since the fluid is initially at rest, the velocity is small and one can neglect the viscous forces compared to the acceleration of the fluid from rest to flowing and the fundamental law of dynamics yields: The relation between the infiltration length and the fluid velocity is plotted in Figure 3.

d

Description of the coupled infiltration and solidification problem

In a one-domain approach, the physical boundary between the tube and the mould is represented by the phase indicator φ that passes from value 0 in the mould to 1 in the tube. Therefore, the mesh in this transition zone should be extremely fine in order to obtain good accuracy.

As mentioned before the VOF method is chosen in order to determine the metal-air interface shape and follow its position as function of time (section 3. presents some limitations in the calculations of the effective thermo-physical properties due to the high contrast between the metal and air volumetric heat capacity and thermal conductivity.

A proposed solution will be presented in the following section. A Poiseuille one phase flow inside a cylindrical mould is simulated, where a pressure difference of 1400 P a is imposed between the mould inlet and outlet. Three mesh refinements were tested and are presented in Table 3.5. The flow is fully developed and permanent after ∼ 4 s, the velocity profiles obtained using the three mesh refinements are compared to the Poiseuille's solution in Figure 3. [START_REF] Cheung | Transient Freezing of Liquids in Tube Flow[END_REF] Figure 3.22: Velocity profiles for different mesh refinements Figure 3.22 shows that the solution is acceptable for mesh 3 (Table 3.5), where the relative error between the numerical and analytical solution is less than 1 %. Thus this mesh will be employed in the simulations of liquid metal infiltration or/and solidification inside a cylindrical mould.

As stated before, a special treatment is required for the calculation of the effective thermophysical properties as function of the metal/air volumetric fraction α. Lets consider the schem- 3. [START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF]. The effective conductivity is close to the air conductivity in this case which is correct, however the volumetric heat capacity affected to this calculation cell where α = 0.1 is overestimated due to the high contrast between the metal and air volumetric heat capacities. To overcome this limitation, a proposed numerical solution is to create a discontinuity of α value at the metal/air interface after the determination of α distribution and the reconstruction of the interface shape as follows: if α > 0.999 then α = 1 (3.56) else α = 0

Note that this new α distribution will be only used to calculate the effective thermo-physical properties. This manipulation prohibited the solidification to start at the metal/air interface, but our α limit choice of 0.999 was so extreme it affected the solution on the metal/mould interface where the metal volumetric fraction distribution was close but less than 1. Thus, the conductivity and volumetric heat capacity on some calculation cells near the metal/mould interface were those of the air, in other words the heat transfer was underestimated in those locations causing some oscillations on the temperature profile along the mould length. These miscalculated thermal flux are on one or two calculation cells near the mould/metal interface. However, we believe that the effect of this numerical error is not drastic and that the results are qualitatively correct. Nevertheless, this conclusion should be validated by changing the α limit and compar- We now consider that the tube is initially filled with air and inserted in a sand mould through which the liquid aluminium is injected and cooled by the mould. The annular cylindrical mould of length 200 mm (unless mentioned differently) has an inner radius of 1.2 mm and outer radius of 5 mm. The external wall of the mould is supposed to be adiabatic (zero thermal flux) and the mould and air are initially at ambient temperature (294 K) (unless mentioned differently).

The liquid aluminium flow is driven by a constant pressure difference of 1400 P a, sufficient to overcome the Laplace pressure at the metal-air interface. We consider that the liquid metal is initially superheated at T = 1045 K (∆T 0 = 112 K above the solidification temperature), which corresponds to an intermediate value of Stefan number of 0.33 built on this temperature difference. This means that the involved sensible heat is lower, but comparable, to the latent heat of the metal. The numerical simulation of the time evolution of the infiltration is shown in Figure 3.24 and the detail of the solidification process in Figure 3. [START_REF] Dairon | Développement d'une nouvelle technique d' élaboration de mousses d'acier par fonderie et caractérisation mécanique[END_REF].

At the beginning of the injection process, the liquid metal infiltrates the tube and expels the air initially present in the tube. The aluminium remains liquid as long as it loses its initial superheat, then solidification starts at the metal-mould contact in the vicinity of the metal-air interface (at t = 0.18 s), a zone which has exchanged heat with the outside on a longer time and then cools down earlier. As the metal flows in the mould, solidification proceeds in both axial and radial directions, until the flow is blocked up (at about t = 0.3 s).The flow is stopped even though the tube section is not closed due to the solidification at the metal/air interface. This behaviour is expected, since the metal close to the air interface is constantly in contact with the cooler mould, so that the temperature decrease in this region is more important. After immobilization of the liquid phase, we can observe diffusive phase change with axial solid growth until the whole In order to observe in more detail the solidification process, a zoom of the solidification zone is shown in Figure 3. [START_REF] Dairon | Développement d'une nouvelle technique d' élaboration de mousses d'acier par fonderie et caractérisation mécanique[END_REF]. One clearly sees that solidification is initiated at the liquid-mould contact a few mm from the air-liquid interface, which slows down the fluid flow and increases the heat exchange with the mould. At the triple contact point (metal/mould/air) the solidified metal blocks the flow of the liquid metal which causes the air to be trapped between the mould and the metal (t = 0.19 s). Then the liquid metal bifurcates around the solidified metal and continue to flow for some additional millimetres leaving a thin layer (2 calculation cells) of air between the mould and metal. The solidified zone progressively expands along the tube and in the radial direction until complete blockage due to the growth of the solid crust in the radial direction and the solidification at the metal/air interface (t = 0.225 s). The liquid flow is then stopped and solidification proceeds in the upstream direction by pure heat diffusion.

It is of interest to observe the isothermal lines map during the solidification process of a cylindrical mould by liquid aluminium. At first in the solid crust formed at the metal/mould interface the isothermal lines t = 0.3 s for T ∈ [932, 933[ K are shown in Figure 3. [START_REF] Dantzig | Solidification[END_REF]. Three main observations can be draw from this map, first that the solid-liquid interface presents some irregularities, second that the thermal flux on the metal/mould vicinity is nearly constant at this time t = 0.3 s as the isothermal lines are parallel and have the same radial slope, at last the thermal flux on the metal/air interface can also be considered as constant if the irregularities were ignored. Also one can clearly sees the thin air layer trapped between the mould and air few millimetres ahead of the metal/air interface, as mentioned earlier. the mould and air in the vicinity of the metal/air interface is presented in Figure 3. [START_REF] Epstein | Complex Freezing-Melting Interfaces in Fluid Flow[END_REF]. The isothermal lines in the air domain are regular and parallel, the heat flux at the tube axis is null which is expected due to the cylindrical symmetry of the problem. The isothermal lines in the mould region near the metal/mould interface present some irregularities similar to the one seen in the solidified metal in Figure 3. [START_REF] Dantzig | Solidification[END_REF]. Another observation is that in order to conserve the thermal flux, the isothermal lines radial slope in the air is steeper than in the mould which is expected since the mould has a greater conductivity than the air. Finally both the mould and the air are still cold as their temperature is about 300 K just few millimetres far from the metal. We chose to have a closer look at the temperature variation in the solidified metal at t = 0.3 s (after the blockage of the liquid flow) and at an axial section x = 75 mm. The temperature radial profile is presented in Figure 3. [START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF]. As expected the temperature is almost constant in the metal region where the thermal conductivity is high (∼ O(10 2 W/(mK)). However, a slope rupture is observed at the mould/metal interface (x = 1.19623 mm), this is due to the numerical error mentioned earlier in section 3. 6.3.1 where the conductivity at the mould/metal interface is falsely taken to be the air conductivity. By changing the scale, the variation of the temperature in the metal region is presented in Figure 3. [START_REF] Fukusako | Solidification of Pure Liquids and Liquid Mixtures Inside Ducts and Over External Bodies[END_REF]. The temperature is constant in the metal liquid phase then at the solid/liquid interface a slope rupture is observed and the temperature decreases quasi-linearly in the solidified metal.

The local heat flux at this axial section (x = 75 mm) is presented in Figure 3. [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]. The heat flux increases radially to reach a peak at the mould/metal interface then decreases in the mould with the radial position. This behaviour is expected as the heat transfer is triggered by the high temperature difference between the metal and mould. Figure 3.31 also shows the underestimated thermal flux value just before the thermal peak near the mould/metal interface due to the numerical error mentioned earlier (section 3. [START_REF] Bejan | Convection Heat Transfer[END_REF]

Influence of the operating conditions

The main focus of our study is to analyse the solidification process during the infiltration of a cylindrical mould by liquid aluminium and to document the influence of the key operational parameters, namely, the imposed pressure difference, the initial superheat of the injected liquid metal and the initial temperature of the mould on the infiltration length and the infiltration time.

The infiltration length is the length in millimetres of the solidified metal inside the mould, while the infiltration time is the flow blockage time after which the solidification proceeds by pure diffusion.

Influence of the pressure difference

We have seen that, under usual operating conditions, the first important parameter is the pressure difference imposed to the system. This difference must be large enough to overcome the Laplace pressure at the meniscus (in our conditions for liquid Al and an inner tube radius of 1. However, Figure 3.33 shows that the infiltration time decreases in a parabolic tendency as function of the pressure difference. This relation was also found in literature [START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF] and could be explained as follows: when increasing the pressure difference the flow velocity increases, and the liquid metal infiltrates further in the mould, this is translated by a greater contact surface between the metal and cold mould which causes the metal to freeze more rapidly. This is a first indication showing that the length of the mould does not need to exceed a given value for a given set of conditions. Of course the process, and particularly the maximum infiltration length, is also driven by the solidification dynamics, which also depends on the thermal conditions. This is the focus of the next sections.

Influence of the liquid superheat

The latter configuration is intended to schematically simulate the infiltration of the pores of a porous mould by a solidifying liquid metal. The thickness of the mould is also an important parameter and here it has been chosen so that the mould may absorb the latent heat contained in the initial liquid without an excessive temperature rise. Besides, the system is supposed to have negligible heat losses towards the ambient. These boundary conditions are not fully optimal, and it would be interesting to consider a heat exchange coefficient towards the environment, but Thus, the flow average velocity in the fluidity tests is higher than the velocities reached due to an imposed pressure difference between the tube inlet and outlet in the numerical simulations. This can explain why there is an important difference between the infiltrated lengths obtained numerically and experimentally (Figure 3.34 & Table 3.4).

In conclusion, in the scoop of this thesis a proper comparison between the numerical results and the experiments measurements was not performed, but we can clearly conclude that the tendency is the same. The validation of the model against the fluidity tests is one of our prior perspective for future work.

Influence of the initial temperature of the mould

Another key parameter of the metal solidification dynamics is the initial thermal state of the mould, since it determines its ability to cool and solidify the liquid metal. We thus consider different values of the mould temperature (294 K, 450 K, 600 K and 800 K) for an initial superheat of 52 K (T = 985 K). The results are presented in Figure 3. [START_REF] Harris | Phase-change phenomena in porous media -a non-local thermal equilibrium model[END_REF] length is very sensitive to the initial mould temperature, and again a linear relation between the infiltration length and time and the mould preheat is detected. This behaviour agrees with the literature (section 2. 1.1.2). The plots show that when the mould is preheated to higher temperature, the liquid metal infiltrates greatly the tube, because at high temperatures the mould cannot easily evacuate the heat from the liquid. A more complete analysis should be done by changing the boundary conditions on the external wall of the mould in order to account for possible heat losses.

Conclusions

In the objective to develop a macroscopic model of the infiltration and solidification of a molten metal inside a porous medium, a local model was first developed where one strut of the metallic foam is considered as a cylindrical tube. Thus the infiltration and solidification of a liquid metal inside a cylindrical mould is studied and presented in this chapter.

Fluidity tests of pure aluminium (Al 99, 5%) for three initial melt superheat temperature (∆T = 112, 172, 222 K) were performed by the CTIF [3]. After treating the surface of the metallic specimen, the solidification modes of the aluminium frozen in the spiral mould is identify by comparing the observations to the literature. Exquiaxed grains were identify around the porosities and at the specimen tip, this is due to the impurities found in Al 99, 5%.

A one domain approach is chosen for modelling the infiltration and solidification of liquid metal inside a mould. Three phase indicators or fractions are introduced: φ the mould phase indicator, α the metal/air volumetric fraction and g l the liquid metal fraction. The effective properties of the equivalent fluid are calculated as weighted averages using these fractions. The metal/air interface shape and position is determined using the VOF method. The solidification of the liquid metal is treated using the enthalpy-porosity method. The conditions of zero velocity field in the mould and solidified metal are ensured using two penalty terms in the momentum conservation equation. The mathematical model is coded in OpenFoam CFD tool.

The model is validated against existing analytical solutions (Poiseuille's flow, internal forced convection, inward cylindrical solidification by diffusion) and good accuracy was found. However when treating the coupled problem of infiltrating and solidification of liquid metal inside a cylindrical mould, a numerical error is detected. This error is due to the miscalculation of the conductivity at some locations (one to two calculation cells) on the metal/mould interface. The effect of this error is believed to not affect the global solution but this information should be validated by further simulations. parametric study gave the relationship between the performance of the process and the relevant operating parameters such as the initial superheat, the imposed pressure difference between the tube inlet and outlet and the initial mould temperature. As recalled in the introduction, the analysis of this coupled heat transfer and fluid flow problem is the first step of the modelling approach, before developing the macroscopic model of the solidification inside a porous matrix under similar conditions, that will be presented in the next chapter.

Macroscopic model of solidification in porous medium

Introduction

The industrial context of the present analysis concerns the elaboration of aluminium metal foams which are characterized by very large porosity values (greater than 80%). These porous structures are obtained after solidification of an infiltrated molten metal in a porous casting mould.

To our knowledge, very few theoretical studies have been devoted to this problem and averaged representation for solidification in porous media still remains an open problem.

We reviewed in Chapter 2 the work on infiltration and phase change in a porous medium. To our knowledge there is no studies of the coupled problem of infiltration and solidification in a porous medium that involves four phases. Most of the existing studies on solidification in a porous medium consider local thermal equilibrium (LTE) between all three phases: liquid/solid metal and the mould. This assumption is false when dealing with the phase change of metals inside a sand mould, because metals have a high fusion temperature that exceeds the burning temperature of the mould. Thus considering LTE between the metallic mixture and the mould returns to consider that the mould is at the metal fusion temperature. Recent studies that consider local thermal non equilibrium (LTNE) between the metallic mixture and the mould couple the two energy conservation equations by calculating a heat exchange coefficient based on existing correlation of the latter [START_REF] Liu | Numerical modeling for solid-liquid phase change phenomena in porous media: Shell-and-tube type latent heat thermal energy storage[END_REF][START_REF] Alomar | Simulation of complete liquidvapour phase change process inside porous evaporator using local thermal non-equilibrium model[END_REF] or the Nusselt number [START_REF] Chen | Experimental and numerical study on melting of phase change materials in metal foams at pore scale[END_REF].

Therefore our objective is to derive a macroscopic model for solidification of a pure molten phase in an homogeneous porous media. The full treatment of this problem would imply to take into account the infiltration by the molten metal, the displacement of the liquid-air interface and therefore to consider four phases. In this first analysis, the complexity will be reduced by assuming that the casting mould is fully infiltrated by the liquid metal (no solidification during infiltration). Under these circumstances the solidification is governed by diffusion involving

Local problem

As previously said, the infiltration of the homogeneous mould by a molten liquid phase is assumed to be completed. In other words, the porous mould is saturated by a motionless liquid phase whose volume fraction is equal to the porosity of the mould. Solidification process is therefore governed by a diffusion process involving three phases: the liquid metal (β-phase), the solidified metal (γ-phase) and the mould (σ-phase) where the local energy conservation equations are given by

∂ ∂t (ρ β H β ) = ∇ • (k β ∇T β ) in the β phase (4.1)
∂ ∂t (ρ γ H γ ) = ∇ • (k γ ∇T γ ) in the γ phase (4.2) ∂ ∂t ((ρ σ c pσ ) T σ ) = ∇ • (k σ ∇T σ ) in the σ phase (4.3)
where H β = c p β T β + L f is the enthalpy of the liquid metal, L f represents the heat latent and

H γ = c pγ T γ is the solid enthalpy. k i (for i = β, γ, σ) is the thermal conductivity of phase i.
The associated boundary conditions are written at the three interfaces locations A βγ , A βσ , and A γσ :

n βγ • (k γ ∇T γ -k β ∇T β ) = n βγ • w βγ (ρ γ H γ -ρ β H β ) at A βγ (4.4) T β = T γ = T f usion at A βγ (4.5) T γ = T σ at A γσ (4.6) 76 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE n γσ • k γ ∇T γ = n γσ • k σ ∇T σ at A γσ (4.7) T β = T σ at A βσ (4.8
)

n βσ • k β ∇T β = n βσ • k σ ∇T σ at A βσ (4.9)
In the next section this local problem will be up-scaled using the volume averaging method [START_REF] Whitaker | The Method of Volume Averaging[END_REF].

Volume averaging

T β = 1 V V β T β (x + y) dV y (4.10)
where V β is the volume of the β-phase contained within the averaging volume and where dV y is used to indicate that the integration is done with respect to the variable y.

T β β = 1 V β V β T β (x + y) dV y (4.11)
The two above-averaged values are related by the following relation:

T β = ε β T β β (4.12)
where ε β is the β-phase volume fraction, known as porosity when the porous medium is saturated by one phase β.

Finally, averaged conservation equations can be obtained by using spatial and temporal partial derivative theorems given by Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF] :

∇T β = ∇ T β + 1 V A βσ n βσ T β dA (4.13) ∂T β ∂t = ∂ T β ∂t - 1 V A βσ n βσ .w βσ T β dA (4.14)
where w βσ is the interface velocity and n βσ is the unit normal vector at the solid-liquid interface After the brief presentation of the volume averaging procedure and theorems for a 2-phases system (section 4. 3.1), let us now apply this method to the 3-phases system presented in Figure 4. 2. We consider the volume V presented in Figure 4.2 chosen to be large enough to contain the three phases β, γ and σ, and yet small enough to respect scale separation l β << L, where l β is the local length scale for phase β and corresponds to the pore diameter of the mould (distance between two cells), and L is the macroscopic scale which corresponds to the mould size. We first consider the averaging of the mass conservation equation for the motionless β-phase (v β = 0)

(Figure 4.1) [86].

Mass conservation

∂ρ β ∂t = ∂ ∂t ε β ρ β β - 1 V A βγ ρ β n βγ . w βγ dA βγ - (4.15) 1 V A βσ ρ β n βσ . w βσ dA βσ w βσ =0 = 0
Since the mould is rigid, the interfacial phase velocity w βσ = 0 and assuming that the intrinsic average of the density ρ β β = ρ β , the above equation reduces to

ρ β ∂ε β ∂t = 1 V A βγ ρ β n βγ . w βγ dA βγ Phase change rate = ṁβ (4.16)
where ṁβ is the phase change rate. Similar average for the γ-phase gives,

ρ γ ∂ε γ ∂t = - 1 V A βγ ρ γ n βγ . w βγ dA βγ = ṁγ (4.17)
Adding Eqs. (4. [START_REF] Calmidi | Forced Convection in High Porosity Metal Foams[END_REF]) and (4.17) leads to

∂ ∂t (ε β ρ β + ε γ ρ γ ) = 1 V A βγ n βγ • w βγ (ρ β -ρ γ )dA βγ shrinkage (4.18)
where the integral on the right hand side represents the shrinkage due to phase change. This shrinkage is generally small since the fraction δ = 1 -ρ β /ρ γ ≤ 10 % for most metals [START_REF] Hawkins | 8 -Thermochemical data[END_REF].

We will assume it is negligible and therefore equation (4.18) becomes

∂ ∂t (ε β ρ β + ε γ ρ γ ) = 0 (4.19)
This result was obviously expected since it confirms that mass conservation is satisfied at the macroscopic scale.

Energy conservation equations

The second step in determining the macroscopic non closed form is to apply the volume averaging method to the heat diffusion equation. Let us develop the macroscopic energy equation for the liquid metal (β-phase). Since the interfacial velocity of the mould/liquid interface is null w βσ = 0, the averaged enthalpy accumulation takes the form:

∂ ∂t (ρ β H β ) = ∂ ∂t ρ β ε β H β β - 1 V A βγ ρ β H β n βγ . w βγ dA βγ (4.20)
while the averaged diffusion term is given by

∇ • (k β ∇T β ) = ∇ • k β ∇T β + 1 V A βγ n βγ • k β ∇T β dA βγ (4.21) + 1 V A βσ n βσ • k β ∇T β dA βσ
In addition, we have

k β ∇T β = k β ∇ T β + 1 V A βγ n βγ T β dA βγ + 1 V A βσ n βσ T β dA βσ (4.22) with ∇ T β = ∇ ε β T β β = ∇ε β T β β + ε β ∇ T β β (4.23)
Let us consider γ β the phase indicator of the β-phase. By definition γ β = ε β , and applying the average operator on the spacial derivative of γ β gives:

∇γ β = ∇ε β + 1 V A βγ n βγ dA βγ + 1 V A βσ n βσ dA βσ = 0 (4.24)
therefore,

∇ε β = - 1 V A βγ n βγ dA βγ + 1 V A βσ n βσ dA βσ (4.
k β ∇T β = k β ε β ∇ T β β + 1 V A βγ n βγ T β dA βγ + 1 V A βσ n βσ T β dA βσ (4.26)
where it has be considered that T β β (r) ∼ T β β (x). Under these circumstances, the nonclosed energy conservation equation for the β-phase takes the form

80 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE ∂ ∂t ρ β ε β H β β - 1 V A βγ ρ β H β n βγ . w βγ dA βγ = ∇ • k β ε β ∇ T β β + (4.27) ∇ • k β V A βγ n βγ T β dA βγ + k β V A βσ n βσ T β dA βσ + 1 V A βγ n βγ • k β ∇T β dA βγ + 1 V A βσ n βσ • k β ∇T β dA βσ
Similarly we obtain the non-closed energy conservation equation for the γ-phase

∂ ∂t (ρ γ ε γ H γ γ ) - 1 V A γβ ρ γ H γ n γβ . w γβ dA γβ = ∇ • (k γ ε γ ∇ T γ γ ) + (4.28) ∇ • k γ V A βγ n γβ T γ dA γβ + k γ V Aγσ n γσ T γ dA γσ + 1 V A γβ n γβ • k γ ∇T γ dA βγ + 1 V Aγσ n γσ • k γ ∇T γ dA βσ
Since the thermo-physical properties of the liquid and solid metal remain close during the solidification process, local thermal equilibrium between β and γ phases is considered in the averaging volume. Therefore, the intrinsic average temperatures for the metal phases are considered to be equal, giving rise to a one-equation model for both βand γ -phases

T β β ∼ T γ γ = T m (4 .29)
This approximation means that the liquid-solid phase change at the macroscopic scale takes place when the average temperature T m = T fusion . Under these circumstances, adding equations (4. [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF]) and (4.28) and taking into account the boundary condition eq.(4.4) leads to

∂ H m ∂t = ∇ • ((ε β k β + ε γ k γ ) ∇ T m ) (4.30)
+ ∇ •   k=γ,σ k β V A βk n βk T β dA βk + k=β,σ k γ V A γk n γk T γ dA γk   + k=β,γ 1 V A kσ n kσ • k k ∇T k dA kσ where H m = ε β ρ β H β β + ε γ ρ γ H γ γ (4.
ε σ ρ σ c p σ ∂ ∂t T σ σ = ∇ • (ε σ k σ ∇ T σ σ ) (4.32) + ∇ •   k=β,γ k σ V A σk n σk T σ dA σk   + k=β,γ 1 V A σk n σk • k k ∇T k dA σk
Equations (4.30) and (4.32) represent the non-closed macroscopic energy conservation of the diffusive phase change of a liquid inside a porous medium. In order to close the problem, the deviations T β , T γ and T σ must be expressed as a function of the averaged temperature fields.

This is done in the next sections through the derivation of the associated closure problems.

Deviation problems

Writing the deviation problems consists in deriving the governing differential equations for the variable deviations, in our case T i for i = β , γ and σ. The conservation equations are obtained by subtracting the non-closed macroscopic energy equations (4.27), (

) after multiplication by the corresponding phase volume fractions and introduction of the spatial decomposition. This derivation is detailed below only for the β-phase and the associated expressions for γand σ-phases will be deduced by analogy.

The introduction of the spatial decomposition into the local equation for the β-phase (4.1) gives

ε β ∂ ∂t ρ β H β β + H β = ε β ∇. k β ∇ T β β + T β (4.33)
Subtracting equation (4.27) from (4.33) gives

ε β ρ β ∂ H β ∂t -ρ β H β β ∂ε β ∂t + 1 V A βγ ρ β H β n βγ . w βγ dA βγ = (4.34) ε β ∇.(k β ∇ T β ) -k β ∇ T β β ∇ε β -∇ •   k=γ,σ k β V A βk n βk T β dA βk   - k=γ,σ 1 V A βk n βk . k β ∇ T β dA βk - k=γ,σ 1 V A βk n βk . k β T β β dA βk
Taking into account Eq. (4. [START_REF] Calmidi | Forced Convection in High Porosity Metal Foams[END_REF]), the LHS of equation (4.34) reduces to

ε β ρ β ∂ H β ∂t + 1 V A βγ ρ β H β n βγ . w βγ dA βγ = RHS (4.
1 V A βk n βk . k β T β β dA βk = -k β ∇ T β β ∇ε β (4.36)
and finally, the deviation equation for the β-phase takes the form

ε β ρ β ∂ H β ∂t + 1 V A βγ ρ β H β n βγ . w βγ dA βγ = ε β ∇ • k β ∇ T β - (4.37) ∇ •   k=γ,σ k β V A βk n βk T β dA βk   Non-local diffusion - k=σ,γ 1 V A βk n βk • k β ∇ T β dA βk
where the second term of the RHS represents a non-local diffusion contribution. The estimate of the two integral terms in the RHS of equation (4.37) are given by

∇ •   k=γ,σ k β V A βk n βk T β dA βk   ∼ O k β T β l β L (4.38) k=σ,γ 1 V A βk n βk • k β ∇ T β dA βk ∼ O k β T β l 2 β ( 4.39)
and due to scale separation l β L (Figure 4.2), we have

O k β T β l β L O k β T β l 2 β (4.40)
In addition, according to Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF], quasi-stationarity also results from scale separation and therefore equation (4.37) reduces to

1 V A βγ ρ β H β n βγ . w βγ dA βγ = ε β ∇ • k β ∇ T β - (4.41) k=σ,γ 1 V A βk n βk • k β ∇ T β dA βk
Let us note that the remaining terms in the RHS have the same order of magnitude.

ε β ∇ • k β ∇ T β ∼ k=σ,γ 1 V A βk n βk • k β ∇ T β dA βk ∼ O k β T β l 2 β (4.42)
while the LHS term can be estimated as

1 V A βγ ρ β H β n βγ . w βγ dA βγ ∼ O w βγ ρ β H β l β (4.
O   w βγ ρ β c p β T β l β   ∼ O k β T β l 2 β (4.45)
This manipulation allows us to estimate the interfacial velocity w βγ as

w βγ ∼ O k β (ρc β ) β l β = α β l β ∼ 10 -2 ms -1 (4.46)
Note that in practical situations, w βγ is generally found to be at least one order of magnitude smaller than the above estimation. Under these circumstances, the interfacial velocity can be neglected, thus eq.(4.41) now reduces to

0 = ε β ∇ • k β ∇ T β - k=σ,γ 1 V A βk n βk • k β ∇ T β dA βk (4.47)
Identical form can be written for the γ and σ-phases. Therefore, the general form for i = β, γ, σ, can be written under the form

0 = ∇ • k i ∇ T i - k =i 1 V i A ik n ik • k i ∇ T i dA ik (4 .48)
The spatial decomposition is now introduced into the boundary conditions. Due to thermal equilibrium we have at the β -γ interface

T β = T γ at A βγ (4.49)
and

n βγ • (k γ -k β ) ∇ T m + n βγ • k γ ∇ T γ -k β ∇ T β = (4.50) -n βγ • w βγ ((ρ β c pβ ) -(ρ γ c pγ )) T m + T β + ρ β L f at A βγ Since ρ β c p β ∼ ρ γ c p γ
, the above equation reduces to

n βγ • (k γ -k β ) ∇ T m + n βγ • k γ ∇ T γ -k β ∇ T β = (4.51) -n βγ • w βγ ρ β L f at A βγ
At this stage, one of the most important difficulties in modelling processes submitted to the evolution of the micro-structure due to phase transformation lies in the treatment of the interfacial velocity w βγ in the deviation problem. In other words, the challenge consists in accounting for the interfacial velocity in the associated closure problems and therefore in the determination 84 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE of the effective transport properties. This would imply a "dynamic" local representation of the micro-structure evolution coupled to the transient evolution of the macroscopic fields (averaged temperature in the present case). One of the answers would be to develop a downscaling procedure in order to go back to the determination of the local fields which are at the origin of the micro-structure evolution [START_REF] Angelini | Premixing-related behavior of steam explosions[END_REF]. This is not the objective of the present analysis where a simplified boundary condition is obtained by making an estimate of the different terms of equation (4.51).

Indeed, due to the estimate of the interfacial velocity w βγ (eq. (4.46)) we deduce that

n βγ • w βγ ρ β L f = O 10 6 (4.52)
while for the two terms of the left hand side of equation (4.51) we have

n βγ • (k γ -k β ) ∇ T m = O (k β -k γ ) T m L (4.53
)

n βγ • k γ ∇ T γ -k β ∇ T β = O k γ T γ l γ - k β T β l β (4.54)
For foam metal moulds, the pore size is about 2 to 3 mm, thus l β ∼ l γ ∼ 10 -3 m while the mould size can vary from a few centimetres to one meter, thus L ∼ 1, 10 -1 m. Therefore, using the thermo-physical properties of aluminium leads to estimate that the LHS is about two orders of magnitude greater than the RHS. Hence, boundary condition eq. (4.51) reduces to

n βγ • (k γ -k β ) ∇ T m source term + n βγ • k γ ∇ Tγ -k β ∇ Tβ = 0 (4.55)
After this simplification the term containing the velocity of the metal solid-liquid interface as well as the latent heat does not figure in the boundary conditions, therefore only one source term remains: ∇ T m . Now at the mould-metal (liquid and solid) interfaces A kσ for k = β, γ, we have

T σ = T k -( T σ σ -T m ) source term at A kσ (4.56) n kσ . k k ∇ T k -k σ ∇ T σ = n kσ .   kσ ∇ T σ σ source term -k k ∇ T m source term    at A kσ (4.57)
In the above boundary value problem two source terms are identified ( T σ σ -∇ T m ) and ∇ T σ σ . Therefore, the temperature deviation for each phase (i = β, γ, σ) can be decomposed under the form

T i = b i .∇ T m + c i .∇ T σ σ + α i ( T m -T σ σ ) (4.

Closure problems

The closure problems for the determination of variables b i , c i and α i are obtained by replacing expression (4.58) in the deviation equation (4.48). Since the macroscopic source terms are constant in the averaging volume, their derivatives are zero and therefore the method of superposition allows to identify, for each phase i = β, γ, σ, three boundary value problems associated to the closures variables b i , c i and α i . Hence, the associated three closure problems phases i = β, γ, σ are:

• Problem I 0 = ∇. (k i ∇b i ) - k =i 1 V i A ik n ik . k i ∇b i dA ik (4.59) b β = b γ at A βγ (4.60
)

n βγ • (k γ -k β ) I + n βγ • (k γ ∇b γ -k β ∇b β ) . I = 0 at A βγ (4.61) b σ = b n at A nσ , n = β, γ (4.62)
n nσ . (k n ∇b n -k σ ∇b σ ) . I = -n nσ . k n I at A nσ , n = β, γ (4.63) b i (r + l j ) = b i (r), j = 1, 2, 3 (4.64)
• Problem II

0 = ∇. (k i ∇c i ) - k =i 1 V i A ik n ik . k i ∇c i dA ik (4.65) c β = c γ at A βγ (4.66
)

n βγ • (k γ ∇c γ -k β ∇c β ) = 0 at A βγ (4.
c i (r + l j ) = c i (r), j = 1, 2, 3 (4.70)
• Problem III [START_REF] Sabatino | A review on the fluidity of Al based alloys[END_REF])

0 = ∇. (k i ∇α i ) - k =i 1 V i A ik n ik . k i ∇α i dA ik (4.
α β = α γ at A βγ (4.72
) [START_REF] Sarler | Stefan ' s work on solid-liquid phase changes[END_REF])

n βγ • (k γ ∇α γ -k β ∇α β ) = 0 at A βγ (4.
α σ = α n + 1 at A nσ , n = β, γ (4.74) n nσ . (k n ∇α n -k σ ∇α σ ) = 0 at A nσ , n = β, γ (4.75
)

α i (r + l j ) = α i (r), j = 1, 2, 3 (4.76)
The resolution of the above three classical Laplacian problems in a local unit cell leads to the determination of the effective properties at a macro scale. In the next section we will present the closed form of the macroscopic energy equations then we will proceed to the numerical resolution.

Averaged equations

In order to obtain the closed form of the average conservation equations, the deviations (4.58) are introduced into the non closed equations (4.30) and (4.32) where intrinsic averaged enthalpies have been defined as: Therefore, the closed form of the averaged equations of the solidification of a metal inside a porous medium takes the form :

H β β = c p β T m + L f (4.
∂ ∂t ((ρc p ) m T m ) + ρ β L f ∂ε β ∂t -v mm .∇ T m -v mσ .∇ T σ σ (4.79) = ∇ • [K mm .∇ T m + K mσ .∇ T σ σ ] -h m ( T m -T σ σ )
where

(ρc p ) m = ρ β ε β c p β + ρ γ ε γ c p γ ∂ ∂t ρ σ ε σ c p σ T σ σ -v σm .∇ T m -v σσ .∇ T σ σ (4.80) = ∇ • [K σσ .∇ T σ σ + K σm .∇ T m ] -h σ ( T σ σ -T m )
The first terms in the LHS of eqs. convective effective transport coefficients that will be presented shortly.

In order to solve the above equations we first have to determine the effective properties. This is done by solving the associated closure problems (section 4. 7.2.1). Taking into consideration the boundary conditions, that will be presented in the next section (section 4. [START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF].2), give the following simplified expressions of the effective properties.

• Metallic mixture effective properties • Mould effective properties

K mm k β = (ε β + κ 1 ε γ ) I + (1 -κ 1 ) V A βγ n βγ b β dA βγ (4.81) + κ 1 V Aγσ n γσ b γ dA γσ + 1 V A βσ n βσ b β dA βσ K mσ k β = (1 -κ 1 ) V A βγ n βγ c β dA βγ + κ 1 V Aγσ n γσ c γ dA γσ + 1 V A βσ n βσ c β dA βσ (4.82) v mm = 1 V A βσ n βσ k β .∇b β dA βσ + 1 V Aγσ n γσ k γ .
K σσ k σ = ε σ I + 1 V A σβ n σβ c σ dA σβ + 1 V Aσγ n σγ c σ dA σγ (4.86) K σm k σ = 1 V A σβ n σβ b σ dA σβ + 1 V Aσγ n σγ b σ dA σγ (4.87) v σσ = 1 V A σβ n σβ . k σ ∇c σ dA σβ + 1 V Aσγ n σγ . k σ ∇c σ dA σγ (4.88) v σm = 1 V A σβ n σβ k σ .∇b σ dA σβ + 1 V Aσγ n σγ k σ .∇b σ dA σγ (4.89) h σ = 1 V A σβ n σβ . k σ ∇α σ dA σβ + 1 V Aσγ n σγ . k σ ∇α σ dA σγ (4.90)
where κ 1 , κ 2 are the conductivity ratios:

κ 1 = k γ k β (4.91) , κ 2 = k σ k γ (4.92)
and K ii is the effective thermal conductivity associated to diffusion in the metallic and mould phases, while K ij represents the coupled conductivity associated to the diffusion known as cross coefficient conductivity. v ij represents the pseudo velocity convecting the average temperatures.

Generally, it has been found that the contribution of the pseudo convective terms is very small compared to the diffusive one (section 4. The multiphysics software Comsol 4.4 is chosen for the numerical resolution of the closure problems. But first, in order to validate the numerical approximation of the different closure problems, we begin by solving the classical Maxwell's closure problem [START_REF] Whitaker | The Method of Volume Averaging[END_REF] by considering a two phase system β (fluid) and σ (solid), and dividing all the terms in the system of equations by the fluid conductivity k β . This leads to:

0 = ∇ 2 b β in β -phase (4.93) 0 = ∇ 2 b σ in σ -phase (4.94) b β = b σ at A βσ (4.95
)

n βσ • (κ -1) + n βσ • (κ∇b σ -∇b β ) = 0 at A βσ (4.96) b β = 0 at A βe (4.97)
where κ = k σ k β represents the conductivity ratio.

For the two phase system presented above, an analytical expression of the effective conductivity for a 2D arrays of non touching cylinders is given by Ochoa-Tapia et al. [START_REF] Ochoa-Tapia | Diffusive transport in two-phase media: spatially periodic models and maxwell's theory for isotropic and anisotropic systems[END_REF] :

K ef f k β = 2κ -ε β (κ -1) 2 + ε β (κ -1) (4.98)
In order to compare the numerical results to the above expression, Chang's unit cell [START_REF] Chang | Multi-scale analysis of effective transport in periodic heterogenous media[END_REF] is chosen (Figure 4.3). Two parametric studies were performed, where for a fixed value of κ = 10 the fluid volume fraction ε β varies from 0.1 to 0.9. Then, for ε β = 0.64 the conductivity ratio κ takes values between 1 and 10 4 . The effective dimensionless conductivity of the 2-phase system is given by:

K ef f k b = (ε β + ε σ κ) I + (1 -κ) V 1 V A βσ n βσ b β dA βσ (4.99)
For the symmetrical Chang's unit cell, the effective conductivity tensor is isotropic and a single component has to be evaluated, therefore the x projection of this tensor is compared to the The effective conductivities obtained for two other phase system geometries were also confronted to the analytical solution Eq.(4.98) since according to Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF] the geometry hardly influences the results.

Square array of squares

Another simple geometry is a square array of squares where the unit cell is presented in Figure 4. [START_REF] Bejan | Convection Heat Transfer[END_REF]. This geometry was also chosen by Nozad et al. [START_REF] Nozad | Heat conduction in multiphase systems-II[END_REF] and their solution was validated against eq.(4.98) [START_REF] Whitaker | The Method of Volume Averaging[END_REF].

For this geometry the boundary condition eq.

Square array of cylinders

Now in order to complete our study , another geometry of a two phase system is tested: a square array of cylinders. The corresponding unit cell is shown in Figure 4. [START_REF] Angelini | Premixing-related behavior of steam explosions[END_REF]. We choose to solve the closure problem on 1/4 th of the geometry as Nozad et al. The conductivity ratio parametric study for ε β = 0.64 (Figure 4.12) shows also a perfect agreement between the numerical and analytical solutions. & (4.121)) the remaining terms can be simplified, at the end we obtain [START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF]:

• Problem I 0 = ∇ 2 b β - 1 V β A βγ n βγ . ∇b β dA βγ - 1 V β A βσ n βσ . ∇b β dA βσ (4.101) in β -phase 0 = ∇ 2 b γ - 1 V γ A γβ n γβ . ∇b γ dA γβ - 1 V γ
v mm + v mσ + v σσ + v σm = 1 V A βσ n βσ k β .∇b β dA βσ + 1 V Aγσ n γσ k γ . ∇b γ dA γσ + 1 V A βσ n βσ k β . ∇c β dA βσ + 1 V Aγσ n γσ k γ . ∇c γ dA γσ (4.133) + 1 V A σβ n σβ k σ . ∇c σ dA σβ + 1 V Aσγ n σγ k σ . ∇c σ dA σγ + 1 V A σβ n σβ k σ .∇b σ dA σβ + 1 V Aσγ n σγ k σ .
v mm + v mσ + v σσ + v σm = 0 (4.134)
Another important simplification is obtained using the boundary conditions (eqs. (4.129) & (4.131)), where it is found that the heat transfer coefficient coupling the 2 energy equations is unique:

h m = -h σ = a v h (4.135)
where a v is the specific area (1/m), and h is the heat transfer coefficient.

Numerical results

Let us consider the unit cell geometry shown in Figure 4. 14 . In section 4. Therefore two parametric studies have been performed. First the volume metal liquid fraction ε β is varied between 0.1 and 0.8 for a fixed value of the conductivity ratio between the liquid and solid metal (κ 1 = 2.96) and a fixed value of the conductivity ratio between the solid metal and the mould (κ 2 = 7.17 e -04 ) (calculated using the properties of the aluminium and Ωsphere mould Table 3.1). Secondly, for a fixed value of the mould volume fraction ε σ and for κ 1 = 2.96, the conductivity ratio between the mould and the metal κ 2 is varied between 10 -4 and 10 4. [START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF] shows that K mm varies quasi-linearly as a function of ε β , and that when ε β → 1, K mm approaches k β . The linear relation is due to the linear variation of the first term of the RHS of eq. (4.82) as a function of ε β that dominate the variation of the remaining terms. As expected the mould effective conductivity K σσ is found to be independent of ε β for a fixed value of ε σ , having the same order of magnitude as the mould conductivity k σ .

On the other hand, Figure 4. [START_REF] Calmidi | Forced Convection in High Porosity Metal Foams[END_REF] shows that the conductive cross coefficients K mσ = K σm are negligible compared to K mm . This result was expected since K mσ is one order of magnitude smaller than min(K mm , K σσ ) [START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF]. As stated earlier K σσ ∼ k σ while the Ω-sphere sand conductivity (k σ = 0.17 W/m 2 K) is small compared to the metal conductivity.

Conductivity ratio κ 2 parametric study large and small values of the mould/metal conductivity ratio (κ 2 < 10 -2 , κ 2 > 10 2 ). This can be explained by the fact that in the two extreme scenarios the mould becomes isothermal and the heat transfer will be independent from the mould conductivity and therefore of the conductivity ratio. It is also noticeable that for κ 2 > 10 2 , the asymptotic value increases with the mould volume fraction ε σ . Figure 4. [START_REF] Chang | Multi-scale analysis of effective transport in periodic heterogenous media[END_REF] shows that the cross coefficients are not negligible for large values of κ 2 where the asymptotic value is one order less than the metal conductivity k β , but in our case the conductivity ratio between the mould and metal is κ 2 = k σ k γ ∼ O(10 -3 , 10 -4 ).

Thus these terms can be neglected from the energy equations (4.79) & (4.80). The mould effective conductivity K σσ is normalized using the mould conductivity k σ . The mould effective conductivity varies slightly as a function of κ 2 . However, K σσ approaches a smaller constant value for large values of the conductivity ratio κ 2 > 1 then for small value of κ 2 < 1. This can be explained by the fact that the increase κ 2 leads to the increase of mould conductivity k σ which is translated by a decrease of the ratio K σσ /k σ . To be noted that the asymptotic values increase with the mould volume fraction ε σ .

Reduced form of the macroscopic model

In section 4. [START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF].2, we showed that the summation of the pseudo convective coefficients is null.

This is validated numerically since the solution of the closure problems using the symmetric unit cell (Figure 4. [START_REF] Brent | The Enthalpy-Porosity Technique For Modeling Convection-Diffusion Phase Change : Application To The Melting Of A Pure Metal[END_REF] gives the following results ( [START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF]):

v mm = v mσ = v σσ = v σm = 0 (4.136)
Therefore under these circumstances, the pseudo-convective coefficients as well as the conductive cross coefficients are neglected and the macroscopic problem can now be written in a reduced form as follows:

∂ ∂t ((ρc p ) m T m ) + ρ β L f ∂ε β ∂t = ∇ • (K mm .∇ T m ) -a v h ( T m -T σ σ ) (4.137) ∂ ∂t ρ σ ε σ c p σ T σ σ = ∇ • (K σσ .∇ T σ σ ) -a v h ( T σ σ -T m ) (4.138)
The simplified problem consists of a classical two energy equations coupled by the heat transfer coefficient h. Where the metal phase change is accounted for in the energy equation of the metallic mixture (second term of the LHS eq. (4.137)). Now all the effective properties are calculated, the closed macroscopic problem can be solved numerically in order to determine the temperature fields T m and T σ σ as well as the time evolution of the liquid metal volumetric fraction ε β and the solid-liquid interface location, this is presented in the next section.

Numerical resolution of the Macro problem

Let us recall the averaged equations of the solidification of a metal inside a porous medium:

∂ ∂t ((ρc p ) m T m ) + ρ β L f ∂ε β ∂t -v mm .∇ T m -v mσ .∇ T σ σ (4.139) = ∇ • [K mm .∇ T m + K mσ .∇ T σ σ ] -h m ( T m -T σ σ )
where

(ρc p ) m = ρ β ε β c p β + ρ γ ε γ c p γ ( 4
∂ ∂t ρ σ ε σ c p σ T σ σ -v σm .∇ T m -v σσ .∇ T σ σ (4.141) = ∇ • [K σσ .∇ T σ σ + K σm .∇ T m ] -h σ ( T σ σ -T m )
The metal foams are characterized by their high porosity values (ε σ = 0.7 -0.9). Thus, the value ε σ = 0.

Numerical accuracy

The numerical solution of the macroscopic model is performed using the multiphysics software Comsol 4.4 (on a 1,7 Ghz frequency, 6 Go RAM personal computer). In order to verify that the solution is independent of the mesh, a mesh refinement study is done where the maximum element size ∆x is changed for a range of values from 0.002 m to 0.006 m. The chosen criterion to determine the mesh independent solution is the initial calculated liquid metal fraction ε β value. Initially the mould porosity is taken to be ε σ = 0.7, and since the metal is at T = 1000 K higher than its fusion temperature then all the metal volume is in the liquid phase and ε β should be equal to 0.3. Figure 4. [START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a horizontal cylinder[END_REF] shows that the calculated initial value of metal liquid fraction becomes independent of the mesh size for ∆x < 0.004 m. Knowing that the calculation time increases when mesh size decreases the value of ∆x = 0.003 m was chosen as the reference mesh size in the scope of this thesis, but another dynamic mesh refinement study should be done using a more powerful calculation machine. is noticeable that the temperature profiles at all times present a slope rupture at the metal fusion temperature T = 933 K. This is due to the latent heat that need to be extracted from the metal in order to solidify, once it is evacuated the metal cools more rapidly as shown by the greater slope for T < T f usion .

Figure 4. [START_REF] Epstein | Complex Freezing-Melting Interfaces in Fluid Flow[END_REF] shows that the variation of the interface position as a function of t 1/2 is quasi linear which agrees with the literature (section 2.2.2) where it was found that the analytical and empirical solutions consider a linear relation between the solidification front position and the square root of time. Now in order to validate our estimation of the solidification front velocity, it is plotted as a function of time in Figure 4. [START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF]. First, due to the important temperature gradient at the wall (x = 0 m) the solidification front velocity is large (v ∼ 27 mm/s), but rapidly the solidification slows down, and after 4 s the front velocity is nearly constant and fluctuates around (v ∼ 5 mm/s). In section 4.4 the liquid-solid interfacial velocity was estimated to have an order of magnitude of O(10 -2 m/s) and was neglected, this estimation is now validated numerically since it is found that the solid/liquid interfacial velocity is O(10 Therefore, hereby another thermal configuration more adequate to the experimental conditions is presented in Figure 4. [START_REF] Fukusako | Solidification of Pure Liquids and Liquid Mixtures Inside Ducts and Over External Bodies[END_REF] where the 4 walls exchange heat with the ambient air through a heat exchange coefficient of h = 25 W/(m 2 K). In this configuration, the mould volume fraction is taken to be ε σ = 0.85, and initially the liquid metal is superheated to T m0 = 1133 K while the mould is kept at ambient temperature T σ σ 0 = 293 K.

The time temperature evolution at the cavity wall is presented in Figure 4. [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]. As we can see, at first the temperature decreases quasi-linearly with time, then at t ∼ 85 s metal phase change starts and the plot shows a temperature plateau until t sol ∼ 300 s, the solidification time of the whole metal domain. For t > t sol the temperature decreases again linearly as a function of time. To be noted that the cooling rate at the beginning of the cooling process in the liquid phase (t < 85 s) is greater than in the solid phase (t > t sol ). This can be explained by the fact that the mould temperature T σ σ has increased after the metal solidification process leading to a decrease of the heat exchange between the metallic mixture and the mould. [START_REF] Guyon | Hydrodynamique physique[END_REF]. The solidification front velocity is almost constant at the beginning of the process and fluctuates around 0.2 mm/s, then after about 250 s the solidification accelerates until complete solidification. This acceleration is due to the small remaining liquid metal volume that need to be solidified. Once again the assumption of neglecting the interfacial velocity in the closure problems (section 4.4) is validated as Figure 4. [START_REF] Guyon | Hydrodynamique physique[END_REF] shows that the solidification front velocity is about 0.15 mm/s much lower than its estimation << O(10 -2 m/s).

Comparison with experiments

The new metal foam manufacturing process via casting consists in infiltrating a sand preform by liquid metal. The goal is to fill the whole mould before solidification starts to avoid defects. The sand preform is formed by sand layers of spherical cells stacked in a cube. Each layer of sand beads is formed separately using the core shown in Figure 4.35-a, then the layers are stacked one over the other to form the sand preform. A schematic of the casting process is shown in Figure 4.36-a: liquid metal, heated to a predefined temperature (here T pouring = 1133 K), is poured in the sand preform by gravity, the obtained metal foam is shown in For the same thermal and morphological conditions, the deviation from one test to another is negligible. However, when the metal layer thickness is changed, so is the solidification time and the duration of the temperature plateau where the latent heat is evacuated. For a metal foam having a 2.5 mm thick metal layer, the solidification starts at about 10 s and finishes at 125 s.

On the other hand, the solidification of the metal foam with metal layer of The effective properties K σσ , a v h are calculated for ε σ = 0.85, while K mm is kept written as a function of ε β . The heat exchange coefficient between the equivalent domain and the ambient is taken to be h = 25 W/(m 2 K) corresponding to natural convection by air. The metal initial temperature T m is taken to be equal to the tests pouring temperature 1133 K, while the mould and air are at ambient temperature 293 K.

The temperature evolution as a function of time at the domain wall represented in Figure 4. [START_REF] Helmig | Efficient modelling of flow and transport in porous media using multi-physics and multi-scale approaches[END_REF] is compared to the temperature measurements in the metal layer of 10 mm thickness. The measured temperature shows that initially the mould is at ambient temperature and when the metal is poured into the preform, the temperature increases but does not reach the given pouring temperature T = 1133 K since the metal has already lost some of its sensible heat and instead T max ∼ 1090 K. The metal then cools down and loose its sensible heat before reaching the fusion temperature where the temperature remains constant as a function of time till complete solidification. The parabolic temperature decrease may be the result of the filling process since the metal volume is not constant as a function of time and the cooling is due to both diffusion and convection in the liquid metal. In our developed macroscopic model, infiltration and fluid flow were not taken into consideration, since initially the mould is considered to be completely filled by liquid metal at rest thus heat transfer occurs only by pure diffusion. Therefore the model validation zone is for t > 80 s, and the main difference between the numerical result and the tests is the solidification starting time and initial cooling rate. However, the final solidification time (here ∼ 280 s) and the cooling slope in the solid phase once the solidification is completed fit perfectly with the measurements.

Conclusions

In this chapter, a macroscopic model of the solidification of a molten metal in a homogeneous porous medium is developed. The model is obtained by applying the volume averaging method on the conservation equations at local scale, LTE is considered between the liquid and solid phase of the metal which is logical since the phase change takes place at the metal fusion The metal foam manufacturing process consists in pouring liquid metal in a porous mould, the goal is to complete the filling process before the solidification of the metal starts to avoid defects. Treating the whole problem will result in dealing with 4 phases model; liquid and solid metal, mould and air. Therefore, in order to simplify the problem the mould is considered to be initially fully saturated by liquid metal at rest, thus the metal solidification is by pure diffusion (no convection). Note that the infiltration and fluid flow were not taken into consideration in developing the macroscopic model.

Numerical resolution of the associated closure problems were done in a representative unit cell figuring the 3 phases (β, γ and σ). The solution in the symmetric unit cell showed that the cross conductive coefficients K σm = K mσ = 0 and the pseudo-convective coefficients v ii = v ij = 0 can be neglected. Therefore, the macroscopic model is written in a reduced form and its numerical resolution showed that the solidification front position varies linearly as a function of √ t, and that the interfacial velocity is ∼ O(10 -3 , 10 -4 ) which validates our assumption in neglecting this term in the closure problems.

A proper numerical validation of the macroscopic model of the diffusion solidification of liquid metal at rest inside a porous mould against metal foam casting tests of a 100 × 100 × 100mm 3 cubes having 85% porosity could not be done because there is a lot of uncertainty around the tests conditions. Even though, the numerical and experimental temperature profiles as a function of time showed some similarities such as the total solidification time and the cooling rate once the solidification is completed.

As stated before the problem was simplified thus the infiltration mechanism that implicates a

Conclusions

The objective of this work is to elaborate a new manufacturing process of metal foams via casting. The process consists in pouring by gravity liquid metal inside a porous sand preform, the liquid must infiltrates the mould before phase change starts to avoid defects. To optimize this process, its operating parameters must be controlled. Thus the aim of this study is to model the infiltration and solidification of liquid metal inside a porous medium. However, due to the complexity of this problem we chose to divide the study into two steps. First, at local scale one strut of the metal foam is considered as a capillary tube and the infiltration and solidification of liquid metal inside a cylindrical mould is studied. Second, a macroscopic model of diffusive solidification is derived using the volume average method.

The local model of infiltration and solidification of liquid metal inside a cylindrical mould, is formulated based on a one domain method and makes use of the volume average method (VOF) to determine and track the liquid/air interface. The enthalpy porosity method is used for modelling the metal phase change and two penalty coefficients were added to the momentum conservation equation to ensure zero velocity solutions in the mould and solidified metal. The four phases present in this process are identified using phase indicators and volume fractions, the effective properties of the equivalent domain are calculated as weighted average using these fractions.

The model is coded in an open source CFD tool (OpenFoam), it was validated against well known analytical solutions (Poiseuille's flow, internal forced convection and inward cylindrical diffusive solidification). Then, the model is used to simulate the coupled problem of infiltration and solidification inside a cylindrical mould. The numerical solution gives a qualitative observation of the infiltration and solidification dynamics: first, the liquid infiltrates the mould and pushes out the air initially present. Then, the liquid flow accelerates and the liquid-air interface is tracked as function of time, until the metal looses all its sensible heat and starts to solidify first at the mould/ metal interface. The solid crust grows and the mould section decreases, then the metal/air interface solidifies as well and the flow is stopped. The solidification proceeds by pure diffusion after flow cessation. The sensitivity of the infiltration time and length on the operating conditions is tested. Three parametric studies prove that the infiltration length (and subsequently time) increases linearly as function of the metal superheat (agrees with literature) and the initial mould temperature. However, the infiltration length increases as function of the pressure difference while the infiltration time decreases which agrees with results found in literature.

The chosen modelling approach of considering a one domain method presents some limitations. The presence of these currents depends on the problem configuration (La > 10 6 ), to minimise their effects a small time step must be used which slows down drastically the numerical calculations. This numerical problem is well known in the CFD community and is not related to the VOF method but to the widely used CSF model. Second, in the aim to spear the liquid/air interface and create a discontinuity in the calculated effective properties a numerical error is found. Actually in some calculation cells near the mould/metal interface the effective properties were those of the air instead of the metal. This error led to an underestimation of the heat transfer between the mould and metal. However, since the miscalculation is located on one or two calculation cells at some location along the mould length, we believe that the effect of this numerical error is not drastic and that the results are qualitatively correct.

The numerical problem that causes the miscalculation of the effective properties is now being solved, the solution of the current simulations will permit the identification of the effect of this error on the results presented in this thesis. However, the author believe that the use of twodomain approach in treating the fluid and mould will be more accurate. Since the contact angle can be imposed directly at the mould/fluid boundary and only the energy conservation equation will be solved in the mould, which will accelerate the numerical calculation. This option was not available using the chosen CFD tool (OpenFoam) where the solid fluid contact is not treated.

However, many existing commercial CFD tools treat the coupled solid/liquid problem.

The second step of our modelling approach is to develop a macroscopic model of the infiltration and solidification of liquid metal inside a porous medium. Such model requires dealing with four phases (the solid and liquid metal, the mould and the air), where the solid/liquid and liquid/air interfaces evolve with time. Therefore we chose to simplify the problem by considering that the liquid metal infiltrated the mould before solidification starts. Thus the porous medium (mould) is saturated by liquid metal at rest and solidification occurs by pure diffusion (no convection).

Meaning that the flow and infiltration processes are not taken into account. The numerical solution of the macroscopic model on a 100 mm 2 cavity showed that the solidification front position varies quasi-linearly as function of the square root of time which agrees with literature, and that the solidification front velocity is of the order of O(10 -3 , 10 -4 m/s).

A proper validation against experimental measurements could not be done because several test conditions are unknown to the author, however some similarities were detected in the time temperature evolutions.

methodAbstract:

  The objective of this work is to elaborate a new manufacturing process of metal foams via casting by modelling the infiltration and solidification of liquid metal inside a porous medium. However, due to the complexity of this problem the study is divided into two steps. First, at local scale one strut of the metal foam is considered as a capillary tube and the infiltration and solidification of liquid metal inside a cylindrical mould is studied. Second, a macroscopic model of diffusive solidification is derived using the volume average method. The local model is coded in an open source CFD tool and three parametric studies were done where the relations between the infiltration length and time as function of the operating parameters are determined. The modelling of the solidification of liquid metal inside a porous medium is simplified by considering that the mould is fully saturated by liquid metal at rest, solidification occurs by pure diffusion. Local thermal equilibrium (LTE) is considered between the solid and liquid phases of the metal while local thermal non equilibrium (LTNE) is retained between the metallic mixture and the mould. The associated closure problems as well as the macroscopic problem were numerically solved. Titre: Modélisation et Simulation Multi-échelle du procédé de fabrication des mousses métallique par voie de fonderie Mots clés: mécanique des fluides, solidification, simulation numérique, milieu poreux, prise de moyenne volumique Résumé: L'objectif est d'élaborer un nouveau procédé de fabrication de mousses métalliques par voie de fonderie en modélisant l'infiltration et la solidification d'un métal liquide dans un milieu poreux. La modélisation est faite en deux étapes. Tout d'abord, à l'échelle locale un brin de la mousse métallique est considéré comme un tube capillaire et l'infiltration et solidification d'un métal liquide dans un moule cylindrique est étudiée. Deuxèmement, le modèle macroscopique de la solidification diffusive d'un métal liquide dans un milieu poreux est obtenu par prise de moyenne volumique. Le modèle local est codée dans un outil CFD open source et trois études paramétriques ont été faites permettant la détermination des relations de la longueur et le temps d'infiltration en fonction de paramètres de fonctionnement. La modélisation de la solidification d'un métal liquide dans un mileu poreux est simplifié en considérant que le moule est complètement saturé par un métal liquide au repos, par suite la solidification se produit par diffusion pure (pas de convection). L'équilibre thermique local (LTE) est considéré entre les phases solide et liquide du métal tandis qu'un non équilibre thermique local (LTNE) est retenue entre la phase métallique et le moule. Les problèmes de fermeture associés ainsi que le problème macroscopique ont été résolus numériquement. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE I would like to express my deep gratitude to my supervisor Professor Benoît Goyeau for his guidance, valuable scientific participation and useful critics of this research work. But, most I would like to thank him for his moral support during the ups and downs of my thesis years.

1 Figure 1 . 1 :Figure 1 . 2 :

 11112 Figure 1.1: An example of metallic foam

4 2 State of the art 2 . 1

 221 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE Chapter At the pore scale

Figure 2 . 1 :

 21 Figure 2.1: Scheme of fluidity tests, a) spiral test, b) vacuum test[START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF]

Figure 2 . 2 :

 22 Figure 2.2: Solidification of a pure metal with no superheat, a) solidification starts at the mould entrance, b) solid grows as the metal continues to flow, c) "chock" at the entrance, cessation of the flow[START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF]

Figure 2 . 3 :

 23 Figure 2.3: Solidification of a superheated pure metal, a) liquid flows in the channel and looses its superheat, b) solidification starts downstream, c) "chock" near the metal-gas interface region, cessation of the flow[START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF]

Figure 2 .

 2 Figure 2.4 illustrates the solidification mechanism of a dilute alloy with no superheat. In case of metal superheat (see Figure 2.5), coarse grains are detected at the channel entrance where the solidification occurs after the cessation of the flow.

Figure 2 .

 2 Figure 2.5 shows the columnar grains obtained from the solidification of a dilute alloy ( Al 99.99%) we can notice that the columnar grains point upstream which is a particularity of columnar grains growing in a fluid flow [29].

Figure 2 . 4 :

 24 Figure 2.4: Solidification of a dilute alloy with no superheat, a) columnar grains start to grow at the channel entrance instantly, b) columnar grains continue to grow as the fluid flow, c) "chocking" occurs at the tube entrance even though the section is not completely solidified[START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF]

Figure 2 .

 2 Figure 2.5: longitudinal macrostructure of a pure aluminium (99.99%) poured with 83 • C superheat [29]

Figure 2 . 6 :

 26 Figure 2.6: Flow and solidification schematic of a rich alloy, a) nucleation starts at the flow tip, nucleation continue while the fluid flows inside the channel, c) flow ceases due to high friction at the flow tip[START_REF] Niesse | Application of theory in understanding fluidity of metals[END_REF]

Figure 2 . 7 :

 27 Figure 2.7: longitudinal macrostructure of an aluminum alloy (Al-5% Sn) poured with 82 • C superheat [29]

Figure 2 . 8 :

 28 Figure 2.8: Melt or pouring temperature effect on the Al-4.5%Cu fluidity for two different metal heads[START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF]

Figure 2 . 9 :

 29 Figure2.9: Fluidity of Al-4.5%Cu as function of metal head[START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF]

Figure 2 . 10 :

 210 Figure 2.10: Fluidity of Al-4.5%Cu in sand moulds. a-effect of moulding materials (1-plaster, 2-green silica, 3-green zircon), b-effect of the sand bond (4-bentonite, 5-linseed oil, 6-sodium silicate), c-effect of mould temperature, d-effect of mould coating (7-hexachloroethane, 8carbon black, 9-uncoated)[START_REF] Flemings | Fluidity of metals: Techniques for Producing Ultra Thin Section Castings[END_REF]

Figure 2 . 11 :

 211 Figure 2.11: Effect of metal mould temperature on casting fluidity of Al and its alloys on cast iron mould [65]

Figure 2 . 12 :

 212 Figure 2.12: Solidification of a liquid flow : the academic problem

  Few studies coupled the energy and hydrodynamic aspect of the problem and solved numerically the mass, momentum and energy conservation equations of the solidification of a fluid flow inside a tube, UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  .1). The sand mould is formed of Ω -sphere beads (special casting sand) having a low conductivity to increase the fluidity of the metal. The casting operation conditions were maintained constant throughout the experiments ( metallo-static height, section dimension, ...). To measure the local flow speed and estimate the infiltration time, Cronofond contacts sensors are implemented every 40 mm along the spiral, and the electric signal at the contact with liquid metal triggers the time measurement with an uncertainty of +/ -2 ms.

Figure 3 . 1 :

 31 Figure 3.1: Fluidity test in a spiral mould

Figure 3 . 2 :

 32 Figure 3.2: Flown distance as function of time

Figure 3 . 3 :

 33 Figure 3.3: Specimen 1, macro-structure of pure aluminium (99.5%) cast in a silica mould; flow direction from left to right and top to bottom.

Figure 3 . 4 :

 34 Figure 3.4: Specimen 2, macro-structure of pure aluminium (99.5%) cast in Ω -sphere mould; flow direction from left to right and top to bottom.

  5) As we can see the expression above has a non zero value only in the transition region where the volume fraction α has a value between ]0, 1[, thus the surface tension force acts only within the interface region. The CFS model presented above is only valid when the assumption of a constant surface tension is made, if not eq.(3.5) should contain an additional term which figures the surface tension gradient. For simplification, we will consider that the surface tension is 32 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE independent of the temperature and that eq. (3.5) is valid. The term κ in eq.(3.5) represents the interface curvature, calculated from κ = -∇. n = -∇. ∇α |∇α| (3.6)

  mentum equation in order to insure a null velocity in the solidified metal, therefore the mo-UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE mentum conservation equation writes: ∂(ρv) ∂t + ∇.(ρ v v) = -∇p + ∇.(µ∇v) + Av (3.11) where A has infinite value in the solidified metal and zero value in the liquid domain so the the typical momentum conservation equation can be recovered. It can have the Carman-Kozeny formulation in the mushy zone for alloy solidification or a simple linear function in case of front place solidification [14].

1 . 2 .

 12 φ the mould indicator: α the metal/air fraction: φ = 1, α = 0 air 1 metal 3. g l the liquid metal fraction:

  3.4) is implemented in a CFD open source software OpenFOAM[5]. The basic source code is written in C++ for Linux operating systems and contains the main modules for solving flow and heat transfer problems using the classical finite volume method[START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF].The numerical method is well-known and the solution procedure of the transient equations can be summarized as follows (see fig.3.5): first a computational mesh is generated and the variable fields are initialized. Then, the time loop starts, the time step is calculated in order to satisfy the Courant-Friedrichs-Lewy (CFL) condition. The transport equation of the metal volume fraction eq.(3.20) is solved first, which enables us to calculate the surface tension force. Subsequently, the metal-gas interface is built-up and the physical properties are updated as functions of the phase fractions. Next, the pressure-velocity coupling is solved using the PIMPLE scheme (a combination of the PISO and SIMPLE algorithms) of OpenFOAM for the continuity and the momentum equations(3.18, 3.19). The convergence of the solution is monitored by calculating the error on the velocity divergence. At this stage, the pressure and velocity fields as well as the metal fraction distribution are known 40 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE for the current time step. Using the velocity field solution, the energy equation (3.24) can now be solved to obtain the temperature field and the liquid metal fraction g l , which is used to update the penalty source term of the momentum equation. Now all the variable fields are known, if the final time is not reached yet, the time loop continues with the calculation of the next time step. The program runs on a 24 processors machine (3,47 Ghz × 19 frequency, 63 Go RAM), the operating system is Ubuntu 12.04 LTS. A typical calculation time on a 6000 x 70 grid and a timestep of 10 -5 s for solidification is of the order of an hour, this time increases drastically when treating infiltration and solidification problem where the calculation time is of the order of couple of weeks. Numerical validation tests have been performed separately to validate the fluid flow and the phase change problems and to verify grid convergence. The simulations presented hereafter have been performed in 2 D, since the natural convection and gravity effects were neglected which make the problem axisymmetric.

Figure 3 . 5 :

 35 Figure 3.5: Program flow chart

  (3.32) the normal vector n at the liquid-air-mould interface is deduced byUNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCEthe rotation of the geometric face normal vector S f by an angle Ω. The defined vectors and corresponding angles at the liquid-gas-mould contact are shown in Figure3.[START_REF] Bejan | Convection Heat Transfer[END_REF]. In 2 D, the angle between the geometric face normal S f and the normal vector n can be obtained from the imposed contact angle θ from the following simple expression

3.5. 1

 1 Convection problem3.5.1.1 Poiseuille's flowThe first step of our code validation is dedicated to simulating an isothermal one-phase flow and comparing it to the Poiseuille solution. The tube is initially filled with liquid aluminium at a temperature higher than its fusion temperature and a zero thermal flux is set at the tube wall. A 44 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE Dirichlet boundary condition of constant uniform velocity profile is imposed at the inlet. The dimensionless Reynolds number is taken equal to 600, which corresponds to an average velocity of 0.1 m/s. The entrance length after which the flow is fully developed is calculated as follows: L m = 0.0576 Re D = 82 mm (3.36) The radial velocity profile in a cross-section beyond the entrance length (x = 100 mm) is in good agreement with the Poiseuille solution, and the relative error on the maximum axial velocity is less than 1% (see fig. 3.7) .

Figure 3 . 7 :

 37 Figure 3.7: Velocity profiles at x = 100 mm

Figure 3 . 8 :

 38 Figure 3.8: Maximum velocity value at the tube axis at t = 0.18 s

Figure 3 .

 3 Figure3.12 shows that the displacing fluid flows forward at the tube centre, and then rolls out at the interface vicinity to flow on backward near the wall. On the other hand, the displaced fluid flows backward near the wall, and rolls in at the interface to flow forward at the tube centre.

•Figure 3 . 11 :Figure 3 . 12 :

 311312 Figure 3.10: Radial velocity profile at x = 60 mm and t = 0.18 s (liquid-air interface)

Figure 3 . 13 :

 313 Figure 3.13: Temperature radial profile at x = 100 mm

Figure 3 .

 3 Figure 3.14: Solid-liquid interface position as function of time

Figure 3 . 15 :Figure 3 . 16 :

 315316 Figure 3.15: Temperature profiles at t = 0.3 s

Figure

  Figure3.16 shows that in case of superheat the solid-liquid interface evolution as function of

Figure 3 . 17 :

 317 Figure 3.17: Solidification process of liquid aluminium at T = 975 K (liquid phase in red, solid phase in orange) initially the tube is completely filled with liquid aluminium, a. (t = 0.07s)aluminium solidification starts at the tube wall at a distance from the entry, b. (t = 0.5s), c. (t = 0.6s) -the solidified aluminium progresses in the radial direction, d. (t = 0.635s) -the solidified aluminium obturates the tube at a distance from the inlet and stops the fluid flow e. (t = 0.65s), f. (t = 0.7s)solidification proceeds by diffusion towards the entry, g. (t = 0.8s)complete solidification is achieved.

Figure 3 .

 3 Figure 3.18: liquid-solid interface radial position and Reynolds number as function of time at x = 25 mm

Figure 3 .

 3 Figure 3.19: Solid-liquid interface position as function of time for different solidification dynamics at x = 25 mm

Figure 3 .

 3 figuring in the definition of the penalty coefficient A. Figure3.20 shows that the evolution of the

  (a) at x = 25 mm (b) at x = 10 mm

Figure 3 . 20 :

 320 Figure 3.20: Maximum velocity at the tube axis as function of time

  [START_REF] Chen | Experimental and numerical study on melting of phase change materials in metal foams at pore scale[END_REF] for liquid Al (ρ = 2380 kg/m 3 ; σ = 0.85 N/m) flowing under pressure differences between 1400 and 1700 P a in a tube with an internal radius of 1.2 mm. The theoretical curves obtained for a contact angle θ = 150°are compared to the numerical simulations obtained with the code in the infiltration phase (before solidification). The comparison between the theoretical model and the numerical results shows a good agreement which confirm our assumption of neglecting the viscous forces compared to the inertial ones . 58 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 3 . 21 :

 321 Figure 3.21: Infiltration length and velocity

  3.2). The method consists in solving a transport equation of a volumetric fraction α. Contrary to a simple phase indicator that can have only two values 0 or 1, the resolution of α transport equation (eq. (3.20)) gives a distribution of values of the volumetric fraction α ∈ [0, 1] in the whole solution domain. This

3. 6 . 3 . 1

 631 Numerical accuracyIn section3.5, a mesh independent solution was easily found in case of flow in a cylindrical tube by comparing the numerical solution to Poiseuille's analytical solution. Another mesh refinement study should be done in case of flow in a cylindrical mould due to the discontinuity change of the porosity at the tube/mould interface. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 3 . 23 :

 323 Figure 3.23: Schematic of the calculation cell in the vicinity of the metal/air interface

Figure 3 . 24 :

 324 Figure 3.24: Liquid aluminium infiltrating a cylindrical mould; liquid phase (red), solid phase (orange), mould (blue), air (cyan)

Figure 3 . 25 :

 325 Figure 3.25: Details of the solidification zone; liquid phase (red), solid phase (orange), mould (blue), air (cyan)

Figure 3 . 26 :Figure 3 . 27 :Figure 3 . 28 :

 326327328 Figure 3.26: Isothermal lines map in the solidified metal at t = 0.3 s for T ∈ [932, 933[ K

Figure 3 . 29 :

 329 Figure 3.29: Temperature radial profile at x = 75 mm and t = 0.3 s

Figure 3 . 30 :

 330 Figure 3.30: Temperature radial profile at x = 75 mm and t = 0.3 s in the liquid and solid phases of the metal

  2 mm : 1100P a), but it cannot be too large because the liquid in the tube would reach large velocities that do not respect the laminar flow assumption. The first set of simulations has been realized with the previous conditions (∆T 0 = 112 K ; T w0 = 294 K). Four values of the imposed pressure difference are taken (1400, 1500, 1600 and 1700 P a) for the same thermal initial and boundary conditions. The results are displayed in terms of the infiltration length and time when flow blockage occurs in the tube.

Figure 3 . 32 :Figure 3 .

 3323 Figure 3.32: Influence of the imposed pressure difference on the infiltration length

Figure 3 . 33 :

 333 Figure 3.33: Influence of the imposed pressure difference on the infiltration time

  70 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE again this introduces one more arbitrary parameter and the present hypothesis makes sense in the context of this study. Under these assumptions, the code is now used to study the influence of one important process parameter, namely the injected liquid initial temperature, on the effective process indicators: the infiltration length and the infiltration time. We are using the same reference hydrodynamic conditions as above (∆p = 1400 P a), the mould and air are considered initially at ambient temperature and only the liquid metal initial superheat above the melting temperature is changed. The values of the injection temperature are those used in the fluidity tests (between 955 and 1155 K), corresponding to initial superheat between 22 K and 222 K receptively. The effect of this parameter on the infiltration length and infiltration time is shown in the plots of Figure 3.34.

Figure 3 . 34 :

 334 Figure 3.34: Influence of the liquid initial superheat on the infiltration length and time

  Figure 3.35: Influence of the mould initial temperature on the infiltration length and time

  The local model allows to test the sensitivity of the infiltration time and length of a liquid metal injected in a capillary mould and undergoing solidification on the operating conditions. A UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  three phases: the liquid and the solid metal and the mould. A volume averaging procedure is used in order to upscale the local conservation equations. The model is derived in the context of local thermal equilibrium assumption (LTE) between the liquid and the solid metallic phases, leading to one energy conservation equation for the equivalent metallic phase. However the local thermal non-equilibrium (LTNE) between the equivalent metallic phase and the mould is UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE retained, giving rise to two coupled energy conservation equations. The associated closure problems are derived and numerically solved allowing the determination of the effective transport properties. Afterwards, the two energy conservation equations are solved numerically in a representative equivalent domain, allowing the determination of the liquid-solid interface position evolution as well as the two temperature distributions in the metallic mixture and in the mould as a function of time.

4. 3 . 1

 31 Definitions and averaging theoremsLet us consider the local temperature T β associated to the β-phase. The superficial volume average of T β in the averaging volume V (Figure4.1) is given by :

Figure 4. 1 :

 1 Figure 4.1: Averaging volume of a 2-phases system

Figure 4 . 2 :

 42 Figure 4.2: Averaging volume of a 3-phases system

25 )

 25 Using equation (4.25) and introducing the Gray decomposition T β = T β β (r)+ T β in equation (4.22), gives

31 )

 31 Applying the same averaging procedure for the σ-phase leads to the following non-closed mac-UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE roscopic energy conservation equation

  ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE In addition, by taking into consideration eq.(4.25), the last term in equation (4.34) can be written under the form k=γ,σ

  58) where b i , c i and α i are the closure variables. The derivation of the associated closure problems UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE is presented in the next section.

  (4.79) & (4.80) represent the accumulation as a function of time, while the second term in the LHS of eq.(4.79) represents the latent heat evacuated by solidification at the solid-liquid interface. The last term in eqs.(4.79) & (4.80) represent the heat exchange that couple the two energy equations. The remaining terms contain both diffusive and

Figure 4 . 3 :

 43 Figure 4.3: Chang's unit cell

Figure 4 . 4 :

 44 Figure 4.4: Longitudinal component of the effective thermal conductivity tensor as a function of ε β , for κ = k σ k β = 10

Figure 4 . 5 :

 45 Figure 4.5: Longitudinal component of the effective thermal conductivity tensor as a function of the conductivities ratio κ = k σ k β for ε β = 0.64

Figure 4 .

 4 Figure 4.4 shows that for a constant value of κ > 1, the effective conductivity decreases in a parabolic way and approaches the fluid conductivity when ε β → 1.Figure 4.5 shows that for

  Figure 4.5 shows that for large values of κ (in other words k σ ), the effective conductivity is independent of the latter since the σ phase becomes isothermal and the temperature field in the β phase will depend only on k β .

  (4.97) is now replaced by a periodic boundary condition:b β (r + l j ) = b β (r), j = 1,2, 3 (4.100) 92 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4 . 6 :

 46 Figure 4.6: Square array of squares unit cell

  [START_REF] Nozad | Heat conduction in multiphase systems-II[END_REF] (Figure4.[START_REF] Bauer | A general analytical approach toward the thermal conductivity of porous media[END_REF], where the boundaries at y = 0 and y = L/4 are considered adiabatic and a Dirichlet boundary condition b β = 0 is applied at x = 0 and x = L/4 . This geometry presents some limitations because the volume liquid fraction cannot take a value less than 0.3, otherwise the cylinder diameter would exceed the square edge length.

Figure 4. 11 shows

 11 the effective conductivity variation as a function of ε β , the numerical result is in perfect agreement with the analytical solution for large values of the liquid volume fraction (ε β > 0.5).

  UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4 . 7 :

 47 Figure 4.7: Longitudinal component of the effective thermal conductivity tensor as a function of ε β , for κ = k σ k β = 10

Figure 4 . 8 :

 48 Figure 4.8: Longitudinal component of the effective thermal conductivity tensor as a function of the conductivities ratio κ = k σ k β for ε β = 0.64 and ε β = 0.19

Figure 4 . 9 :

 49 Figure 4.9: Square array of cylinders unit cell

Figure 4 . 10 :

 410 Figure 4.10: Square array of cylinders unit cell

Figure 4 . 11 :

 411 Figure 4.11: Effective thermal conductivity as a function of ε β

Figure 4 . 12 :

 412 Figure 4.12: Effective thermal conductivity as a function of the conductivity ratio κ = kσ k β for ε β = 0.64

Figure 4 . 13 :

 413 Figure 4.13: Effective thermal conductivity as a function of the conductivity ratio κ = kσ k β for ε β = 0.64 for two arrays with different particle geometries

4. 7 . 2

 72 Numerical resolution of the closure problems of heat conduction with solidification in a 3-phase system Let us recall and solve numerically the three closure problems of the 3-phase system (β,γ and σ) in order to evaluate the corresponding closure variables and the effective properties present in eqs. (4.79) & (4.80). Re-writing eqs.(4.59-4.76) for each phase (β,γ and σ) gives the following detailed closure problems:

  ∇b σ dA σγ Boundary conditions (eqs. (4.107) & (4.109)) lead to the simplification of the two first and last terms of the RHS of eq.(4.133). On the other hand using the boundary conditions (eqs. (4.119)

7 . 1

 71 we have shown that the geometry of the particle array hardly affects the solution. Thus to simplify the numerical solution, we choose this symmetrical unit cell instead of an array of spheres that may be more adequate since the sand preform used in the manufacturing of metal foams consists of an arrangement of Kelvin cells layers. 100 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCEIn order to simplify the problem, we consider that the σ-phase and the β-phase are never in contact. This means that initially the phase change has started and the solid phase is already formed at the contact with the mould. It is of interest to follow the evolution of the effective properties as a function of the conductivity ratio between the mould and metal, as well as the liquid metal volume fraction ε β , e.g. as a function of the metal phase change.

Figure 4 . 14 :

 414 Figure 4.14: Averaging volume for 3-phases solidification configuration

  4 .The solution of the closure problems for a 2 D of square arrays gives an isotropic effective thermal conductivity tensor and therefore only one component is necessary to determine the latter.Volume metal liquid fraction ε β parametric studyLet us first examine the effect of the variation of the liquid metal volume fraction ε β on the effective conductivities. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4 . 15 :

 415 Figure 4.15: Effective thermal conductivity as a function of the liquid metal volume fraction ε β , where κ 1 = 2.96, κ 2 = 7.17e -04 and ε σ = 0.1

Figure 4 . 16 :

 416 Figure 4.16: The effective conductivity K mσ = K σm as a function of the liquid metal volume fraction ε β , where κ 1 = 2.96, κ 2 = 7.17e -04 and ε σ = 0.1

Figure

  Figure4.[START_REF] Flemings | Fluidity of Aluminium alloys. An experimental and quantitative evaluation[END_REF] shows that K mm varies quasi-linearly as a function of ε β , and that when ε β → 1, K mm approaches k β . The linear relation is due to the linear variation of the first term of the

Figure 4 . 17 :

 417 Figure 4.17: Effective thermal conductivity as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1.

  Figure4.17 shows that the metal effective conductivity K mm is nearly independent of κ 2 especially for small values of ε σ . To be noted that the effective conductivity increases when the mould volume fraction ε σ decreases, in other words, when the metal volume fraction increases. This behaviour is expected since the metal effective conductivity is related to the metal volumetric fraction.

Figure 4 .

 4 Figure 4.18 shows that the conductive cross coefficients approach two asymptotic values for

Figure 4. 18 :

 18 Figure 4.18: Effective thermal conductivity K mσ = K σm as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1

Figure 4 . 19 :

 419 Figure 4.19: Log-Log plot of the effective thermal conductivity K mσ = K σm as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1

Figure 4. 20 :

 20 Figure 4.20: Effective thermal conductivity of the mould as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1

Figure 4. 21 :

 21 Figure 4.21: Effective heat transfer coefficient as a function of the metal-mould conductivities ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1

Figure 4 . 22 :

 422 Figure 4.22: Log-Log plot of the effective heat transfer coefficient as a function of the metalmould conductivities ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1

A 2 D

 2 cavity (Figure 4.23) of 0.1 m length is considered as the equivalent continuous medium of the porous medium and metallic mixture. Initially the metal is at T = 1000 K, higher than its fusion temperature (T = 933 K), while the mould is at the ambient temperature T = 293 K. A Dirichlet condition of T = 293 K is imposed at x = 0 m and the three remaining boundaries are considered adiabatic. Therefore the metal will solidify by loosing heat to the mould and to the cold wall.

Figure 4. 23 : 2 D

 232 Figure 4.23: 2 D square equivalent continuous medium

Figure 4 . 24 :

 424 Figure 4.24: Initial calculated value of the liquid volume fraction ε β as a function of 1/∆x and the corresponding relative error % for different mesh refinement values

  Figure 4.25 shows that the solution is nearly independent of the latter, therefore the value of ∆t = 10 -1 s is chosen for the future simulations. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4 . 25 :

 425 Figure 4.25: Time evolution of the liquid volume fraction ε β for three time step: ∆t = 10 -1 , 10 -2 , 10 -3 s

Figure 4 .

 4 Figure 4.27 shows the temperature profiles at y = 0.05 m for different time. At t = 0 s the initial and boundary conditions cause a discontinuity at x = 0 m where the temperature passes

Figure 4 . 26 :

 426 Figure 4.26: The time evolution of the liquid volume fraction ε β obtained using two models : (Model eqs. 4.79 & 4.80 and Simplified Model eqs. 4.137 & 4.138)

Figure 4 . 27 :

 427 Figure 4.27: Metallic temperature distribution at y = 0.05 m for different time

Figure 4 . 29 :

 429 Figure 4.29: Solidification front velocity as a function of time

Figure 4 .

 4 Figure 4.32 shows the time evolution of the metal liquid fraction, initially ε β = 0.15 after almost

Figure 4 . 32 :

 432 Figure 4.32: The evolution of the metal liquid fraction ε β as a function of time

  & 4.

  Figure 4.36-b where a metal layer is formed above the metal foam. The thickness of this layer can be controlled, and two series of casting tests were performed for metal thickness of 2.5 mm and 10 mm. The metal foams obtained from these casting tests have the following dimensions: 100 × 100 × 100 mm 3 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4 . 33 :

 433 Figure 4.33: Solidification front position as a function of the square root of time

Figure 4. 35 :

 35 Figure 4.35: Core used to create the sand layers formed of spherical cells

Figure 4. 37 :

 37 Figure 4.37: The evolution of temperature as a function of time in the metal layer of 2.5 mm thickness [3]

  temperature T β β = T γ γ = T f usion , and the thermo-physical properties of the liquid and solid metal are close. However, local thermal equilibrium cannot be considered between the mould and the metal mixture since T f usion >> T σ σ and there is a high contrast between their thermophysical properties. The originality of this work lies in considering LTNE between the metallic mixture and the mould, giving rise to two energy conservation equations coupled by an effective heat transfer coefficient a v h calculated by solving the associated closure problems. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Figure 4. 39 :

 39 Figure 4.39: The evolution of temperature as a function of time at the equivalent domain wall (Numerical solution) and in the metal layer of 10 mm thickness (tests 4, 8, 12 and 16)

  120 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE forth phase (air) was not taken into consideration. Therefore, future work should be dedicated to the treatment of the coupled infiltration and solidification process of liquid metal inside a porous medium. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE Chapter 5

First

  , the use of the CSF model to model the surface tension is found to create spurious currents. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE 123

  figure in the reduced form of the macroscopic model.

  The diffusive solidification macroscopic model is a first step in modelling the manufacturing process of metal foams. Future work should be dedicated to the derivation of a macroscopic model of diffusive and convective solidification coupled to the infiltration taken into account the presence of the four phases. Such a model can be obtained by applying the average operator on the mass, momentum and energy equations. However, treating the interfaces and their time evolution is a challenging task left for the reader. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE 3.31 Local thermal flux at x = 75 mm and t = 0.3 s in the vicinity of the metal/mould interface .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.32Influence of the imposed pressure difference on the infiltration length . . . . .3.33Influence of the imposed pressure difference on the infiltration time . . . . . .3.34 Influence of the liquid initial superheat on the infiltration length and time . .

  .

  .

3. 35

 35 Influence of the mould initial temperature on the infiltration length and time .

  .

4. 1 4. 8 4. 9

 189 Averaging volume of a 2-phases system . .. . . . . . . . . . . . . . . . . . .4.2 Averaging volume of a 3-phasessystem . . . . . . . . . . . . . . . . . . . . .4.3 Chang's unit cell .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.4 Longitudinal component of the effective thermal conductivity tensor as a function of ε β , for κ = k σ k β = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Longitudinal component of the effective thermal conductivity tensor as a function of the conductivities ratio κ = k σ k β for ε β = 0.64 . . . . . . . . . . . . . . 4.6 Square array of squares unit cell . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Longitudinal component of the effective thermal conductivity tensor as a function of ε β , for κ = k σ k β = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . Longitudinal component of the effective thermal conductivity tensor as a function of the conductivities ratio κ = k σ k β for ε β = 0.64 and ε β = 0.19 . . . . . . Square array of cylinders unit cell . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Square array of cylinders unit cell . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Effective thermal conductivity as a function of ε β . . . . . . . . . . . . . . . . 4.12 Effective thermal conductivity as a function of the conductivity ratio κ = kσ k β for ε β = 0.64. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.13 Effective thermal conductivity as a function of the conductivity ratio κ = kσ k β for ε β = 0.64 for two arrays with different particle geometries . . . . . . . . . . . 4.14 Averaging volume for 3-phases solidification configuration . . . . . . . . . . . 4.15 Effective thermal conductivity as a function of the liquid metal volume fraction ε β , where κ 1 = 2.96, κ 2 = 7.17e -04 and ε σ = 0.1 . . . . . . . . . . . . . . . 4.16 The effective conductivity K mσ = K σm as a function of the liquid metal volume fraction ε β , where κ 1 = 2.96, κ 2 = 7.17e -04 and ε σ = 0.1 . . . . . . . . . . 4.17 Effective thermal conductivity as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 4.18 Effective thermal conductivity K mσ = K σm as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE 4.19 Log-Log plot of the effective thermal conductivity K mσ = K σm as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1 . . . . . . . . . . . . . . . . . . . 4.20 Effective thermal conductivity of the mould as a function of the metal-mould conductivity ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Effective heat transfer coefficient as a function of the metal-mould conductivities ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Log-Log plot of the effective heat transfer coefficient as a function of the metalmould conductivities ratio κ 2 = k σ k γ for different mould volume fractions ε σ , where κ 1 = 2.96 and ε β = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 2 D square equivalent continuous medium . . . . . . . . . . . . . . . . . . . 4.24 Initial calculated value of the liquid volume fraction ε β as a function of 1/∆x and the corresponding relative error % for different mesh refinement values . .4.25 Time evolution of the liquid volume fraction ε β for three time step: ∆t = 10 -1 , 10 -2 , 10 -3s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.26 The time evolution of the liquid volume fraction ε β obtained using two models : (Model eqs. 4.79 & 4.80 and Simplified Model eqs. 4.137 & 4.138) . . . . . 4.27 Metallic temperature distribution at y = 0.05 m for different time . . . . . . . 4.28 Solidification front position as a function of √ t . . . . . . . . . . . . . . . . . 4.29 Solidification front velocity as a function of time . . . . . . . . . . . . . . . . 4.30 2 D cavity with a mixed boundary condition h = 25 W/(m 2 K) . . . . . . . . 4.31 Temperature evolution at the cavity wall as a function of time . . . . . . . . . . 4.32 The evolution of the metal liquid fraction ε β as a function of time . . . . . . . 4.33 Solidification front position as a function of the square root of time . . . . . . . 4.34 Solidification front velocity as a function of time . . . . . . . . . . . . . . . . 4.35 Core used to create the sand layers formed of spherical cells . . . . . . . . . . 4.36 Foam casting tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.37 The evolution of temperature as a function of time in the metal layer of 2.5 mm thickness [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.38 The evolution of temperature as a function of time in the metal layer of 10 mm thickness [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.39 The evolution of temperature as a function of time at the equivalent domain wall (Numerical solution) and in the metal layer of 10 mm thickness (tests 4, 8, 12 and 16) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Table 3 . 1 :

 31 Properties of the sand mould[START_REF] Dairon | Développement d'une nouvelle technique d' élaboration de mousses d'acier par fonderie et caractérisation mécanique[END_REF]

	Property	Silica	Ω -sphere
			(aluminosilicate: Al 2 O 3 /SiO 2 )
	Density (kg/m 3 )	1500	430
	Conductivity (W/mK)	0.86	0.17

Specific heat capacity (J/kgK) 996 780 24 UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Table 3 . 2 :

 32 Properties of liquid aluminium[START_REF] Hatch | Aluminum: Properties and Physical Metallurgy[END_REF] -sphere moulds were implemented with Cronofond contacts sensors because of the complexity of the implementation process. From which, one test failed to indicate the infiltration time due to some technical problems. Therefore, the infiltration lengths as function of the time is obtained for only one test. The results as well as the calculated local velocity and the Reynold's number calculated base on the hydraulic diameter (D H = 1.73 mm) are presented in Table3.3. We notice that the local Re number is less than 2400 thus the flow can be considered as laminar. The infiltration length varies quasi-linearly as function of time as shown in Figure3.2.

	Property	liquid	solid
	Density (kg/m 3 )	2380	2380
	kinematic viscosity (m 2 /s)	4.2 10 -7	-
	Conductivity (W/mK)	80	237
	Specific heat capacity (J/kgK)	1180	920
	Surface tension (N/m)	0.85	-
	Contact angle (°) Al-Silica	150	-
	Melting temperature (K)	933	933
	Latent heat (J/kg)	398 000	-
	Pr	0.015	-

UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE Only the Ω

Table 3 . 3

 33

	Re

: flown distance as function of time obtained from the fluidity test of liquid aluminium at T = 1123 K distance (mm) time (s) fluid mean velocity (m/s)

Table 3 .

 3 4: Experimental results of aluminum fluidity tests for various superheat tempeartures T pouring (K) ∆T superheat (K) Ste

	test	L(mm) time (s)

  Calculating the temporal derivative of the liquid fraction in order to replace it by its expression in eq. (3.24):

1 K is chosen. UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Table 3 .

 3 . 5: Radial mesh refinement for an axial cell size of ∆x = 0.05 mm

	Mesh quality	Mesh 1	Mesh 2	Mesh 3
		Tube mould Tube mould Tube mould
	∆r min	0.02	0.02	0.01	0.01	0.005 0.005
	∆r max	0.08	0.53	0.1	0.65	0.1	0.65
	Regression	1.06	1.11	1.078	1.14	1.09	1.15
	Cells / Radius	25	30	30	30	35	35

Table 3 .

 3 result, but due to lack of time this will not be presented in the scope of this thesis.

	6: Properties of air at ambient temperature
	Density (kg/m 3 )	1
	kinematic viscosity (m 2 /s)	1.48 10 -5
	Conductivity (W/mK)	0.0257
	Specific heat capacity (J/kgK)	1000
	Table 3.7: Properties of the equivalent fluid in the calculation cell
	Conductivity (W/mK)	0.0285
	Specific heat capacity (J/m 3 K) 281622
	3.6.3.2 Numerical result of the metal liquid infiltration and solidification inside a mould

UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE ing the numerical

  ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE n nσ . (k n ∇c n -k σ ∇c σ ) . I = n nσ . k σ I at A nσ , n = β, γ (4.69)

			67)
	c σ = c n	at A nσ , n = β, γ	(4.68)
	86	UNIVERSITÉ PARIS-SACLAY
		ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY
	ROUTE DE L'

  A βσ n βσ . k β ∇α β dA βσ + 1 V Aγσ n γσ . k γ ∇α γ dA γσ (4.85)

	h m =	1 V
						∇b γ dA γσ	(4.83)
	v mσ =	1 V A βσ	n βσ . k β ∇c β dA βσ +	1 V Aγσ	n γσ . k γ ∇c γ dA γσ	(4.84)
	88					UNIVERSITÉ PARIS-SACLAY
				ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY
	ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  7.2). Finally, h i are the heat transfer terms coupling the

	4.7 Numerical resolution of the closure problems
	4.7.1 Numerical validation: heat conduction without solidification in a two-phase
	system
	2 energy equations.
	UNIVERSITÉ PARIS-SACLAY
	ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY
	ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  Aγσ n γσ . ∇b γ dA γσ (4.102)

	α σ = α γ + 1	at A γσ	(4.130)
	n γσ . (∇α γ -κ 2 ∇α σ ) . I = 0	at A γσ	(4.131)
	α β (r + l j ) = α β (r)	for j = 1, 2, 3	(4.132)
	Using the boundary conditions cited above, some simplifications of the effective transport coef-
	ficients presented in eqs. (4.79) & (4.80) can be obtained. For instance, adding the pseudo
	convective transport coefficients gives:
	in γ -phase
	UNIVERSITÉ PARIS-SACLAY
	ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY
	ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  The heat transfer coefficient h is normalized byk β /a v l β 2 ,where k β is the metal conductivity, a v is the specific surface area and l β is the β-phase length ∼ O(10 -3 m). The evolution of a v hl β 2 /k b as a function of the conductivity ratio κ 2 is plotted in Figure4.[START_REF] Chen | Experimental and numerical study on melting of phase change materials in metal foams at pore scale[END_REF]. The heat transfer coefficient approaches asymptotic values for small and large values of the mould/metal

conductivity ratio (κ 2 < 10 -1 , κ 2 > 10 2 ). The asymptotic values increase with the mould volume fraction ε σ . To be noted that in our case where κ 2 ∼ O(10 -3 , 10 -4 ), the heat transfer UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

  7 is chosen as a reference value to validate our model and calculate the effective properties by solving the closure problems, as shown earlier in section 4.[START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF].2.1. The effective properties shown to be independent of the liquid metal volume fraction ε β (section 4.[START_REF] Alexiades | Mathematical Modeling Of Melting And Freezing Processes[END_REF].2.1), were evaluated for ε σ = 0.7 (K σσ , v σσ , v σm , v mm , a v h ). While the remaining effective properties (K mm , K σm ) are kept written as a function of the liquid metal volume fraction ε β as their values vary during the solidification process.

  10 mm thickness starts later at about 50 s and lasts for more than 200 s which is expected since the metal volume in this case is larger.

A lot of unknown conditions prevent us from comparing the test temperature measurements to UNIVERSITÉ PARIS-SACLAY ESPACE TECHNOLOGIQUE / IMMEUBLE DISCOVERY ROUTE DE L'ORME AUX MERISIERS RD 128 / 91190 SAINT-AUBIN, FRANCE

Acknowledgements

)

where