with my main contributions in green.

Édouard is someone you will instantly appreciate once you meet him. He is the kindest advisor and one of the smartest people I know. I hope that I could be as lively, enthusiastic and energetic as Édouard. Thank you Édouard for all the support, care and help... Thank you for being my family in France.

Special dedications to the SAGE team and the administration people at INRIA, specially for Fabienne Cuyollaa and Sophie Viaud.

My lab friends! I am grateful for all the nice moments we shared: Souhila, Mestapha, Louis-Bernard, Osama, Mona, Meriam and Shadi. Thank you Shadi for your care and the nice moments I spent with your family. You all made my stay in the lab nice and calm.

My lovely friends: Imen, Rana, Layal, Zainab and Nahla. Thanks for being always a source of love and support.

Rennes gave me two precious gifts: Areej and Marwen. You are indeed my guarding angels and the closest to my heart. Marwen, I can't express how grateful I am, we Les écoulements de uide et les transferts de chaleur et de masse dans les milieux poreux constituent un axe de recherche, théorique et appliquée, qui a eu une attention de plus en plus importante durant ces dernières années. Cela est dû à l'importance de ce domaine de recherche dans de nombreuses applications scientiques (chimiques, environnementaux, géologiques, mécaniques, ...). Une des applications du transfert de chaleur en milieu poreux est l'étude des feux préhistoriques an de comprendre le comportement des anciens humains ; ce sujet a été étudié par Mohamad Muhieddine dans sa thèse [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF]. Le présent ouvrage constitue une continuation de cette dernière thèse.

Le travail de M. Muhieddine [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] est constitué essentiellement de deux grandes parties. Dans la première, il a résolu l'équation de diusion de la chaleur lors de changement de phase (eau/glace) en utilisant la méthode de l'accumulation de la chaleur latente (LHA: Latent heat accumulation) et la méthode de la capacité thermique apparente (AHC: Apparent heat capacity). Une technique de maillage adaptatif a été appliqué sur la méthode LHA en une seule dimension (1D). Il a aussi étudié le couplage entre la diusion thermique et l'écoulement de la vapeur d'eau produite par l'évaporation en milieu poreux initialement saturé. Dans la deuxième partie, il a présenté une méthode inverse 1D dans le but d'estimer les propriétés thermo-physiques d'un sol saturé et subissant un changement de phase en se basant sur des données synthétiques.

Cadre General

Mon travail de thèse se situe dans un projet interdisciplinaire qui a commencé il y a plus de dix ans. Le projet ARPHYMAT [Cana] est consacré à l'étude des feux préhistoriques, en regroupant trois domaines complémentaires ; ARchéologie, PHYsique et MAThématiques. Il vise à estimer la durée des feux des anciens rassemblements hu- mains, an de comprendre les raisons derrière cette utilisation (cuisson, chauage, poteries, ...) et le comportement social de ces groupes.

L'ensemble des tâches (T ) qui décrivent ce projet sont présentés dans le diagramme de la page suivante.

Certaines expériences ont été réalisées (T 2 ) dans des situations réelles (en brûlant des morceaux de bois), tandis que d'autres en eu lieu en laboratoire (le feu a été remplacée par une plaque électrique chauante). L'objectif principal de ces expériences est la détermination des transferts de chaleur dans le sol an de comprendre le phénomène de séchage et de ses conséquences sur un sol réel. En outre, ces expériences nous permettent, en utilisant des méthodes appropriées (T 3 et T 5 (B)), d'estimer certaines propriétés Résumé Détaillé Estimation:

• Durée des feux.

• Propriétés thermo physiques du milieu poreux. Projet interdisciplinaire ARPHYMAT.

Étude des feux préhistoriques http://arphymat.univ-rennes1.fr/ Caractéristiques:

• Fort réchauffement du sol (milieu poreux).

• Existence de l'eau. Models physiques:

• models actuels: Expériences:

• Réplication (terrain).

• Expériences de laboratoire : determination de ρ𝐶 𝑒 , 𝜆 𝑒 (conductivité thermique effective).

Modèle mathématique (cas sec seulement):

• Laloy & Massard.

• Laplace Transform.

Diffusivité α

Méthodes Numériques & simulations:

• problème direct:  Calcul théorique de 𝜆 𝑒 :

• calcul exact de la forme du ménisque liquide entre deux grains sphériques. Un schéma récapitulatif représentant le cadre de ma thèse avec mes principales contributions en vert.

thermo-physiques du sol. Enn, elles nous aident à quantier l'inuence de la quantité d'eau existante dans le sol (rappelons que le changement de phase absorbe une partie très importante de l'énergie thermique).

En plus des expériences (T 2 ), des simulations numériques (T 5 ) sont nécessaires an de mieux comprendre l'inuence de certains paramètres. Ces simulations sont basées sur des modèles physiques (T 4 ).

Les codes numériques utilisés dans T 5 sont des codes de laboratoire, développés par É. Canot et M. Muhieddine à partir de 2009. La raison de ce choix est d'avoir un accès facile au code an de pouvoir modier les méthodes numériques et introduire de nouveaux modèles physiques qui sont parfois non conventionnels (nous verrons plus tard que, pour un milieu poreux contenant de faibles quantités d'eau liquide, la conductivité thermique eective présente un comportement d'hystérésis).

Les modèles physiques en T 4 présentent diérents degrés de diculté. Le modèle M 1 a été utilisé par É. Canot et M. Muhieddine depuis plusieurs années, mais il comporte plusieurs lacunes. Pour cette raison, le modèle M 2 a été développé pour traiter l'alternance des phases chauage/refroidissement avec changement de phase, mais sans utiliser la méthode AHC. Un troisième modèle plus réaliste (M 3 ) est proposé, mais il n'est pas encore implémenté.

Une autre partie du travail, T 6 , concerne la détermination de la conductivité thermique eective d'un milieu poreux humide. Bien que de nombreux résultats concernant le cas de sol sec puissent être trouvés dans la littérature, le cas de trois composantes (solide, eau et gaz) n'est pas bien décrit. C'est pourquoi nous avons fait un premier modèle simple qui montre un comportement d'hystérésis. T 6 permet le calcul de la conductivité thermique eective, qui à son tour est nécessaire pour le futur modèle M 3 .

Ma thèse est divisée principalement en six chapitres.

Chapitre 1

Le premier chapitre est une description générale d'un milieu granulaire où nous dénissons les propriétés physiques liées à notre étude : conductivité thermique (λ), la masse volumique (ρ), la capacité thermique (C ), la tension supercielle (forces capillaires), l'angle de contact d'une interface gaz-liquide en contact avec un solide et les schémas de répartition de l'eau dans un milieu granulaire en particulier le régime pendulaire.

Chapitre 2

Dans le deuxième chapitre, nous cherchons à déterminer la conductivité thermique effective à l'échelle macroscopique d'un problème à partir de ce qui se passe à l'échelle microscopique. Ainsi, nous considérons un milieu granulaire formé de trois composants : solide, l'air, et l'eau liquide, cette dernière étant répartie en régime pendulaire. Nous calculons la conductivité thermique eective (λ e ) de ce milieu en utilisant la loi de Fourier. Pour cette raison, nous considérons un modèle simple à l'échelle microscopique formé de deux grains sphériques identiques et nous nous concentrons sur le transfert de chaleur entre eux en présence d'air sec et de l'eau liquide qui est attachée aux grains solides sans gravité sous la forme d'un ménisque. Nous distinguons deux états diérents pour le ménisque : un état 'avec contact' lorsque le liquide fait un véritable pont entre les deux sphères et un état 'sans contact' lorsque le liquide est divisé en deux gouttes diérentes, séparées par une couche d'air mince. Tout d'abord, nous avons calculé la forme exacte de la surface du ménisque liquide en considérant qu'il est en équilibre mécanique, donc avec une courbure totale constante.

La forme géométrique du ménisque est obtenue en intégrant un système algébrique différentiel qui nous permet de calculer le volume d'eau liquide (ce volume servira plus tard pour calculer la fraction liquide du milieu à l'échelle macroscopique). Deuxièmement, nous calculons le ux de chaleur totale qui traverse notre domaine de calcul, en résolvant numériquement l'équation de la chaleur à l'équilibre avec les conditions aux limites adéquates. La discrétisation spatiale est eectuée en utilisant la méthode de Volumes Finis résultant dans un système linéaire qui est résolu par le solveur linéaire creux UMFPACK. Connaissant la température au niveau de chaque n÷ud de chaque cellule, on peut calculer le ux de chaleur. Une fois, le ux de chaleur obtenu, nous pouvons en déduire la conductivité thermique eective en utilisant la loi de Fourier.

Nous présentons la variation du ux thermique (proportionnelle à la conductivité thermique eective par la loi de Fourier) en fonction de la fraction volumique de liquide (c'est-à-dire l'humidité) où l'on observe un comportement d'hystérésis et un saut de la valeur du ux de chaleur en raison du passage entre les deux états possible du ménisque, avec ou sans contact.

Ce travail a été publié dans [CDM + 14] et accepté pour publication dans la revue internationale 'ASME J. of Heat Transfer' en 2015. Chapitre 3 Le chapitre trois est une continuation du travail eectué par M. Muhieddine dans sa thèse. Il a résolu l'équation de diusion de la chaleur lors de la fusion de la glace par la méthode LHA [START_REF] Prapainop | Simulation of ice formation by the nite volume method[END_REF] avec un schéma d'Euler explicite et par la méthode AHC [START_REF] Bonacina | Numerical solution of phase-change problems[END_REF][START_REF] Yao | An alternative formulation of the apparent heat capacity method for phase-change problems[END_REF] en utilisant un solveur ODE basé sur un schéma BDF et il a appliqué une technique de maillage adaptatif à la méthode LHA seulement. J'ai également résolu l'équation de diusion de la chaleur avec changement de phase (glace / eau) par la méthode LHA avec un schéma d'Euler explicite et la méthode AHC avec un schéma d'Euler explicite et la méthode BDF (Backward Dierentiation Formula) qui est basée sur un schéma implicite. Comme le schéma implicite BDF nécessite le calcul d'une matrice jacobienne, nous l'avons calculée et générée par le calcul formel (Maple), puis stockée au format creux, ce qui n'était pas le cas dans [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] où la jacobienne a été codé manuellement sous une forme approchée. L'équation de conduction de la chaleur utilisée pour modéliser le transfert de chaleur avec changement de phase est donnée par:

(ρC) e ∂T ∂t = div λ e gradT

(1)

Dans la méthode AHC, la chaleur latente de changement de phase est intégrée dans la capacité thermique C par l'introduction d'une distribution de Dirac, elle-même approchée par une gaussienne (voir gure 3.6). En utilisant la méthode AHC, le choix de l'intervalle de changement de température de phase ∆T aecte la précision de la solution de l'équation de la chaleur. Certains essais ont montré qu'il est raisonnable de considérer ∆T proportionnelle à h [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] (∆T = kh) où h est la taille des mailles et k est une constante choisie de manière à obtenir une bonne précision avec peu de uctuations pour la température en fonction du temps. La valeur de ∆T qui assure la meilleure précision de la solution avec quelques uctuations est appelée ∆T optimum . Pour mieux comprendre l'eet de ce paramètre sur l'histoire de la température, nous avons considéré plusieurs valeurs de ∆T : ∆T = ∆T optimum 4

, ∆T = ∆T optimum et ∆T = ∆T optimum × 4. de subdivisions autour de la zone de changement de phase (cellules mixtes dans la méthode LHA et dans l'intervalle de température de changement de phase pour la méthode AHC). Le cas de la fusion de la glace est choisi an de comparer la solution numérique avec la solution analytique. Pour chaque méthode utilisée, LHA (schémas explicites) et AHC (les schémas explicites et implicites), nous avons calculé l'erreur L 2 à la fois pour le prol de la température au temps nal et l'historique de la température à une position spécique et nous l'avons représentée comme variable en fonction de temps de calcul CPU. Nous avons conclu que LHA est la méthode la plus performante (Erreur / CPU), puis AHC utilisant un solveur ODE (implicite), enn AHC utilisant un schéma d'Euler explicite. L'ordre de convergence des LHA est de O(h 1.22 ), celui de AHC implicite est de O(h 0.91 ) et nalement celui de AHC explicite est de O(h 0.67 ).

Dans la deuxième partie de ce chapitre, nous présentons Zohour, un algorithme de maillage adaptatif 2D développé par É. Canot [Canc] Chapitre 4

Le chapitre quatre est consacré au problème inverse et est divisé en deux parties principales. Dans la première partie, nous avons amélioré le problème inverse 1D présenté dans [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] pour l'estimation des propriétés physiques de la matrice solide du sol (s) : conductivité (λ s ), la porosité (φ) et la capacité thermique volumique ((ρC ) s ) dans un milieu poreux initialement saturé et qui subit une évaporation. Le problème direct (diusion de la chaleur) résolu par la méthode AHC est utilisé pour créer les données synthétiques. Nous utilisons le critère des moindres carrés pour résoudre le problème inverse dans laquelle les coecients de sensibilité apparaissent et où nous essayons de minimiser la fonction coût qui est exprimée comme la norme de la diérence entre la température synthétique et les données numériques obtenues par notre modèle approché.

En utilisant le même nombre de mailles pour créer les données synthétique (à tous les n÷uds, mais au temps de calcul nal) et dans le problème inverse, M. Muhieddine a obtenu des valeurs proches mais pas exactes pour λ s et φ alors que (ρC) s a stagné à son estimation initiale. Nous avons amélioré le problème inverse en mettant à la même échelle tous les paramètres, en plus du calcul des données synthétique à des positions spéciques, mais à tous les pas de temps. Les valeurs de convergences obtenues sont exactement les valeurs requises (voir tableau 1). Dans le cas où le problème inverse ne réussit à converger à cause d'un mauvais choix des valeurs initiales, nous présentons une nouvelle stratégie pour le résoudre en deux étapes. Tout d'abord, nous calculons le problème inverse en utilisant un intervalle de changement de phase de la température (∆T ) supérieur à (∆T ) optimum utilisé dans le problème direct de manière à obtenir un historique de température lisse avec un minimum de uctuations. Ensuite nous exécutons de nouveau le problème inverse (∆T ) optimum en utilisant les valeurs de convergence obtenues à l'étape 1 comme estimation initiale, nous appelons cette technique "enchaînement de problèmes inverses".

Dans la deuxième partie de ce chapitre, nous présentons le même problème inverse, mais dans un système de coordonnées 3D-axisymétrique. Le système d'équations comme celui dans le cas 1D est composé de l'équation de la chaleur avec trois équations de sensibilité résultant de la dérivée de l'équation de diusion de la chaleur par rapport aux trois paramètres inconnus. Le système de quatre équations couplées non linéaire forme un système d'équations aux dérivées partielles. La discrétisation est eectuée en utilisant la méthode des lignes où les discrétisations spatiales et temporelles sont faites séparément. La discrétisation spatiale est eectuée en utilisant la méthode de Volumes Finis, ce qui permet d'obtenir un système semi-discret d'équations qui est résolu par un solveur ODE automatique basé sur un schéma BDF dont la Jacobienne est générée par Maple et stockée en format creux.

Une analyse de sensibilité de paramètres a été eectuée à la n de cette seconde partie et montre que la porosité est le paramètre le plus identiable, puis vient la conductivité et enn la capacité thermique volumique.

Ce travail a été publié dans [START_REF] Mansour | Identication of thermophysical properties of the soil in 3d-axisymmetric coordinate system using inverse problem[END_REF] et et accepté pour publication dans la revue internationale 'ASME J. of Heat Transfer' en 2015.

Chapitre 5

Le Chapitre 5 est essentiellement l'application d'un problème inverse 3D-axisymétrique pour des données expérimentales. Les premières données expérimentales ont été réalisées par Szubert [START_REF] Szubert | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] et Bergonzi [START_REF] Bergonzi | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] et correspondent au cas d'un milieu saturé en eau, tandis que les secondes ont été réalisées par Cordero [START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF] correspondent au cas d'un milieu sec. Le problème inverse n'a pas réussi à converger dans le premier cas; cela est dû au fait que le modèle physique utilisé n'est pas assez réaliste et a besoin d'être amélioré :

aucun écoulement de vapeur dans le milieu poreux (mais il pourrait être ajouté et testé).

les forces de gravité et capillaires sont négligées (la gravité ne peut pas être prise en compte à cause d'une impossibilité dans la méthode AHC).

le rayonnement est négligé bien qu'un chauage fort est utilisé (300 à 600 degrés Celsius). On sait que la silice est semi-transparente à certaines longueurs d'onde.

Par conséquent, bien que le transfert conductif soit dominant, il y a aussi un transfert radiatif.

un modèle d'ébullition (température > 100 degrés Celsius) n'est pas utilisé bien qu'il y ait certainement une zone d'ébullition constituée de bulles de vapeur (processus non stationnaire); un capteur à l'intérieur de cette zone resterait à 100 degrés Celsius.

Par conséquent, un modèle plus complexe est nécessaire; il sera discuté dans le dernier chapitre.

Le problème inverse a donné de bons résultats pour le deuxième ensemble de données expérimentales (cas sec), cependant il faut préciser que seule la diusivité thermique α est identiable.

Par ailleurs, d'autres méthodes connues (comme celle de Laloy & Massard [START_REF] Laloy | Nouvelle méthod thermique d'étude des foyers préhistoriques[END_REF] et celle de la transformée de Laplace [START_REF] Carslaw | Conduction of heat in solids[END_REF]) ont été utilisées pour les expériences en sec : la diusivité estimée est α = 2 -3 × 10 -7 m/s 2 pour le sable de Fontainebleau.

En utilisant notre problème inverse, où le problème des moindres carrés est résolu par la méthode LMA nous avons obtenu de mauvais résultats α = 3 × 10 -5 m/s 2 : nous enregistrons une erreur d'un facteur 100 par rapport aux valeurs de diusivité pour le sable de Fontainebleau dans la littérature et ceux obtenus par transformation de Laplace.

Nous avons supposé que l'erreur obtenue provient, principalement, de la méthode pour mesurer les positions des capteurs. Nous avons donc reformulé le problème inverse en recherchant α en même temps que les positions des capteurs, avec les contraintes additionnelles (linéaires) suivantes : non-négativité des diérents paramètres recherchés et somme des erreurs de la position des capteurs nulle. Cette dernière contrainte suppose toutefois que les erreurs sur la mesure de la position des capteurs obéit à la loi normale, ce qui est sans doute le cas pour des expériences utilisant un grand nombre de capteurs.

Avec cette nouvelle reformulation, les résultats obtenus sont plus réalistes car plus proches de ceux trouvés dans la littérature α = 1.73 × 10 -7 m/s 2 pour la première expérience et α = 2.42 × 10 -7 m/s 2 pour la deuxième.

Ce travail a été présenté à la 7ème Conférence Internationale 'Porous Media & Annual Meeting (Interpore), May 18 -21, 2015, Padova, Italy'.

Chapitre 6

Dans ce chapitre qui pose les bases de futurs modèles physiques, nous prenons en compte la présence d'air qui vient se mélanger par diusion à la vapeur d'eau. Deux cas totalement diérents sont présentés, dans lesquels le milieu est alors formé de trois constituants : solide (matrice poreuse), liquide (eau) et de gaz (vapeur d'eau + air).

Dans une première partie, nous présentons un modèle à phase séparées où nous traitons le transfert de chaleur et de masse dans un domaine 1D en présence d'air le but visé est de pouvoir traiter l'alternance de phases de chauage et de refroidissement (avec changement de phase). Le système d'équations comporte : l'équation de diusion Mes contributions les plus importantes et originales sont les suivantes : Amélioration d'une méthode numérique calculant la conduction de chaleur en milieu poreux saturé avec changement de phase :

-Étude de l'inuence de l'intervalle de température ∆T (paramètre de régularisation) sur la précision de la méthode AHC.

Pour traiter l'alternance des phases de chauage/refroidissement en milieu poreux saturé avec changement de phase :

-Établissement d'un nouveau modèle physique à zones séparées, qui prend en compte la diusion de l'air dans la vapeur d'eau ;

-Présentation des premiers résultats numériques, avec une méthode numérique (qui ne pouvait pas être AHC).

Maillage adaptatif :

-Utilisation systématique de maillage adaptatif (1D) pour le problème de conduction avec le changement de phase en utilisant trois méthodes diérentes (LHA, AHC 'explicite', AHC 'implicite' en utilisant solveur ODE automatique) an de comparer leurs performances respectives (plus précisément, le rapport performance/coût) ;

-Introduction d'un maillage adaptatif 2D pour résoudre le transfert de chaleur par conduction.

Problème inverse :

-Introduction de la mise à l'échelle des paramètres recherchés et Étude de l'eet de ∆T sur la convergence de paramètres (enchaînement des problèmes inverses : résoudre le problème inverse en deux étapes) ; -Application du problème inverse à des données expérimentales où nous avons ajouté la position des capteurs comme inconnues, avec des contraintes additionnelles pour fermer le système, de type statistique (si les capteurs sont en grand nombre, on suppose que l'erreur de leur position suit une loi Normale) et de type géométrique (déplacement possible des capteurs uniquement dans la direction de plus grande sensibilité).

Calcul de la conductivité thermique eective (λ e ) en présence de petite quantité d'eau liquide entre deux sphères ; il s'agit d'un modèle simple pour un domaine à trois composants (solide, liquide et gaz) :

-Calcul de la forme exacte du ménisque liquide par intégration d'un système DAE (la plupart des études existantes utilisent un modèle toroïdal) ; -Le comportement d'hystérésis du ux de chaleur en fonction de la fraction liquide, c'est-à-dire de l'humidité.

3.1

The Latent Heat Accumulation Method (LHA) . . . . . . . . . . . . . .

3.2

The Apparent Heat Capacity Method (AHC) . . . . . . . . . . . . . . . The work in this thesis takes place in an interdisciplinary project which started more than ten years ago. The ARPHYMAT project [Cana] is devoted to the study of prehistoric res, by grouping the three complementary domains, ARcheology, PHYsics and MAThematics. It aims to estimate the duration of ancient human groups' res, in order to understand the reason behind using the re (cooking, heating, potteries,...)

Conclusion Bibliography

and also the social behavior of these groups.

All the tasks (T ) that describe this project are shown in the diagram on the next page.

Some experiments have been done (T 2 ), under real situations (burning wood pieces) and in laboratory (replacing re by electric heating plate). The role of these experiments is mainly the determination of the way heat transfers in the soil in order to understand the drying phenomenon and its consequences on real soil. Moreover, these

Introduction

Try to estimate:

• Duration of fire.

• Thermophysical properties of the porous medium.

ARPHYMAT interdisciplinary project. Study of prehistorical fires http://arphymat.univ-rennes1.fr/ Characteristics:

• Strong heating of the soil (porous medium).

• Presence of water. Physical models:

• Current models: Theoretical calculation of 𝜆 𝑒 :


• Exact calculation of the shape of liquid meniscus between two spherical grains.

•

Homogenization

• Hysteresis behavior of the heat flux. flow experiments allow us to estimate some thermophysical properties of the soil using some appropriate methods (T 3 and T 5 (B)). Finally, they help us in quantifying the inuence of the quantity of water in soil (recall that the phase change absorbs a very important part of the thermal energy).

Besides the experiments (T 2 ), numerical simulations (T 5 ) are needed in order to better understand the inuence of some parameters. These simulations are based on physical models (T 4 ).

The numerical codes used in T 5 are laboratory codes, done by É. Canot and M.

Muhieddine (2009). The reason behind this choice is the capability to have easy access and thus modifying the numerical methods and introducing new physical models that are sometimes not conventional (we will see later on that the eective conductivity, in presence of small quantity of liquid water, exhibits a hysteresis behavior).

The physical models in T 4 present dierent degrees of diculty. The M 1 model has been used by É. Canot and M. Muhieddine for several years but it suers from several drawbacks. For this reason, the M 2 model has been developed to treat alternating heating/cooling stages with phase change but without using the AHC method. A third, more realistic, model (M 3 ) is proposed but it is not yet implemented.

Another part of the work, T 6 , concerns the determination of the eective thermal conductivity of a wet porous medium. Although many results concerning the dry case can be found in literature, the case with three components (solid, liquid water and gas)

does not exist. For this reason, we made a rst simple model which shows a hysteresis behavior. T 6 permits the calculation of the eective thermal conductivity which in its turn is needed for the future model M 3 .

In this thesis, we are interested in the problem of heat transfer in a porous medium (saturated and unsaturated). In the rst part, we compare dierent methods that solves the phase change problem with or without the use of an adaptive mesh technique in a saturated porous medium. In the second part, we present a model to deal with phase change in an unsaturated porous medium paying a special attention to the determination of the eective thermal conductivity λ e as function of liquid content. The third part concerns the determination of the thermophysical properties of the soil by inverse problem.

Outline

In the rst chapter, we present a general description of a granular medium where we dene the physical properties related to our study: thermal conductivity, density, capacity, surface tension (capillary forces), the contact angle of a liquid-gas interface in contact with a solid and the water distribution layouts in granular medium especially the pendular regime.

Chapter two explores the possibility to nd the eective thermal conductivity of a granular medium in pendular regime formed of three components: solid, liquid water and dry air where we aim to calculate the eective thermal conductivity (λ e ) of this medium using Fourier law. For this reason, we consider a simple model at a microscale level formed of two identical grains and we focus on heat transfer between them in presence of dry air and liquid water which is attached to the solid grains without gravity in the form of a meniscus. We rst nd the exact shape of the liquid meniscus by integrating a dierential algebraic system which allows us to calculate the volume of liquid water later on. Then, we calculate the total heat ux which ows through our computational domain by solving the steady state heat equation with adequate boundary conditions. Once the heat ux is obtained, we can deduce the eective thermal conductivity using Fourier law. We assume that the medium (at macroscopic scale) is periodic with a known microstructure and thus an approximative homogenization technique can be applied (upscaling method) [START_REF] Auriault | Homogenization of Coupled Phenomena in Heterogeneous Media[END_REF].

This work was published in [CDM + 14] and accepted for publication as a research paper in ASME J. of Heat Transfer.

Chapter three is a continuation of the work performed by M. Muhieddine in his thesis. He solved the heat diusion equation during fusion of ice by LHA [START_REF] Prapainop | Simulation of ice formation by the nite volume method[END_REF] with Euler explicit scheme and by AHC method [BC73, YC93a] using an ODE solver based on a BDF scheme where he applied an adaptive mesh technique to LHA method only.

In the rst part of this chapter, we also solve the heat diusion equation with phase change (ice/water) by LHA with explicit Euler scheme and AHC with Explicit Euler scheme and backward dierentiation formula (BDF) implicit scheme.

We give a special attention to the choice of the phase change temperature interval ∆T which aects the accuracy of the solution of the phase change problem using AHC method. An adaptive mesh renement is applied to both LHA and AHC methods to improve their performance.

In the second part of this chapter, we present 'Zohour', a node-based original adaptive 2D mesh algorithm developed by É. Canot [Canc]. This adaptive mesh is applied to the AHC method to solve the fusion problem in a 2D coordinate system. The advantage of using this algorithm to any other existing adaptive mesh algorithm is the fact that it is a home-made algorithm and any needed modication can be easily applied.

Chapter four is devoted to the inverse problem and is divided into two main parts. In the rst part, we improved the 1D inverse problem presented in [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] for the estimation of physical properties of the soil solid matrix: conductivity, porosity and volumetric heat capacity in a porous medium during evaporation. This improvement was done by adapting a parameter's scaling technique and adding a special technique in solving the inverse problem using dierent values of ∆T over two or more steps to insure convergence (we call this technique: chaining of inverse problems). The code validation stage was based on the comparison between numerical results and synthetic data.

In the second part of this chapter, we present the same inverse problem but in a 3D axisymmetric coordinate system. In both parts, the system of equations is composed of the heat equation together with three sensitivity equations resulting from dierentiating the heat diusion equation with respect to the three unknown parameters. The system of four coupled equations forms a non-linear system of partial dierential equations. We discretize the equations using the method of lines where spatial and time discretizations are done separately. The spatial discretization is done using FVM and thereafter we end up having a semi-discrete system of ODEs which is solved by an automatic ODE solver based on a BDF scheme.

The work of this chapter is already published in [START_REF] Mansour | Identication of thermophysical properties of the soil in 3d-axisymmetric coordinate system using inverse problem[END_REF] and accepted for publication as a research paper in ASME J. of Heat Transfer.

Chapter ve is basically the application of a 3D axisymmetric inverse problem on experimental data. The rst kind of used experimental data were those obtained for a water saturated case performed by Szubert [START_REF] Szubert | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] and Bergonzi [START_REF] Bergonzi | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] in a joint work, while the second were for a dry case performed by Cordero [START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF]. The inverse problem failed to converge in the rst case due to the fact that the physical model used is not close enough to reality and need to be improved. Whereas, the inverse problem provided good results for the second set of experimental data (dry case). The inverse problem in the second case is restricted to the identication of the thermal diusivity α together with the sensors' positions, adding appropriate additional constraints. Note that archaeologists are interested in such a problem because the ancient hearths are made up of a sequence of heating and cooling stages.

This work was presented in the 7th

In the second part, we present, without attacking, a new physical model that can be used to deal with the phase change phenomenon in a wet unsaturated porous medium with low humidity (small quantity of liquid water: pendular regime). This model is welladapted to the ARPHYMAT experiments and will be the subject of a future work.

Chapter 1

General Physical Description

In this section we present a general description of the porous medium and the physical parameters involved in our studies.

Porous medium

A porous medium is a body composed of a persistent solid part, called solid matrix, and the remaining void space (or pore space) that can be lled with one or more uids ( e.g. water, oil and gas).

Typical examples of a porous medium are soil, sand, cemented sandstone, karstic limestone, foam rubber, lungs and kidneys.

We can dierentiate between:

A single phase system where the void space of the porous medium is lled by a single uid (e.g. water) or by several uids completely miscible with each other (e.g. fresh water and salty water, or air and water vapor).

A multiphase system where the void space is lled by two or more uids that are immiscible with each other i.e. they maintain a distinct mobile boundary between them (e.g. water and oil or water and air ).

In all our study, the solid matrix is considered as non-deformable.

Continuum Approach

The important feature in modeling porous media ow is the consideration of dierent length scales.

Macroscopic scale → Microscopic scale → Molecular scale.

The macroscopic scale is of order of 10 meters, if we zoom into a scale of about 10 -3 m we arrive at the microscopic scale where individual sand grains and pore channels are visible see gure 1.1. In order to derive a mathematical model on the macroscopic level another continuum is considered. Each point in the continuum on the macroscopic level is assigned average values over elementary volumes of quantities on the microscopic level. This process leads to macroscopic equations that do not need an exact description of the microscopic conguration [START_REF] Vafai | Handbook of Porous Media[END_REF]. Anisotropy typically comes from the conguration at microscopic scale.

Some Physical Parameters Porosity

Porosity is the ratio of void volume (pores volume) to the total volume:

φ = P ores V olume T otal V olume (1.1)
The following table represents typical values of porosities for dierent substances [START_REF] Kaviany | Principles of Heat Transfer in Porous Media[END_REF]:

Saturation

Saturation is the measure of the uid volume present in the pore volume of a porous medium. By denition, the saturation of a uid is the ratio of the uid volume to the pore volume.

In our study, α represents the degree of saturation of liquids:

α = V olume of liquid V olume of pores (1.2)
The porous medium, in this work, is considered saturated when the liquid fraction in the pores α is equal to 1.

Permeability

Permeability is the measure of the ability of a porous medium to permit uids to ow through its pores or voids. It depends only on the geometry of the solid matrix, in particular, pore size. The greater the void ratio, the higher the value of the permeability.

Usually, the permeability K is determined by experimental measurements through Darcy's law given by equation:

v f = - K µ f ∇P (1.3)
v f is the ltration velocity of Darcy for the uid phase, µ f is the viscosity of the uid and P is the pressure of the uid.

The permeability of some porous mediums are shown in the following table :  chapter 1 1.3 × 10 -14 -5.1 × 10 -14

Gravel 10 -7 -10 -9

Specic heat capacity

The specic heat capacity (C ) of a substance is the amount of heat per unit mass required to change its temperature by one degree. The units of C is J/kgK.

Density

The density (ρ), or more precisely, the volumetric mass density, of a substance is its mass per unit volume (kg/m 3 ).

Volumetric heat capacity

Volumetric heat capacity (ρC ) is the quantity of heat required to change the temperature of a unit volume by one degree (SI unit: J/m 3 K).

Thermal conductivity

Thermal conductivity (λ) is the property of a material to conduct heat. Heat transfer occurs at a lower rate across materials of low thermal conductivity than across materials of high thermal conductivity. It is measured in Watts per meter Kelvin (W/mK ).

Many factors aect the soil thermal conductivity (λ) like soil water content, mineral composition, temperature, porosity, ... Fourier's Law, the basic rate equation of the conduction process, given by the equation below, provides the denition of the thermal conductivity and forms the basis of many methods of determining its value.

q = -λ∇T

(1.4)
q is the local heat ux density in W/m 2 , λ is the thermal conductivity in W/mK and ∇T is the temperature gradient in K/m.

Surface tension

The surface tension γ is the elastic tendency of liquids which makes them acquire the least surface area possible. At a typical liquid-air interface, surface tension results from the greater attraction of water molecules to each other (due to cohesion) than to the molecules in the air (due to adhesion). Because of the relatively high attraction of water molecules for each other, water has a high surface tension (72.8 × 10 -3 N/m at 20 • C) compared to that of most other liquids. Surface tension is an important factor in the phenomenon of capillarity.

Contact angle

Consider thermodynamic equilibrium between three phases: the liquid phase (l), the solid phase (s) and the gas/vapor phase (g). If the solid-gas interfacial energy is denoted by γ sg , the solid liquid interfacial energy by γ sl and the liquid-gas by γ lg (i.e. the surface tension) then the equilibrium contact angle θ is determined by Young's equation [START_REF] Quéré | Rough ideas on wetting[END_REF] (see gure 1.2):

cosθ = γ sg -γ sl γ lg (1.5)
chapter 1

Figure 1.2: The contact angle between the liquid interface and the solid

Porous medium properties

The transient heat conduction equation in a porous medium is given by:

∂ ∂t [(ρC) e T ] = div (λ e ∇T ) (1.6)
where the subscript e represents the eective property.

To model the heat transfer in porous medium, we need to calculate the eective physical properties:

Eective volumetric heat capacity

The eective volumetric heat capacity (ρC) e which is proportional to the volume fraction and is given by:

(ρC) e = φ(ρC) f + (1 -φ)(ρC) s (1.7)
where ρ is the density, C is the specic heat capacity and φ is the porosity. The subscripts e, f and s indicates the eective parameters of the medium, the properties of the uid and porous matrix properties respectively.

The simple relation presented in equation 1.7 is due to the fact that the product (ρC) is an extensive variable with respect to the volume which is not the case for the specic heat capacity.

Eective thermal conductivity

The models used to evaluate the eective thermal conductivity λ e are based on the identication of the porous medium as a regular simple geometric structure. The conductivity is then calculated by rigorous theoretical resolution, or by numerical calculation using simplifying assumptions [GLB + 13]. In the latter case, we get explicit formulations for λ e as function of the conductivities of the present phases and porosity (see table 1.5).

Table 1.5: Models used to calculate the eective thermal conductivity in a granular porous medium [START_REF] Bories | Transferts de chaleur dans les milieux poreux[END_REF].

Model Expression Remarks

Series (harmonic mean)

1 λ e = φ λ f + 1 -φ λ s Maxwell (lower law) λ e λ f = 2φ + (3 -φ) λs λ f 3 -φ + φ λs λ f
A medium consisting of spheres dispersed in a continuous uid phase where no inuence between particles.

Geometric It should be emphasized that the use of the presented models must be carefully done where in all circumstances one must ensure that they are applicable to the situations and mediums considered, not to mention that all the models previously presented take chapter 1 into account only two components and there are three in our case (solid, liquid and gas).

λ e = λ φ f λ 1-φ s Kunii & Smith λ e = λ f φ + (1-φ)β
λ e λ f = 2 λs λ f 2 (1 -φ) + (1 + 2φ) λs λ f (2 + φ) λs λ f + 1 -φ Medium formed of
The geometric structure of the medium plays an important role in the calculation of λ e which cannot be uniquely determined from the knowledge of the porosity and the thermal conductivities of each phase.

As a consequence, the use of experimental determination of λ e is generally essential for a truly accurate value. In chapter 2, we present a procedure to nd the eective thermal conductivity of granular medium formed of three constituents (solid matrix, liquid water and air) at a macroscopic scale out of a problem at microscopic scale.

1.4 Water content regimes in a wet granular medium Granular materials are collections of macroscopic particles, like glass beads or sand, which are visible by the naked eye. In a granular medium, water molecules are preferentially attracted to the surface of the grains, they can be adsorbed from the vapor present in the pore air. Thus, small amounts of water are always present in the form of thin lms covering the surface of the solid skeleton. As the amount of water in a porous medium increases, it is attracted to the water adsorbed at the grain's surface by the cohesion forces, but on the other hand tends to minimize the area of the water-air interface. This form of water is known as the capillary water, and occurs initially in acute corners of the pores forming pendular rings around the contact points (see gure 1.4 a.). This water-air conguration is known as the pendular regime [ILHE00, NCJ58].

Pendular water is less tightly bound by the solid phase, but it occurs in isolated regions, and does not form continuous ow paths, so it can be considered as macroscopically immobile [START_REF] Pinder | Subsurface Hydrology[END_REF].

If more water is added to the system, the regions occupied by pendular water coalesce and continuous thicker lms are formed along the pore walls. At this stage, known as the funicular stage, the ow of liquid water is possible (see gure 1.4 b.).

The saturated case (see gure 1.4 c.) is very dicult to be obtained unless special care is taken (industrial process). • No liquid water.

• Full of air and/or water vapor.

Dry

Cohesion between grains is negligible.

• Small quantity of liquid water.

• Air and/or water vapor.

Pendular

Wet

• Liquid bridges are formed at the contact points of the grains.

• Cohesive forces act through the liquid bridges.

• Liquid water doesn't move due to capillary forces.

• Middle quantity of liquid water.

• Air and/or water vapor.

Funicular

• Liquid bridges around the contact points and liquid-lled pores coexist.

• Both give rise to cohesion between particles.

• Liquid water can move due to pressure gradient, gravity and capillary forces.

• Full of liquid water.

• No gas is present.

Saturated

• No cohesive interaction appears between particles.

• Liquid water move due to pressure gradient only. The eective thermal conductivity is one of the most important parameters characterizing the energy transport properties of porous media and has been studied extensively by using both theoretical and experimental approaches. As it is well known, the eective thermal conductivity of porous media depends not only on the thermal property and the volume fraction of each constitute component, but on the micro-structure (i. e. the spatial distribution of all components) of the media as well. On the other hand, it is known that the presence of water aects a lot the thermal properties because the thermal conductivity of water is about 25 times that of dry air. This has been mentioned by Chen [Che08] who emphasized that a small increment in water contents around the contact points will lead to rapid increment in the eective thermal conductivity. Some papers [ULMF06, YS08] deal with the prediction of thermal conductivity in terrestrial soil media in the presence of liquid water as a component but without specifying its spatial localization. The work of Mitarai and Nakanishi [START_REF] Mitarai | Simple model for wet granular materials with liquid clusters[END_REF] concerns the pendular regime (see gure 2.1) in which the liquid content is small and where the water menisci form bridges at each contact point (between the grains) thus inducing two-body cohesive force so the study was restricted to the mechanical properties and not the thermal ones. In this chapter, we explain a simple procedure to nd the eective thermal conductivity at a macroscopic scale (having a periodic pattern) out of a problem at microscopic scale. This called upscaling or homogenization (we will come back to this procedure at the end of this chapter). Assuming a periodic pattern is classically used by some authors (Luikov et Slurry (Droplet) regimes. We are concerned with the pendular regime where we have a small liquid content and where the liquid bridges are formed at the contact points of grains due to the surface tension responsible for tendency of liquids to minimize their surface area, giving rise to a membrane-like surface [START_REF] Tro | Chemistry: A Molecular Approach[END_REF]. We dierentiate between two states, contacting state when the liquid water makes a true bridge between the two spheres and non-contacting state when the liquid split in two dierent drops separated by a thin air layer.

Figure 2.1 represents a real granular medium in pendular regime. This three component medium has no theoretical model to calculate the eective thermal conductivity and it is dicult to be discovered by experimental results. However, a numerical value of this eective property can be found in two stages, that is to say rst nding the equilibrium shape of the liquid meniscus for a given liquid volume, second solving the steady-state heat equation over a small domain (supposed to be representative enough of the whole domain) and deduce the eective thermal conductivity, via the Fourier law, by using the thermal ux through the domain.

q = -λ e .∇T

(2.1)

where q is the local heat ux density in W/m 2 , λ e is the eective thermal conductivity in W/mK and ∇T is the temperature gradient in K/m. Fourier law shows that the time rate of heat transfer through a material is proportional to the negative gradient in temperature.

In this chapter, we focus on the heat transfer between two spherical solid grains of same radius in presence of dry air and where liquid water is attached to the solid grains, without gravity. Our model takes into account the role of the contact angle between the liquid-gas interface in addition to the surface tension (energy required to increase the surface area of a liquid by a unit area) which enable the calculation of the exact shape of the liquid meniscus. A simplication of the grains' assembly is represented in gure 2.2. We suppose that the whole medium presents a periodic pattern. The present model implies many assumptions; there is no doubt that some important characteristics of a real granular medium (especially the randomness of both shape and position of the grains, and the roughness of their surface) should aect the results obtained in this work. On the other hand, we can state that gravity doesn't aect a lot the shape of the liquid meniscus (the Bond number, based on a reference length chosen as the height of the meniscus, is of the order of one thousandth for a grain of radius 1 mm). Nevertheless, we think that our simple model has the worthiness to reveal a very interesting behavior, that is the hysteresis during the change in the liquid volume.

Lastly, although the true temperature eld is 3D, we suppose that it is well described by the local axisymmetric eld around one symmetry line joining two spheres.

Micro-scale Model

We want to solve the steady-state heat transfer in an elementary cell containing only two spherical solid grains, a few quantity of liquid water, and dry air. In a real situ-chapter 2 ation, the liquid water should be in thermodynamic equilibrium with its vapor, so the surrounding gas should be a binary mixture of dry air and water vapor. In the present study, the water vapor is neglected and we assume that there is no adsorbed water on the grain surface. The computational domain (see gure 2.3) can be found using the symmetry planes of the problem knowing the main direction of the heat ux -for example, the bottom symmetry plane comes from the fact that the gravity is neglected. Top and bottom sides are kept at constant (but dierent) temperature whereas the vertical ones are supposed to be isolated (this is related to the model used by Yun [START_REF] Yun | Fundamental study of thermal conduction in dry soils[END_REF] in their study of thermal conduction in dry soils). As stated before, the liquid meniscus is centered along the axis between the two spheres leading to an axi-symmetric problem, so we can reduce the geometrical independent variables to the (r, z) cylindrical coordinates. Each component (solid, liquid water and dry air) has a constant thermal conductivity and the thermal contact between the components is supposed to be perfect.

Because it is in mechanical equilibrium (actually, in a quasi-static state because in a real situation a change of volume can be expected due to many reasons), the shape of the liquid meniscus must have both a constant total curvature (it is a minimum surface) and a prescribed contact angle with the solid. The percentage of liquid water, i. e. the ratio of its volume to that of the available space between the spheres is called humidity.

The input parameters are: R 0 : the radius of the two solid spheres : half the gap between the spheres V : the volume of liquid water γ : the surface tension between liquid water and air θ : the contact angle between the liquid interface .

and the solid In this study, the radius R 0 is taken as the reference length. The output variable of the problem is the (vertical) ux q across the domain which is related to the eective thermal conductivity of our composite medium (by Fourier law). Note that in a complete model, the liquid volume V should be related to the gas pressure (in a real situation, the total mass of the water is xed and then the repartition between gas and liquid depends on the local pressure) but this is out of the scope of this study.

Surface tension γ and contact angle θ are usually temperature dependent and this dependence (which can be found elsewhere in literature) is not so negligible. However, at the scale of our elementary cell, we can state that the temperature dierence is small so that these two parameters can be treated as constants.

Equilibrium shape of the liquid meniscus

In physics, the Young-Laplace equation relates the pressure dierence to the shape of the interface [START_REF] Guyon | Physical Hydrodynamics[END_REF]:

∆P = γ κ (2.2)
where ∆P is the pressure dierence across the uid interface, γ is the surface tension and κ is the total curvature. At micro-scale, ∆P has a constant value and this could be easily demonstrated using the hydrostatic equation (2.3).

∆P = -ρg∆z

(2.3)
where ρ is the density of uid, g is the gravitational constant and ∆z is the change in elevation.

Moreover, γ is supposed to remain constant so according to equation (2.2), κ is constant as well.

κ = κ 1 + κ 2 = constant (2.4)
where κ is the total curvature of the liquid meniscus curve (C) and κ 1 and κ 2 are respectively its plane and axi-symmetric curvatures. Let {r(s), z(s)} be any parameterization of the curve (C) in the axi-symmetric cylindrical coordinates then the total curvature is given by:

κ = r z -z r (r 2 + z 2 ) 3 2 + z r (2.5)
chapter 2

where z = dz/ds, r = dr/ds, z = d 2 z/ds 2 and r = d 2 r/ds 2 . If we choose s such that it corresponds to the arc length parameterization of curve (C) then:

r 2 + z 2 = 1 (2.6)
and thus equation (2.5) becomes:

r z -z r + z r = κ (2.7)
The geometrical shape of the liquid meniscus is obtained by integrating the dierential algebraic system formed of equations (2.6) and (2.7) (indeed, after transforming the system in a rst order type, the former equation appears under an algebraic form).

Regarding the initial conditions for this integration, the following explanation should be emphasized: we already mentioned in section 2.1 that both the volume V and the contact angle θ were chosen as input parameters; actually, it is dicult to integrate the meniscus shape (equations (2.6-2.7)) under these two constraints. We have preferred instead to dene a point in the (r, z) plane (point A in gure 2.4) which is the beginning of the integration; the integration is then stopped when the meniscus curve crosses the solid boundary, (point B in gure 2.4). Please note that the point A is located on the r-axis only for the contacting case; for the non-contacting case, the point A must be located on the z-axis. The constraint concerning the prescribed contact angle at B, on the contrary, is kept and the integration process is repeated by trial and error until obtaining the good contact angle.

To summarize, the initial conditions for the integration of (2.6-2.7) are as follows:

contacting case:

   r = r 0 z = 0 z = 0 non-contacting case:    r = 0 z = z 0 r = 0
It is worth noting that since the contact angle is prescribed, the total curvature κ and the meniscus volume V are conjugates in the sense where the knowledge of one of them leads to that of the other. If we choose a value of V , the dierential system (2.6-2.7) will be dicult to be integrated, because of the presence of two constraints: the contact angle θ and the liquid volume V . To simplify the integration, we prefer to choose a value of the total curvature κ and deduce the volume V (the computation is detailed afterwards). Anyway, the nal curves of the thermal uxes will be plotted with respect to the volume V , as a control parameter. (2.8) Equation (2.8) is actually added to the system of equations (2.6-2.7) and therefore V b is a result of the integration process. To nd the volume of liquid water V we must subtract the volume of the spherical cap V c (in brown) which is given by: V c = πh 2 3 (3R 0 -h)

The volume of the liquid water V

where h = z Band R 0 = 1 (radius of the grain). Finally, V = V b -V c. A is an arbitrary point dened as the beginning of the integration process.

2.2.2

The contact angle θ between liquid and solid interface and the intersection point between the liquid water and the solid grain, we can easily calculate θ. Let ψ be the angle formed by the tangent to the liquid meniscus at B then ψ = arctan(z B /r B ) (the prime denotes the derivative w.r.t. the curvilinear abscissa s) and assume that ϕ is the angle formed by the tangent to the solid grain at B then ϕ = arcsin(r B /R 0 ). Simply θ = ψ -ϕ. The solution of our system of equations (2.6-2.7) is obtained numerically by a shooting method over the starting position of the integration (r 0 in the contacting state or z 0 in the non-contacting state); this starting position is varied until the prescribed contact angle θ is found. The family of the curves obtained when we vary the total curvature κ, is drawn in gures 2.6 and 2.7 for the two states. chapter 2

In the non-contacting state (no gravity), the shapes of the liquid menisci correspond exactly to spherical caps. We could calculate these shapes mathematically without the need to integrate a DAE system.

Numerical Computation of the heat ux

Solving our problem numerically appears unavoidable: some authors (resp. [START_REF] Batchelor | Thermal or electrical conduction through a granular material[END_REF][START_REF] Garrett | Compressive pressure dependent anisotropic eective thermal conductivity of granular beds[END_REF]) have used analytical solutions but they are respectively restricted to an asymptotic behavior (so, adapted to a local geometrical zone) or to too crude approximations.

Consider the steady state heat equation in 3D axisymmetric cylindrical coordinates:

1 r ∂ ∂r λ(r, z) r ∂T ∂r + ∂ ∂z λ(r, z) ∂T ∂z = 0
(2.9)

with homogeneous Dirichlet boundary conditions on top and bottom sides and zero

Neumann boundary conditions (i. e. isolated conditions) on the external vertical sides.

This steady state equation is solved using a Finite Volume scheme applied on a regular structure mesh of rectangular cells. The detailed discretization over an interior cell with node (i,j) is explained below:

i+ 1 2 i-1 2 j+ 1 2 j-1 2 1 r ∂ ∂r λ(r, z) r ∂T ∂r r dr dz = λ i,j + λ i+1,j 2
r i + r i+1 2 T i+1,j -T i,j ∆r - λ i,j + λ i-1,j 2
r i + r i-1 2 T i,j -T i-1,j ∆r ∆z and i+ 1 2 i-1 2 j+ 1 2 j-1 2 ∂ ∂z λ(r, z) ∂T ∂z r dr dz = λ i,j + λ i,j+1 2
T i,j+1 -T i,j ∆z - λ i,j + λ i,j-1 2 T i,j -T i,j-1 ∆z r i ∆r
The whole equation becomes:

(λ i,j + λ i-1,j ) (r i + r i-1 ) r i ∆ 2 r T i-1,j + (λ i,j + λ i+1,j ) (r i + r i+1 ) r i ∆ 2 r T i+1,j - (λ i,j + λ i+1,j ) (r i + r i+1 ) r i ∆ 2 r - (λ i,j + λ i-1,j ) (r i + r i-1 ) r i ∆ 2 r + 2 (λ i,j + λ i,j+1 ) ∆ 2 z + 2 (λ i,j + λ i,j-1 ) ∆ 2 z T i,j + 2 (λ i,j + λ i,j-1 ) ∆ 2 z T i,j-1 + 2 (λ i,j + λ i,j+1 ) ∆ 2 z T i,j+1
(2.10)

where ∆r and ∆z represents the distance between two consecutive nodes in the r and z directions respectively. λ i,j is the value of the conductivity at the node of the specied cell a geometric test is done to specify the nature of the cell (solid, water or air) to specify the value of the conductivity at the requested cell.

After writing the discretization on the other boundary cells, the equations will be written as a linear system and the matrix of this linear system is stored in a sparse way and the UMFPack linear sparse solver is used. Some tries have shown that a 500 by 500 mesh is required to obtain a good accuracy (see gure 2.9). Each numerical computation takes about few seconds on a laptop (Intel Core i7 @ 2.7 GHz).

Now we have obtained by numerical simulations the value of the temperature at the node of each cell. To calculate the heat ux that ows through our domain, we need to sum up the ux that ows through each cell on the top boundary (for example). Let (i, j) be the node of a non-corner top boundary cell (see gure 2.8).

Figure 2.8: A representation of a top (non-corner) cell with its neighbor cells.

The temperature gradient at cell (i, 1) (top boundary cell for which j = 1) is given by:

∂T ∂z i = T i,1 -T i,2 z 2 -z 1 Recall that r i+ 1 2 = r i + r i+1 2 and r i-1 2 = r i + r i-1 2
so the ux that ows through the cell with node (i,1) is:

q i = π λ i,1 + λ i,2 2
∂T ∂z i r 2 i+ 1 2 -r 2 i-1 2 (2.11)
Finally, to calculate the total ux which ows through our computational domain we add up the ux at each top-border cell.

Heat ux with respect to water liquid volume

During the evaporation/condensation of water in a wet granular medium, the liquid volume in each elementary cell changes with time. It is interesting to know the variation of the eective thermal conductivity with respect to the liquid volume. By using the Fourier law on our computational domain, the eective thermal conductivity is proportional to the heat ux. This heat ux depends of course on the water liquid volume because, as stated in earlier, the thermal conductivity of the water is much more important than that of the dry air. Further, we expect a jump in the heat ux curve because there are two possible geometrical congurations for the liquid menisci (the contacting state and the non-contacting one, as described in gures 2.6 and 2.7).

The numerical values are the following:

thermal conduct. of water λ l = 0.61 W m -1 K -1 thermal conduct. of dry air λ a = 0.024 W m -1 K -1 thermal conduct. of solid λ s = 1.52 W m -1 K -1 radius of the grains R 0 = 1 × 10 -3 m surface tension of water γ = 72 × 10 -3 N m -1
The computed heat ux is presented in gure 2.10, for humidity (or liquid fraction, which is the ratio of the liquid volume divided by the void volume) ranged from 0 to 0.18. /R 0 = 0.1, θ = 30 deg. Note the great increase in the heat ux due to the presence of water: about a factor 5 when humidity is only 0.12.

As expected, a hysteresis behavior is highlighted by our numerical computations and is more visible in gure 2.11 which presents an enlarged view of gure 2.10. The small interior part of the with contact curve is never attained because it is physically unstable (see gure 2.12 for explanation). The smaller surface area corresponds to the stable meniscus because during dynamic evolution, liquid water tends to minimize its surface area. Therefore, the upper part of the curve should never be attained.

Once the heat ux through the computational domain has been obtained, it could be very easy to deduce the eective thermal conductivity λ e by using the Fourier law applied in a global manner to our domain. Actually λ e is proportional to the heat ux so it is not necessary to present the numerical results for λ e , as the main objective of this chapter is to focus on the hysteresis behavior.

The role of the parameters and θ is summarized in table 2.1. Referring to gure 2.11, we characterize the hysteresis by dening the two following quantities:

the Hysteresis Length which is the dimension of the hysteresis cycle along the x-axis; the Flux jump which is the other dimension along the y-axis.

Table 2.1 shows that the hysteresis is more and more pronounced as the distance increases. Moreover, it is more pronounced for small contact angles. In the extreme geometrical situation where the distance is zero (the solid spheres are contacted themselves) the hysteresis disappears because we obtain always the contacting state. chapter 2

Conclusion

A strong hysteresis behavior for the eective thermal conductivity has been revealed when changing the humidity of a granular medium. It is due to a switch between two dierent geometrical congurations of the liquid meniscus attached to the two solid spheres as a bridge. Presented as is, the eective conductivity law λ e is not isotropic (whereas it should be) because the meniscus' axis is supposed to be always aligned with the heat ux. In real situations, not only the distance between grains vary, but also the angular direction. DEM (Discrete Element Method which explicitly describes the dynamics of assemblies of particles) simulations could be used with benet, because they can give us important information about the statistic distributions of gap and contact angular direction in spheres' packings. By using a procedure close to that of Usowicz [START_REF] Usowicz | Thermal conductivity modelling of terrestrial soil media -a comparative study[END_REF], based on thermal resistance of regularly distributed spheres of dierent properties, we could combine the results obtained in the current chapter, parameterized in terms of (i) the gaps between the grains and (ii) the angular direction θ of the contact line between the grains, in order to predict the thermal conductivity at the macroscopic level.

In this study, the heat transfer has been restricted to the conduction phenomenon. It is evident that, according to the problem treated, other types of transfer should be taken into account, such as e. g. latent heat due to evaporation/condensation at the gas-liquid interface, or radiation. A simple version of a Stefan problem is the melting of a semi-innite slab of ice, initially at a temperature (T 0 < T f ), T f is the fusion temperature, the boundary which is raised initially at t=0 to a temperature above T f , at which it is subsequently maintained (T = T w ). An interface on which melting occurs, moves from the boundary to the slab and separates a region of water from one of ice at T f . In the literature, various techniques have been applied to treat the latent heat release for phase-change problems [START_REF] Voller | General source-based method for solidication phase change[END_REF][START_REF] Crank | Free and Moving Boundary Problems[END_REF]. Among them, two main classes of methods can be distinguished, the front-tracking methods and the xed-grid enthalpy methods. In the enthalpy methods, the phase-change conditions are implicitly accounted for in the denition of enthalpy.

Therefore, the solid/liquid interface does not need to be simultaneously tracked and is usually recovered afterward from the computation of isotherm. Compared to the fronttracking methods, the enthalpy methods are attractive due to the simple (implicit) representation of latent heat evolution at the interface, the ease of implementation, the capability of handling mushy phase change problems, and the wide range of applicability for multidimensional problems with complicated interface shapes.

This chapter is divided into two main parts. In the rst part, we present two wellknown enthalpy methods, the Latent Heat Accumulation method (LHA) and the Apparent Heat Capacity method (AHC) to solve fusion (melting of ice) problem in 1D coordinate system that was already presented in [START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] and we completed further on this work. The importance of such problem is the presence of an exact mathematical chapter 3

solution which was analyzed in many studies [CJ59, Ozi80, KC90]; the analytical solution presented below (see equations 3.1 and 3.2) for the case of fusion of ice assuming that ice and water have the same density ρ s = ρ l and that the physical properties in each phase are constant (heat capacity C and conductivity λ):

T l = T w + (T f -T w ) erf (x * ) erf x * sl where 0 < x * < x * sl (3.1) T s = T 0 + (T f -T 0 ) erfc α l αsx * erfc α l αsx * sl where x * sl < x * < ∞ (3.2)
where T 0 is the initial temperature, T w is boundary temperature, T f is the fusion or melting temperature,

x * = x 2 √ α l t
is the dimensionless position and α = λ ρC is the thermal diusivity which are calculated from corresponding thermal properties for both solid α s and liquid α l values. The values of dimensionless position of solid-liquid interface x * sl are obtained via the nonlinear algebraic equation: In his PhD thesis, M. Muhieddine presented an adaptive mesh in 1D coordinate system and applied it to the LHA method. In this work, we adapt this adaptive mesh to the AHC method in order to compare the performance and accuracy of the two numerical methods (LHA and AHC) with respect to the analytical solution.

T f -T 0 T f -T w λ s λ l α l α s exp -α l αs (x * sl ) 2 erfc α l αs x * sl + exp -(x * sl ) 2 erf x * sl - √ πx * sl L C l (T f -T w ) = 0 (3.
The second part of this chapter is devoted to Zohour, a node-based original adaptive 2D mesh algorithm developed by É. Canot. We use this adaptive mesh together with AHC method to solve the phase change problem (fusion) in 2D coordinate system.

Preliminary results are presented at the end of this chapter.

It is important to mention that all the results obtained in this chapter are obtained using:

Latent heat of fusion of ice: L = 338 × 10 3 J/kg The, so called, enthalpy approach is based on the fact that the energy conservation law, expressed in terms of energy (enthalpy) and temperature, together with the equation of state containing all the physical information needed to determine the evolution of phases. In the basic enthalpy scheme, enthalpy is used as the primary variable and the temperature is calculated from a previously dened enthalpy-temperature relation. As mentioned in [START_REF] Prapainop | Simulation of ice formation by the nite volume method[END_REF][START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF], the enthalpy-temperature relation is dened by:

T 0 = -10 • C and T w = 20 • C λ l = 0.556 W/
H = ρC s T T < T f (solid) ρC s T f + ρC l (T -T f ) + ρL T ≥ T f (liquid) (3.4)
where H is enthalpy, T is the temperature, T f is the temperature of fusion, ρ is the density of both liquid and solid (assumed equal), C s and C l represent the thermal capacity of the solid and liquid phases respectively and L is the latent heat of fusion.

In a solid to liquid phase transformation (i.e., fusion), the conservation of energy can be written down in terms of the mixture enthalpy H as:

∂H ∂t = ∇.(λ∇T ) (3.5)
where λ is the thermal conductivity.

Note that ice and water are assumed to have equal densities and thus there is no need to take into account the convection term because the displacement of the material is always null.

In LHA, the physical domain is divided into nite number of cells (control volumes)

where we apply the energy conservation to each control volume. The accumulated latent heat Q is initialized to zero and at each time step, we check the phase and temperature of each control volume to obtain a heat balance. If the control volume V i is in solid phase and its temperature T i is greater than the fusion temperature T f , the cell will be marked as a mushy cell [START_REF] Prapainop | Simulation of ice formation by the nite volume method[END_REF]. The temperature of the current time step is calculated using an explicit scheme, the temperature of the control volume is set to the fusion temperature and latent heat increment is calculated:

∆Q i = ρC m (T i -T f ) V i , where C m = (1 -θ)C l + θC s with θ is the fraction of solid ice in the mushy cell. The latent heat increment ∆Q i is added to Q i for subsequent time steps until Q i = Q total i where Q total i = ρLV i (total latent heat of cell i).
The cell becomes liquid and latent heat increment is no longer calculated. LHA implicitly contains the moving interface condition within the numerical model which renders it more exible than the front tracking methods.

In equation 3.5,

∂H ∂t

can be replaced by

∂H ∂T ∂T ∂t

and thus equation 3.5 can be written as:

ρC G ∂T ∂t = ∇.( λ∇T
) (3.6)
The Latent Heat Accumulation Method (LHA)

The spatial discretization of equation 3.6 is done using the nite volume method (see equation (3.7)) where the subscripts G, E and W represents respectively the current cell, and its east and west neighbors (see gure 3.3).

∂T ∂t = - λ W + 2λ G + λ E 2 (ρC) G (∆x) 2 T G + λ G + λ E 2 (ρC) G (∆x) 2 T E + λ G + λ W 2 (ρC) G (∆x) 2 T W (3.7)
y n+1 = y n + hf (t n , y n ) (3.9)
The spatial temperature prole presented in gure 3.4 taking ∆x = 6.25 × 10 -3 m presents the analytical solution in red and the numerical solution in blue. We notice that the dierence between the two solutions all over the domain is more noticeable close to the phase change area upon zooming. In gure 3.5, we represent the variation of the temperature as function of time at xed position (x = 5 cm). We notice some uctuations in the numerical solution especially at the interface which is normally due to the fact that the temperature is blocked in the mixed cell during a certain time and then released.

In our computation, the accuracy of the numerical solutions is measured against the analytical solution. It is known that any measure of accuracy involves a choice of particular error norms. To measure the error in temperature prole at nal time (i.e. t is xed), we dene the following L 2 -average norm error:

e x = i [T i -T exact (x i )] 2 ∆x i 1 2 l (3.10)
Here T i is the numerical solution of temperature at grid point x i at t = t f inal = 50 h, T exact (x i ) is the exact analytical solution at x i and t f inal and l is the length of the computational domain. Analytical solution in red and numerical solution in blue using LHA method. The average L 2 relative error over the whole domain is 1.82 × 10 -2 %. In this gure we represent 1 8 of the whole computational domain.

To measure the error in temperature history at a precise position x (here we take x = 5 cm), we dene the following L 2 -average norm error:

e t = k T k -T exact (t k ) 2 ∆t k 1 2 t max (3.11)
where T k is the temperature at grid point x at t = t k , T exact (t k ) is the exact analytical solution at x and at t k and t max is the total time of the simulation .

Some authors present a global L 2 error: e

= i k T k i -T exact (x i , t k ) 2 ∆t k ∆x i 1 2
.

By separating e into two dierent parts, i.e. e x and e t , we focus on the inuence of time and space discretization. Of course, we will check that both errors e x and e t approaches zero as ∆x and ∆t become smaller and smaller.

The Apparent Heat Capacity Method (AHC)

As the most popular and perhaps the simplest enthalpy method, the so-called Apparent Heat Capacity method [START_REF] Bonacina | Numerical solution of phase-change problems[END_REF] (AHC, sometimes also called eective specic heat method) has some special computational attributes. In the AHC method the enthalpy function is used for introducing the apparent heat capacity only, but it is not present in the resulted governing equation explicitly. Therefore the AHC method formulation of the heat transport with phase change has the same form as that when there is no phase change; and the only dierence is in the denition of the apparent heat capacity, which includes the Dirac function and is usually approximated by some smoothed functions.

Furthermore, the AHC method is the easiest to implement into existing general heat transfer programs. It is known that the computation of the apparent heat capacity is of key importance in the implementation of the AHC method.

σ = χ (T -T f ) (3.12)
is representative of the volumetric fraction of the solid phase (ice). χ is a step function whose value is zero when T < T f but equals one otherwise. Its derivative, i.e. the variation of the initial phase fraction with temperature, is a Dirac delta function [START_REF] Snieder | A Guided Tour of Mathematical Methods for the Physical Sciences[END_REF]:

dσ dT = δ (T -T f ) (3.13)

Numerical Strategy

To alleviate the singularity of the Dirac delta function, it can be approximated by a rectangular function or a normal distribution function over a nite temperature interval. Using a normal distribution function, the Dirac can be transformed to a smooth variation:

dσ dT = 1 ∆T 2 π exp   - √ 2 ∆T 2 (T -T f ) 2   (3.14)
in which ∆T is the assumed phase change interval. Integrating equation (3.14) provides an approximation of the volumetric fraction of the solid phase σ:

σ(T ) = 1 2 1 + erf √ 2 ∆T (T -T f ) (3.15)
As a consequence of all what preceded, the physical properties (C and λ) provided by Bonacina could be regularized as mentioned in [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] (see gure 3.6 continuous lines and equations 3.16 and 3.17).

C = C s + (C l -C s )σ + L ∂σ ∂T (3.16) λ = λ s + (λ l -λ s )σ (3.17)
The heat conduction equation used to model the heat transfer with phase change is given by Bonacina et al.:

(ρC) e (T ) ∂T ∂t = div (λ e ∇T ) (3.18)
where (ρC) e is the eective volumetric heat capacity and λ e is the eective thermal conductivity.

We use the method of lines to solve the heat diusion equation (PDE). The spatial discretization is performed using the vertex-centered nite volume method where the spatial variable is discretized into N discretization points and each state variable T is transformed into N variables corresponding to its value at each discretization point.

After spatial discretization, we endup having a system of ordinary dierential equations (see equation (3.19)) which could be solved by a simple Euler explicit scheme or by an implicit automatic ODE solver. We choose to apply both techniques in order to compare their eciency.

dT dt = A(T )T (3.19)
The ODE coecient matrix A(T ) has a tridiagonal structure due to the 1-D Laplacian discretization. We choose an ODE solver based on a BDF implicit scheme which possesses the property of stability and therefore does not suer from the stability step size constraint [START_REF] Hindmarsh | Description and use of LSODE, the Livermore solver for ordinary dierential equations[END_REF][START_REF] Yao | Application of the homographic approximation in the enthalpy method for phase-change problems[END_REF]. The BDF implicit scheme requires the calculation of a Jacobian matrix which is calculated and generated by a Computer Algebra System (CAS, Maple or Maxima) which generates automatically an optimized Fortran code.

The Jacobian matrix is then stored in a sparse format. Note that the numerical calculation is performed with ddebdf routine of the SLATEC Fortran library which was modied to use the UMFpack sparse linear solver [START_REF] Davis | Algorithm 832: Umfpack, an unsymmetric-pattern multifrontal method[END_REF]. The ODE solver performs time integration by adjusting automatically the time step in the BDF scheme (order is also automatic and varies between 1 and 5) and all these primary numerical libraries are grouped in the easy-to-use MUESLI library [Canb]. Note that some very ecient linear solver can be found for a tridiagonal system but the UMFpack solver will be required for more complex matrices in the coming chapters.

Choice of ∆T

The choice of the value of the phase change temperature interval ∆T has a big inuence on the accuracy of the solution of the heat equation using AHC method. We plot the value of ∆T which corresponds to the minimum quality factor as function of the mesh size h. This value of ∆T is called ∆T optimum (see gure 3.8).

To understand the eect of this parameter on the temperature history, we consider 3 dierent values of ∆T : ∆T = ∆T optimum 4

, ∆T = ∆T optimum and ∆T = ∆T optimum × 4. Figure 3.9 represents the temperature history with these 3 dierent values of ∆T together with the analytical solution. In gure 3.10, we notice that the numerical solutions with ∆T optimum and

∆T

Results using AHC method with uniform mesh

The spatial temperature prole presented in gure 3.13 taking N = 640 mesh cells presents the analytical solution in red and the numerical solution (using Euler explicit scheme) in blue. We notice that there is a slight dierence between the two solutions all over the domain especially around the phase change region. In gure 3.14, we represent the variation of the temperature as function of time at xed position (x = 5cm) which could be considered as a typical representative position of phase change phenomenon.

We notice that there are less uctuations in the numerical solution compared to those present using LHA with same mesh size. Using AHC method with an implicit ODE solver, we obtain a temperature history with slightly less uctuations than that using AHC method with Euler explicit scheme and with slightly smaller L 2 relative errors in both the temperature prole and temperature history (see gures 3.15 and 3.16). On the other hand, the calculation cost is higher using an ODE solver. In order to get rid of the uctuations in the temperature history using both methods and to reduce the error between the analytical and numerical solutions, we need to use an adaptive mesh renement technique which will be applied to both LHA and AHC methods.

Adaptive Mesh Renement

Gridding techniques fall into two mains categories, xed gridding technique and adaptive gridding technique.

The rst technique involves cells of equal size and if renement is needed to reach high precision, global renement technique is used in a way to preserve the uniform-mesh criterion. Using ne grid over the entire domain can be computationally intensive both in terms of CPU time and memory requirements.

The latter falls into two broad classes, adaptive mesh redistribution (moving mesh method) and adaptive mesh renement. Moving mesh schemes employ a xed number of grid points and attempt to dynamically move the grid points to areas of high error in the mesh. Moving mesh schemes can't be easily implemented and suer from many deciencies. They can be unstable and can result in mesh tangling and local degradation of the solution. They can never reduce errors below an arbitrary limit, and they often fail when time-dependent boundary conditions are enforced, as they are incapable of handling the migration of regions containing irregularities or singularities in the solution as it evolves in time [START_REF] Oden | Advances in adaptive methods in computational uid mechanics[END_REF].

Adaptive mesh renement is an intelligent alternative to a global renement approach because it requires fewer resources to achieve the same level of solution accuracy. Con-versely, better accuracy can be achieved with the same resources by adaptively rening the mesh in only the important regions. Such schemes present a dicult data management problem, since they involve a dynamic regeneration of the mesh, renumbering of grid points, cells or elements, and element connectivities as the mesh is rened which is inconvenient for many advanced solvers. On the other hand, they are ecient in terms of accuracy/CPU time and can be also used to coarsen a mesh i.e. use less number of mesh cells and therefore reduce the number of unknowns when the error is lower than an assigned lower-bound tolerance [START_REF] Oden | Advances in adaptive methods in computational uid mechanics[END_REF]. Figure 3.17 summarize all these information.

Gridding Techniques

Fixed gridding technique

Global Refinement: High computational cost (CPU/memory), excessive number of interpolations.

Adaptive gridding techniques

Moving Mesh: Local refinement, fixed number of cells, unstable, easily implemented.

Adaptive refinement mesh: local refinement, Efficient (precision/CPU), few interpolations. Taking into consideration all what preceded, the adaptive mesh renement has many advantages compared to other mesh techniques. An adaptive mesh could be applied over a cell-based non-conforming scheme, a cell-based conforming scheme or a node-based (Voronoi) scheme.

Cell-based schemes store solution unknowns at cell centroids and impose conservation on the cell itself whereas node-based schemes store solution unknowns at the cell vertices and impose conservation on the cell-duals.

A mesh is said to be conforming if adjacent cells intersect only at either a common edge or a common vertex or an entire face (see gure 3.18 for a non-conforming mesh which does not exist in 1 D). Rening a non-conforming cell-based mesh is easier than a chapter 3 conforming cell-based mesh but a numerical scheme is easier to be applied to the latter than the former [SLT + 12].

The other type of adaptive mesh renement is a node-based adaptive mesh based on Voronoi tessellation (see gure 3.19), the renement technique in this kind of mesh is of moderate diculty whereas a numerical scheme can be easily applied. The adaptive mesh that we chose to apply to our numerical methods (LHA and AHC) in 1D is of node-base type where we add and remove nodes for the initial mesh.

The starting point is a xed basic uniform mesh over the whole domain, we rene the basic mesh by a number of subdivisions near the phase change interface in a way that at each level of subdivision we add a node at the middle of two successive nodes. As the interface moves with time to another 'cell', the previously added nodes are removed and new nodes are added at the interface node. Our rened zone is dened as a small interval around the phase change temperature (0 • C in this case) so the mesh points need to be concentrated around this interval and dynamically adjusted to follow the change in temperature as it propagates in time. In this adaptive mesh technique, during time evolution all nodes are xed in position and this guarantees a good accuracy due to small number of needed interpolations (interpolation is needed only when a new node is inserted).

LHA with our adaptive mesh

The simplest way to avoid the presence of non-physical behavior (uctuations) resulting from the use of the enthalpy formulation is to reduce the spatial step size. However, if this is done uniformly over the whole domain then the overall method would be computationally expensive. The mesh spacing needs only to be rened around the position of the moving phase change boundary which suggests that some form of adaptive mesh algorithm would be useful. We chose to adapt our adaptive mesh to LHA method in 1D so rst the melting front is initially located in the rst cell and can move anywhere along the discretized domain. When the phase-change front transfer to a new cell (mushy cell), the node added in a previous step is removed and a new node will be added to the new element. This technique was implemented as a double linked list in Fortran 95 as explained in [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] where each item of this list contains information about geometry of the cell and physical variables (temperature, Latent heat accumulation necessary in LHA method to detect the phase change cell and a tag for the phase state). The user chooses the number of recursive subdivisions at the beginning of the program depending on the needed precision. chapter 3

Results and consequences

The use of an adaptive mesh in LHA method has enhanced the accuracy of the solution.

Figures 3.21 and 3.22 represent the numerical results obtained using 80 nodes in the basic mesh and 3 subdivision levels. The calculated errors for both spatial temperature and temperature history are respectively 1.66 × 10 -2 % and 4.05 × 10 -1 % which are approximately equal to those obtained by a uniform mesh with 640 nodes where as the CPU time is reduced by approximately a factor of 5. In gures 3.23 and 3.24, we notice that for a constant number of subdivisions, the error between the numerical and analytical solutions decrease as the total number of mesh cells increase which insures the consistency of our method. The error in spatial temperature prole (e x ) decreases from 6.37 × 10 -2 % to 1.18×10 -3 % as the basic number of mesh cells (number of subdivisions = 2) increase from 80 to 1280 whereas it decreases from 1.66 × 10 -2 % to 7.16 × 10 -4 % for number of subdivisions = 3. The error in temperature history (e t ) at x = 5 cm decreases from 9.33 × 10 -1 % to 3.21 × 10 -2 % as the basic number of mesh cells (number of subdivisions = 2) increase from 80 to 1280 whereas it decreases from 4.05 × 10 -1 % to 1.41 × 10 -2 % number of subdivisions = 3.

In general, if we take into consideration the precision (error) and CPU time (see gures 3.25 and 3.26), we notice that if we seek low cost (small CPU time), we must choose small number of basic mesh cells (N = 40 green curve) with few subdivisions.

On the other hand, if we seek high precision, we must choose relatively large number of basic mesh cells (N = 640 magenta curve) with few subdivisions.

Actually, the error using an adaptive mesh is given by:

error ≤ k|h max | a
where h max is the biggest mesh size, k is a constant independent of h max but dependent on the numerical function. a is a constant which depends on the regularity of the function. h max remains constant as the number of subdivisions varies so it is not worth using large number of subdivisions to get high accuracy.

Adaptive mesh applied to AHC method

We adapt also our adaptive mesh to the AHC method in 1D. As in LHA, we need to rene the mesh around the position of the moving phase change boundary.

If an explicit scheme is used to solve the ODE system, then the application is easy, we search the rst node that attains a temperature less than melting temperature and we start rening the cell just before it by adding new nodes while we keep checking the temperature of the nodes. If the temperature of the rened zone is less than the current node temperature then the rened zone is shifted to the right (i.e. we delete added nodes from previous cell and add new nodes to the current cell).

If an automatic ODE solver is used, we use the same procedure to detect the zone to be rened but we added a stopping condition well adapted to our phasechange problem: the time integration is stopped when the temperature of the cell next to the right of the 'rened zone' lies in a small interval around the melting temperature then we start a new time integration and we continue in the same way until the nal time is attained.

It is important to mention that the multiple restart of the ODE solver results in high CPU time.

Results and consequences

The use of an adaptive mesh in AHC method has enhanced the accuracy of the solution.

Figures 3.27, 3.28, 3.29 and 3.30 represent the numerical results obtained using 80 nodes in the basic mesh and 3 subdivisions. We use an explicit scheme and an ODE solver based on a BDF implicit scheme respectively. Using Euler explicit scheme, the calculated errors for both spatial temperature and temperature history are respectively 1.92 × 10 -2 % and 3.84 × 10 -1 % which are approximately equal to those obtained by a uniform mesh with 640 nodes where as the CPU time is reduced by approximately a factor of 5. Using an ODE solver based on a BDF scheme, the calculated errors for both spatial temperature and temperature history are respectively 2.06 × 10 -2 % and 3.47 × 10 -1 % which are approximately equal to those obtained by a uniform mesh with In gures 3.31, 3.32, 3.33 and 3.34, we notice that for a constant number of subdivisions, the error between the numerical and analytical solutions decrease as the number of basic mesh cells increase which insures the consistency of our method. As in LHA with adaptive mesh, we take into consideration the precision (error) and CPU time to check the eciency of AHC method using an explicit Euler scheme and an ODE solver based on a BDF implicit scheme respectively (see gures 3.35, 3.36, 3.37 and 3.38). If high precision is requested then we must choose relatively big number of basic mesh cells (blue or magenta curves) with few subdivisions whereas if low cost (CPU time) is requested then we can choose small number of basic mesh cells (green curve) with few subdivisions.

3.3.3

Comparing the performance of the three methods using adaptive mesh renement Figures 3.39 and 3.40 represents the variation of L 2 error in temperature history (e t ) at x = 5 cm as function of CPU time for LHA method with explicit scheme, AHC method with explicit scheme and AHC method with an automatic ODE solver (BDF implicit scheme). For dierent number of basic mesh cells and dierent subdivisions, LHA ended up to be the most performant method concerning both precision and cost. Concerning AHC method, as expected, an explicit scheme is more performant than an implicit BDF scheme for small number of mesh cells because ∆x is relatively big and thus the time step also can be chosen relatively big. On the other hand, AHC method with implicit BDF scheme is more ecient than that with explicit scheme for large number of mesh cells. This last gure shows that, for constant level of subdivisions, LHA is basically the most ecient method except when N basic = 80 where AHC with ODE solver (implicit scheme) is more performant if precision/cost is taken into consideration.

Part II: Phase change problem in 2D coordinate system

In this part, we present an adaptive 2D mesh Zohour and adapt it to the AHC method to solve the phase change problem (fusion) in 2D coordinate system.

The physical problem consists of melting a rectangular slab of ice, of dimension l 0 ×l 0 such that l 0 = 25 cm, initially at a temperature T 0 < T f . We raise the temperature of the left top boundary to T w (T w > T f ), at which it is subsequently maintained.

The energy equation in 2D coordinate system is given by:

(ρC) e ∂T ∂t = ∇.(λ e ∇T ) (3.21)
with the following initial and boundary conditions: 2 cm, 4 cm and 6 cm respectively for a 500 × 500 mesh.

T (x, y, 0) = T 0 for 0 ≤ x ≤ l 0 and 0 ≤ y ≤ l 0 T (x, y, t) = T w at y = 0 for 0 ≤ x ≤ l 0 2 T (x, y, t) = T 0 at x = l 0 for 0 ≤ y ≤ l 0 and at y = l 0 for 0 ≤ x ≤ l 0 ∂T ∂y = 0 at y = 0 for l 0 2 < x ≤ l 0 ∂T ∂x = 0 at x = 0 for 0 < y ≤ l 0

Zohour

Zohour is a node-based adaptive 2D mesh algorithm written and developed by É. Canot [Canc] in Fortran. This algorithm denes, at the interior of a rectangular domain, a set of nodes (identied by an ID and coordinates (x, y)). The mesh is then created by building the Voronoi set of cells.

The initial mesh is composed of regular squared set, at each level of subdivision, In fact, at each level of subdivision, only nodes are added or removed according to some prescribed rules and verifying some criteria in a way that we conserve maximum information i.e. the initial mesh remains xed and thus very few interpolations are needed which reduces the error.

new
There exists a function f which is associated to a certain real value at each point of the domain. This function can be the result of the solution of a PDE (it corresponds to temperature T in our case). The mesh once constructed can be updated according to one or more criteria: The main criteria is the Hessian matrix based on the curvature of the function.

As the Hessian becomes stronger in a certain zone, a higher level of renement is needed in that zone (i.e. the creation of smaller mesh cells).

Another criteria could be the phase change zone in a problem related to phase change. A higher level of renement is needed close to the interface position.

Here we are interested in Finite Volume method of node-based type. The common approach of all nite volume methods is to approximate integral of normal component of the ux q on a set of nite (control) volumes V i , i=1,2,..,n.

dV i (q.n) ds ≈ j∈Σ(i) q i,j
Here, the nodes contain the value of a function T (x, y) (temperature). From these nodes, we construct a Voronoi mesh (partitioning the 2D plane by all the perpendicular bisectors of two neighboring nodes). If the nodes are distributed on equidistant horizontal and vertical lines, the Voronoi diagram is a square grid whose sides are horizontal and vertical. A Voronoi mesh has a very interesting property: each side of a mesh cell is located on the perpendicular bisector of two neighboring nodes, thus the calculation of the gradient of T (x, y) normal to the side can be easily performed in an ecient and accurate manner regardless the shape of the cell.

The local renement technique presented in Zohour is similar to that of HOMARD [EDF] but the advantages of Zohour are: the eciency of Hessian estimation which is composed of two parts: nding the position of the neighboring nodes and its calculation.

the fact that unlike most adaptive algorithms which have a renement factor of 2, Zohour meanwhile has a factor of √ 2, so it is more progressive.

Zohour and AHC method

We use the AHC method to solve the heat diusion equation in 2D coordinate system.

The implementation of the AHC method is already described in the previous part together with numerical technique used the method of lines. The spatial discretization is performed using a vertex-centered nite volume scheme where the spatial variable is discretized into N × N discretization points and each state variable T is transformed into N 2 variables corresponding to its value at each discretization point.

We end up with a system of ODEs that we chose to solve (for simplicity) with an Euler explicit scheme. The explicit scheme has the disadvantage to be conditionally stable (∆t should not be greater than a critical value ∆t c (see equation 3.22) ) but it can be easily implemented (no matrix resolution).

∆t c = 1 2 ρ s C s ∆x 2 ∆y 2 λ s (∆x 2 + ∆y 2 ) (3.22) chapter 3
Suppose that T n i is the temperature of cell i at time t n , t 0 = 0 , t n+1 = t n + ∆t and T 0 i is given, then the explicit scheme to solve the heat diusion equation can be written as:

T (n+1) i = T (n) i + ∆t (ρC) e,i A i C λ e,i ∂T i ∂n dl (n) (3.23)
where A i is the area of the cell i, (ρC) e,i and λ e,i are the eective volumetric heat capacity and the eective conductivity calculated using the AHC method at the cell i.

C

∂T i ∂n dl is the line integral representing the sum of uxes on the sides of cell i.

Highlighting the eciency and accuracy of Zohour

In this section, we compare the results obtained by solving the fusion problem using a regular mesh and using Zohour.

In Zohour, we can control whether renement is done at each time step ∆t or each 20, 50, ... time steps. Numerical trials has showed that the overall CPU time does not vary a lot if re-meshing is performed each 20 or 50 or 75 time steps due to the fact that the CPU time of the explicit update of the vector is negligible compared to that required to calculate the Hessian function. We represent in gures 3.46 and 3.47, the variation of the L 2 error in both temperature prole at the nal time and the temperature history at a xed position as function of the CPU time. In both gures, the level of subdivision chosen is two and re-meshing is performed at each, 20, 50 and 75 time steps. The accuracy/cost curves showed that it is not worth it to perform re-meshing at each time step or for large number of time steps. In the rst case, the CPU time is greater and in the second case we risk to lose accuracy. For all what follows re-meshing will be done each 20 time steps.

Recall that in our computation, the accuracy of the numerical solutions obtained by Zohour is measured against the reference solution. To measure the error in temperature prole at nal time (i.e. t is xed), we dene the following L 2 error:

e x =   i j [T i,j -T ref erence (x i , y j )] 2 A i,j   1 2 (3.24)
where T i,j is the numerical solution of temperature at grid point (x i , y j ) at t = t f inal = 20 h, T exact (x i , y j ) is the reference solution at (x i , y j ) and t f inal and A i,j is the area of the cell whose node is located at (x i , y j ). To measure the error in temperature history at a precise position (x, y) (here we chapter 3 take x = 0 cm and y = -2 cm), we dene the following L 2 error:

e t = k T k -T ref erence (t k ) 2 ∆t k 1 2 (3.25)
where T k is the temperature at grid point (x, y) at t = t k , T ref erence (t k ) is the reference solution at (x, y) and at t k . Finally, ∆t k = t k+1 -t k .

As we mentioned in the previous section, the renement takes place next to the phase change area and in the region where the Hessian attains a large value. This could be clearly seen in gures 3. We notice that for same number of mesh cells, the error decrease as the level of subdivisions increase. In gures 3.52 and 3.53, we consider the temperature history for three dierent positions ((0, -2 cm), (0, -4 cm), (0, -6 cm)), the continuous lines represent the reference solution where as the + represent the numerical solution taking same number of basic mesh cells (70 × 70) for respectively one level and three levels of subdivisions.

Conclusion

In the rst part of this chapter we tackled the phase change problem (Fusion) in 1D

coordinate system where we presented two methods (LHA and AHC) to deal with the phase change. Moreover, an adaptive mesh was adapted and tested over the two methods to conclude nally that the LHA method is more competitive than the AHC method.

In the second part, we presented a new 2D adaptive mesh algorithm Zohour and used it to solve the 2D phase change (fusion) problem. The algorithm is tested against a reference solution obtained using a regular ne rectangular mesh. Results obtained are very interesting and the algorithm is quite ecient. Nevertheless, the results obtained are preliminary and further work is required.

1. Improve the performance of this new adaptive 2D mesh algorithm.

2. Study the eect of the variation of the parameters related to the use of Zohour:

The Hessian coecient and the thickness of the phase change zone.

Chapter 4

Estimation of the thermophysical properties of the soil by Inverse problem: theoretical aspects

The work presented in this chapter is motivated by the studies of agricultural and archaeological soils. The investigation of the thermal properties of the soil can have signicant practical consequences such as evaluation of optimum conditions for plant growth and development and can be utilized for the control of thermal-moisture regime of soil in the eld [START_REF] Olandunjoye | Variability of soil thermal properties of a seasonally cultivated agricultural teaching and research farm[END_REF]. These properties inuence how energy is partitioned in the soil prole so the ability to monitor them is a tool to manage the soil temperature regime that aects seed germination and growth. It can also provide information about the use of re by ancient civilizations whether for cooking, heating or self-defense ... [START_REF] Eckmeier | Detecting prehistoric re-based farming using biogeochemical markers[END_REF]. A systematic application of numerical modeling in a particular eld of agriculture and archeology which is the study of seed germination and archaeological hearths is presented.

Our porous media of interest, the soil, is composed of a solid matrix (s) which is composed of diverse minerals and organic components and the void space (pores) lled by uids (f ) which is may be water liquid or gas or both. If all soil pores are lled with water the soil is said to be saturated (there is no air left in the soil). Another case is that of dry soil when the pores are lled with gas (i.e. no liquid water is present). The physical properties of the uid: liquid water(l), water vapor (v) and air (a) are known: volumetric heat capacity ((ρC ) l , (ρC) v and (ρC) a ) and thermal conductivity (λ l , λ v and λ a ). The unknowns are those of the solid matrix ((ρC ) s and λ s ).

In this chapter, we introduce a numerical strategy in both 1D and 3D-axisymmetric coordinate systems to estimate the thermophysical properties of the soil solid matrix (volumetric heat capacity (ρC) s , thermal conductivity λ s and porosity φ) of a saturated porous medium where a phase change problem (liquid/vapor) appears due to intense heating from above. Usually φ is the true porosity, however when the soil is not saturated (which should concern most cases), φ may be taken equal to the part of water in the pores. This is of course an approximation which is correct for the energy balance; however, a case too far from the fully saturated conguration should be avoided because chapter 4

capillary and gravity forces are neglected in our current model. An unsaturated model limited to water in pendular regime will be outlined at the end of this work (see chapter 6).

The inverse problem, presented here, consists of the estimation of thermophysical properties of the soil knowing the heating history curves at selected points of the altered soil [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF]. In general, the mathematical formulation of inverse problems leads to models that are typically ill-posed [START_REF] Engl | Nonlinear inverse problems: theoretical aspects and some industrial applications[END_REF]. In such problems, we usually minimize a discrepancy between some experimental data and some model data [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF]. In our problem, we use the least square criterion in which the sensitivity coecients appear and where we try to minimize the discrepancy function which is expressed as the norm of the dierence between the experimental temperature and the numerical data obtained by our approximated model [START_REF] Björck | Numerical methods for least squares problems[END_REF]. The system composed of the energy equation together with other boundary initial problems resulting from dierentiating the basic energy equation with respect to the unknown parameters must be solved [START_REF] Majchrzak | Identication of substitute thermal capacity of solidifying alloy[END_REF].

At the stage of numerical computations, the Damped Gauss Newton method and the Levenberg Marquardt algorithm are used to minimize the least square criterion; that requires the solution of a system of highly nonlinear ordinary dierential equations. We propose a global approach similar to that presented by Muhieddine et al. [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF] using also the apparent heat capacity method (AHC see chapter 3) to deal with the phase change problem. It is important to note that in our new conguration, the solution is reached after taking into consideration the temperature history at selected points of the domain and at dierent time steps which was not the case in [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF] where the authors reached the solution by taking the temperature history at the nal time only and at all the points in the computational domain. This new conguration is more close to the experimental available data where the temperature is recorded at few sensors every 10 seconds for the whole time duration (typically few hours). Also, we used a scaling technique which sounds to be inevitable since our unknown parameters have dierent order of magnitude. We used the method of lines, where time and space discretizations are considered separately. The space discretization is performed using a vertex-centered nite volume method; the discretization in time is performed via an ODE solver that uses a backward dierentiation method (BDF) which is an implicit method for numerical integration of dierential equation that requires the calculation of the Jacobian matrix.

The code verication stage is based on the comparison between the numerical values of the parameters obtained by the inverse problem and those used to create the synthetic data. The advantage of our conguration to that presented in [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF] is that we propose a model which is more realistic and closer to the experimental setup i.e. our synthetic data consists of the calculation of the temperature at few sensors (around 5)

during the whole heating duration.

Identication of the thermophysical properties of a saturated soil during phase change in 1D

4.1 Forward problem (1D)

The physical problem consists of heating a saturated soil by a re. To model this problem, we replace the soil by a perfect porous medium in 1D nite domain of length l = 10 cm with constant and uniform properties heated from above by a constant temperature T c (temperature of the re between 300 • C and 700 • C). T c must be greater than T v (the evaporation of water or phase change temperature which is normally 100 • C).

There is a notable dierence in the temperature history when modeling heat transfer only on one hand and heat transfer coupled with steam ow on the other hand [START_REF] Muhieddine | Coupling heat conduction and water-steam ow in a saturated porous medium[END_REF]. In spite of this dierence and for simplicity reasons, we use the energy equation neglecting the convection term to model heat conduction transfer in the soil:

(ρC) e ∂T ∂t = div (λ e ∇T ) (4.1)
with the following initial and boundary conditions:

T (x, 0) = T 0 in Ω At x = 0 : T (x, t) = T c for t ∈ (0, t end ] (Dirichlet) At x = l : ∂T ∂x = 0 for t ∈ (0, t end ] (Neumann)
where T represents the temperature, T 0 is the initial temperature at t 0 = 0 (T 0 = 20 • C), ρ is the density, C is the specic heat capacity, λ is the conductivity, φ is the porosity (φ = 0.2 in all what follows), the subscripts e, f and s indicate the equivalent parameters of the medium, the properties of the uid and the porous matrix properties respectively.

Note that the thermophysical properties of the uid are temperature dependent and that is why the problem is nonlinear.

The eective volumetric heat capacity and the eective conductivity are dened by the equations (4.2) and (4.3):

(ρC) e = φ(ρC and Φ = 0.058. To avoid the tracking of the interface of the phase change problem (liquid/vapor) which appears when the water existing in the soil turns into gas, the Apparent Heat Capacity (AHC) method is used because it allows a continuous treatment of a system involving phase transfer. In this method [BC73, MCM09], the latent heat is taken into account by integrating the heat capacity over the temperature and the computational domain is considered as one region (see chapter 3). The Dirac delta function, representing the equivalent heat capacity, can be approximated by the normal distribution:

) f + (1 -φ)(ρC) s (4.2) λ e = λ f φ + (1 -φ) β λ f λs + Φ (4.
dσ dT = 1 ∆T 2 π exp   - √ 2 ∆T 2 (T -T v ) 2   (4.4)
where T v is the phase change temperature. The integration of equation (4.4) yields the error function approximations for the initial phase fraction.

σ(T ) = 1 2 1 + erf √ 2 ∆T (T -T v ) (4.5)
The functions dened by equations (4.4) and ( 4

C f = C l + (C v -C l )σ + L dσ dT (4.6) λ f = λ l + (λ v -λ l )σ (4.7) ρ f = ρ l + (ρ v -ρ l )σ (4.8)
where L is the latent heat of evaporation of liquid water (L = 2.256 × 10 6 J/kg) and the subscripts l and v indicate respectively the properties of the liquid water and that of the water vapor at 100 • C (ρ l = 1000 kg/m 3 , C l = 4000 J/kgK, λ l = 0.6J/kg, ρ v = 0.8 kg/m 3 , C v = 2000 J/kgK and λ v = 2.5 × 10 -2 J/kg).

Numerical method

We need to solve the heat diusion equation (PDE) so we choose the method of lines which is a way of approximating PDEs by ODEs where space and time discretizations are considered separately. The spatial discretization is performed using the vertex-centered nite volume method which conserves the mass locally and preserves continuity of uxes.

To apply the spatial discretization, the computational domain is divided into a nite volume grid or mesh with equal length h = ∆x as shown in gure 4.2. ) are computed as exactly the middle of two consecutive nodes, i.e.

x i+ 1 2 = 1 2 (x i + x i+1 ). The spatial derivatives are approximated by using a nite volume formula on three points so we end up with a semi-discrete system of ODEs which can be written in the form:

dT dt = B(T )T (4.9)
The ODE coecient matrix B(T ) has a tridiagonal structure due to the 1-D Laplacian discretization. As mentioned in chapter 3, the precision of computation in AHC method is sensitive to the value of phase change temperature interval ∆T chosen to approximate the dirac distribution. In fact, our phase change problem becomes more and more sti as the value of ∆T approaches zero. The diculty with sti problems is the prohibitive amount of computer time required for their solution by classical ODE solution methods.

The best scheme to solve such sti problems is the BDF scheme which possesses the property of stability and therefore does not suer from the stability step size constraint [START_REF] Hindmarsh | Description and use of LSODE, the Livermore solver for ordinary dierential equations[END_REF]. The BDF implicit scheme requires the calculation of a Jacobian matrix which is calculated and generated by a Computer Algebra System (CAS, Maple or Maxima) which generates automatically an optimized Fortran code. The Jacobian matrix is then stored in a sparse format. Note that the numerical calculation is performed with ddebdf routine of the SLATEC Fortran library which was modied to use the UMFPACK chapter 4

sparse linear solver [START_REF] Davis | Algorithm 832: Umfpack, an unsymmetric-pattern multifrontal method[END_REF]. The ODE solver performs time integration by adjusting automatically the time step in the BDF scheme (order is automatic and varies between 1 and 5) and all these primary libraries are grouped in the easy-to-use MUESLI library [Canb].

Inverse problem (1D)

In order to solve the parametric inverse problem consisting of nding the volumetric heat capacity (ρC) s , the conductivity λ s and the porosity φ of the solid matrix of a saturated soil, it is necessary to know the values of temperature T f gi at selected points (sensors) of the porous medium domain for times t f (external time-steps dierent from internal time-steps calculated automatically by the ODE solver): T f gi = T g (x i , t f ) where i = 1, 2, ..., M and f = 1, 2, ..., F . M and F are the total number of sensors and time steps respectively. We use the least squares criterion to solve this inverse problem so we try to nd the soil parameters that minimize the cost function which is dened by:

S((ρC) s , φ, λ s ) = 1 2 T f i -T f gi 2 2 (4.10)
where T f i = T (x i , t f ) are the temperatures being the solution of the direct problem for the assumed set of parameters at the point x i , i = 1, 2, ..., M for the time t f , f = 1, 2, ..., F and T f gi is the measured temperature at the same point x i for time t f .

It is important to mention that the authors in [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF] calculated the temperature at the nal time only and at all the points of the domain.

Parameter scaling

In [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF], the authors mentioned that the heat equation is not sensitive to the heat capacity in comparison to the other parameters and that it stagnate at its initial guess.

In reality, this is due to the fact that the parameters we are investigating are of very dierent magnitudes so it is necessary to perform parameter scaling or otherwise many searches would not converge. Gradient search techniques generally require parameter scaling to obtain ecient search convergence [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF].

The rst basic rule of scaling is that the variables of the scaled problem should be of similar magnitude and of order unity in the region of interest. If typical values of the variables are known, a problem can be transformed so that the variables are all of the same order of magnitude. The most commonly used transformation is of the form :

p = Dp (4.11)
where p is the vector of original variables p j , p is the vector of scaled variables pj and D is a constant diagonal matrix whose diagonal elements are set to be equal to the order of magnitude of its corresponding variable. We have to keep in mind that when the variables are scaled then the derivatives of the objective function are also scaled [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF].

Method of resolution

To illustrate the method of resolution, we dene the following vectors:

T g =                  T 1 g1 ... T F g1 T 1 g2 ... T F g2
where T g is the vector of measured temperatures, g(p (k) ) is the vector of temperatures being the solution of the direct problem for the assumed set of parameters and r(p (k) ) is the residual vector at the iteration k and N = M × F . In [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF], the authors used the Gauss-Newton method to solve the nonlinear least square problem which fails to converge in our case when the temperature is calculated at few sensors only but for the whole simulation time due to some lack of information. Moreover, Gauss-Newton method is not locally convergent on problems that are very non-linear or have very large residuals which is the case in our problem. Since the performance of the Gauss-Newton method is strongly dependent on the residual size, we adopted the use of the Damped Gauss Newton method which is an improved version of the Gauss-Newton algorithm [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF]. Damped Gauss-Newton method is known to be locally convergent on almost all nonlinear least squares problems including large residual or very nonlinear problems [START_REF] Björck | Numerical methods for least squares problems[END_REF].

The cost function S((ρC) s , φ, λ s ) dened by equation (4.10) can be re-written as:

S(p

(k) ) = 1 2 r T (p (k) )r(p (k) ) (4.12)
Such necessary condition for the minimization of S(p (k) ) can be represented in equation J(p k ) is dened by:

J(p (k) ) =                   W 1,(k) 1 R 1,(k) 1 Z 1,(k) 1 ... ... ... W F,(k) 1 R F,(k) 1 Z F,(k) 1 W 1,(k) 2 R 1,(k) 2 Z 1,(k) 2 ... ... ... W F,(k) 2 R F,(k) 2 Z F,(k) 2 ... ... ... W 1,(k) M R 1,(k) M Z 1,(k) M ... ... ... W F,(k) M R F,(k) M Z F,(k) M                   (4.14)
The elements of the sensitivity matrix are called the Sensitivity Coecients. The sensitivity coecient J f i,j is thus dened as the rst derivative of the estimated temper- ature at position i and time f with respect to the unknown parameter p j [ÖO00], that is,

J f i,j = ∂T f i ∂p j (4.15) where W f,(k) i = ∂T f i ∂(ρC) s | (ρC)s=(ρC) (k) s , R f,(k) i = ∂T f i ∂φ | φ=φ (k) and Z f,(k) i = ∂T f i ∂λ s | λs=λ (k) s .
The Damped Gauss Newton algorithm iteratively nds the minimum of S. Starting with an initial guess p (0) for the minimum, the method proceeds by the iterations:

p (k+1) = p (k) + m (k) (4.16)
m (k) is called the increment vector and is dened by:

m (k) = -d J(p (k) ) T J(p (k) ) -1 J(p (k) ) T r(p k ) (4.17)
d is the damping parameter (0 < d ≤ 1). An optimal value of d could be obtained using a line search algorithm [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF]; in our case, we used trial and error to nd a suitable constant damping parameter.

Governing Equations

In the following, we present the heat equation together with the three sensitivity equations resulting from the dierentiation of the heat diusion equation (4.1) with respect to the soil parameters p j . (In equation (4.18) the divergence operator and the partial derivative with respect to the parameters have been interchanged because the former depends on the spatial derivatives and the latter do not). where U j = ∂T /∂p j . The general sensitivity equation is accompanied with the following boundary and initial conditions:

∂ ∂p

At t = 0 : U j (x, 0) = U j0 = 0 in Ω At x = 0 : U j (x, t) = 0 for t ∈ (0, t end ] (Dirichlet) At x = l : ∂U j (x, t) ∂x = 0 for t ∈ (0, t end ] (Neumann)
The derivative of (ρC) e and λ e with respect to each soil parameter is given by:

d(ρC) e d(ρC) s = d(φ(ρC) f ) dT . dT d(ρC) s + (1 -φ) = φρ f dC f dT . dT d(ρC) s + φC f dρ f dT dT d(ρC) s + (1 -φ) = φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 . dT d(ρC) s + φC f (ρ v -ρ l ) dσ dT . dT d(ρC) s + (1 -φ) = φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 .W + φC f (ρ v -ρ l ) dσ dT .W + (1 -φ) (4.20) d(ρC) e dφ = d(φ(ρC) f dφ + d((1 -φ)(ρC) s ) dφ = (ρC) f + φ d((ρC) f ) dT dT dφ -(ρC) s = (ρC) f -(ρC) s + φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 Z + φC f (ρ v -ρ l ) dσ dT .Z (4.21) d(ρC) e dλ s = d(ρ C ) e dT dT dλ s = d(φ(ρC) f ) dT . dT dλ s = φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 R + φC f (ρ v -ρ l ) dσ dT .R (
(φλ s + (1 -φ)λ f ) 2 .W (4.23) dλ e dφ = E [φλ s + (1 -φ)λ f ] 2 (4.24) where E = [λ s dλ f dT dT dφ ][φλ s + (1 -φ)λ f ] -[λ s -λ f + (1 -φ) dλ f dT dT dφ ]λ f λ s dλ e dλ s = F [φλ s + (1 -φ)λ f ] 2 (4.25) where F = [λ s dλ f dT dT dλ s + λ f ][φλ s + (1 -φ)λ f ] -[φ + (1 -φ) dλ f dT dT dλ s ]λ f λ s and dλ f dT = (λ v -λ l ) dσ dT

Elimination of the approximation used in [MCM12]

In order to determine the sensitivity coecients (W, R and Z) appearing in the sensitivity matrix, we must solve the three sensitivity equations without using the approximation:

div (λ e ∇T ) ≈ λ e div (∇T )

(4.26)
used by Muhieddine et al. [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF] and which allowed the authors to write:

(ρC) e λ e ∂T ∂t = div (∇T )

(4.27)
This approximation leads to an approximate sensitivity matrix (Jacobian) and thus the problem will not converge to an exact solution ( maybe also due to lack of parameter's scaling mentioned before and studied later). Dierentiating the energy equation with respect to (ρC) s , φ and λ s respectively leads to:

∂W ∂t + 1 (ρC) e φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 W + φC f (ρ v -ρ l ) dσ dT W + (1 -φ) ∂T ∂t - 1 (ρC) e ∂ ∂x λ e ∂W ∂x - 1 (ρC) e ∂ ∂x φλ 2 s dλ f dT (φλ s + (1 -φ)λ f ) 2 W ∂T ∂x = 0 (4.28) ∂R ∂t + 1 (ρC) e (ρC) f -(ρC) s + φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 R + φC f (ρ v -ρ l ) dσ dT R ∂T ∂t - 1 (ρC) e ∂ ∂x λ e ∂R ∂x - 1 (ρC) e ∂ ∂x E (φλ s + (1 -φ)λ f ) 2 ∂T ∂x = 0 (4.29) ∂Z ∂t + 1 (ρC) e φρ f (C v -C l ) dσ dT + L d 2 σ dT 2 Z + φC f (ρ v -ρ l ) dσ dT Z ∂T ∂t - 1 (ρC) e ∂ ∂x λ e ∂Z ∂x - 1 (ρC) e ∂ ∂x F [φλ s + (1 -φ)λ f ] 2 ∂T ∂x = 0 (4.30)
These three sensitivity equations (4.28), (4.29) and (4.30) are completed with adequate initial and boundary conditions. W , R and Z are the unknowns of the sensitivity equations and T is the temperature. The obtained system of coupled equations (heat diusion equation + 3 sensitivity equations) is a nonlinear system of partial dierential equations. To solve this system, we use the same numerical strategy used in the forward problem (method of lines + nite volume method). After spatial discretization, the system of coupled equations can be written in the form:

F (t, Y, Y ) = 0 with Y (t 0 ) = Y 0 (4.31)
where Y = (T W R Z) t . The system in equation (4.31) can be solved by an ODE solver as in the forward problem.

Algorithm

The aim of the inverse problem is the calculation of the vector parameters p that minimizes the cost function S presented in equation (4.10). The Damped Gauss Newton (DGN) algorithm that we chose to apply to our nonlinear least square problem is as follows:

1-Choose a constant damping parameter d (0 < d ≤ 1).

2-Choose an initial value p (0) ; initialize the iteration k = 0 3-Perform the parameters' scaling to obtain p;

4-Solve the system (heat equation with phase change + sensitivity equations) using p(k) to dene the parameters of the soil. The equivalent parameters of the system are calculated by the apparent heat capacity method (AHC); 6-Solve the linear system J(p

• deduce T f,(k) i , W f,(k) i , Rf, (k)
(k) ) T J(p (k) )p (k+1) = J(p (k) ) T J(p (k) )p (k) -d J(p (k) ) t r (k)
for p(k+1) .

7-If the criteria of convergence are reached, end.

Calculate the original parameters' vector p (k+1) = D.p (k+1) .

If not, iterate:

p(k) ← p(k+1) and go to 4.

Stopping criteria

Classically, there are three convergence tests used in the algorithms for nonlinear least square problems (e.g [START_REF] Moré | User guide for MINPACK-1[END_REF]). We chose to apply only two of them. The rst test is the X-convergence which is based on an estimate of the distance between the current approximation x and the previous solution x* of the problem. If D is the current scaling matrix, then this convergence test attempts to guarantee that:

D(x -x * ) ≤ XT OL. Dx * (4.32)
where XTOL is a user supplied tolerance (we used XTOL = 10 -6 ). The second test, the main convergence test, is based on an estimate of the distance between the Euclidean norm F (x) of the residuals at the current approximation x and the previous value F (x * ) at the previous solution x* of the problem. This convergence test (Fconvergence) attempts to guarantee that:

F (x) ≤ (1 + F T OL). F (x * ) (4.33)
where FTOL is another user-supplied tolerance (we used FTOL = 10 -6 ).

Code verication

The code verication is based on choosing a plausible example where the soil parameters {(ρC) s , φ, λ s } are given constant values. These values are used by the forward problem to calculate the temperature at 5 dierent positions of the domain of length 10 cm (5 sensors at the depths x = 1 cm, x = 2 cm, x = 3 cm, x = 4 cm and x = 5 cm respectively) and are intended to show the eect of parameter scaling on the convergence of calculations. These temperatures are recorded every 24 seconds for 4 hours. In table 4.1, we used the same number of mesh cells in the forward problem (to create the synthetic data) and in the inverse problem. We also removed the approximation 4.26.

The results presented in this table prove that scaling is an important factor to obtain the desirable results. The target of our work is to perform a numerical simulation that is the closest possible to the real experimental case. For this reason, we generate the synthetic data using a very large number of mesh cells (around 6000) to obtain accurate results. So, these data play the role of the experimental data in the inverse problem which is run using small number of mesh cells (40, 80, 120, 160 ...). Figure 4.3 represents the convergence of the conductivity using 120 mesh cells in the inverse problem (the gures representing the convergence of the volumetric heat capacity and porosity are similar).

We notice that convergence is achieved after few tens of iterations (37 in this case).

Role of ∆T

As we have seen in chapter 3, the choice of the phase change temperature interval ∆T in the AHC method aects the temperature prole. Recall that ∆T is proportional to h [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] (∆T = kh) where h is the mesh size and k is a constant chosen in a way to obtain good accuracy with fewer uctuations in the temperature prole (we chose k = 1000). As a consequence, the value of ∆T plays an important role in the results of the inverse problem. If the initial values of parameters are far from the exact solution then the damped Gauss Newton method might not converge using the optimal value of ∆T (see subsection 3.2.2). To study the eect of this important parameter, we chose the example chosen earlier {(ρC) s = 1.95 × 10 6 , φ = 0.2, λ s = 0.756} and we run the inverse problem using dierent values of ∆T with 120 as number of mesh cells.

We notice that the inverse program fails to converge for ∆T kh = 1 or 2 and when ∆T kh ≥ 11 whereas it converges for 3 ≤ ∆T kh ≤ 10. We notice that the values of (ρC) J doesn't have full column rank and nally is that when the Gauss-Newton step is too long, the Levenberg Marquardt step is close to being in the steepest-descent direction -J T r and is often superior to the DGN step. We use the LMDER1 Minpack subroutine for numerical solution of nonlinear least square problems. LMDER1 is based on Moré's LMA version where the user must provide a subroutine to calculate the functions r 1 , r 2 , ... r m and the Jacobian matrix ∂r i (p) ∂p j

. LMDER1 follows the convergence criteria mentioned in section 4.2.5.

Algorithm

The aim of the inverse problem is the calculation of the vector parameters p that minimizes the cost function S presented in equation (4.10). The Levenberg Marquardt algorithm is summarized below:

1. Given ∆ (k) > 0, nd d (k) ≥ 0 such that if: (J (k) ) T J (k) + d (k) (D (k) ) T D (k) m (k) = -(J (k) ) T r (k) then either d (k) = 0 and D (k) .m (k) 2 ≤ ∆ (k) or d (k) > 0 and D (k) .m (k) 2 = ∆ (k) . 2. If r(p (k) + m (k) ) 2 < r(p (k) ) 2 set p (k+1) = p (k) + m (k)
and evaluate J (k+1) ; otherwise set p (k+1) = p (k) and J (k+1) = J (k) . 3. Choose ∆ (k) and D (k+1) .

The choice of ∆ and D is explained in details in [START_REF] Moré | The Levenberg-Marquardt algorithm: implementation and theory[END_REF].

Applying LMA to our Inverse Problem: Results

Using the LMDER1 Minpack subroutine (which is embedded in the easy-to-use MUESLI library [Canb]) and providing the Jacobian matrix, we obtain the results summarized in table 4.4. The Jacobian matrix is calculated using Maple. Also, using LMA, the inverse problem might not converge for a certain set of initial guess so the use of the technique of chaining the inverse problem is inevitable.

Identication of the thermophysical properties of a saturated soil during phase change in a 3D axisymmetric coordinate system

The inverse problem in a 1D coordinate system is not well-adapted to the experiments performed in the laboratory so it is more convenient to perform the inverse problem in a 3D axisymmetric coordinate system.

Forward problem (3D with axial symmetry)

The physical problem consists of heating the soil by a re. To model this problem, we replace the soil by a perfect porous medium in 3D axisymmetric domain of dimensions r = 30 cm and z = 50 cm with constant and uniform properties heated from above by a constant temperature T c (temperature of the re between 300 • C and 700 • C). T c must be greater than T v (the evaporation or phase change temperature which is normally 100 • C). with the following initial and boundary conditions:

At t 0 = 0,      T (r, z, t 0 ) = T c if z = 0 and 0 ≤ r ≤ l 1 2 T (r, z, t 0 ) = T 0 otherwise For t > t 0                            T (r, z, t) = T c if z = 0 and 0 ≤ r ≤ l 1 2 T (r, z, t) = T i if z = l 2 and r = l 1 ∂T ∂z = 0 if z = 0 and l 1 2 < r ≤ l 1 ∂T ∂r = 0 if r = 0
where T represents the temperature, T 0 is the initial temperature at t 0 = 0, T c is the re temperature; l 1 is the length of the domain along r axis , l 2 is the length of the domain along z axis ; ρ is the density, C is the specic heat capacity, λ is the conductivity, φ is the porosity, the subscripts e, f and s indicate the equivalent parameters of the medium, the properties of the uid and the porous matrix properties respectively. The computational domain is represented in gure 4.9 where we use Dirichlet and Neumann boundary conditions. In reality, all the boundary conditions except that of the re should be of Robin type because heat exchange always exists with the external medium (air). Some

Numerical method

As in 1D, we need to solve the heat diusion equation (PDE) so we choose the method of lines which is a way of approximating PDEs by ODEs where space and time discretizations are considered separately. The spatial discretization is done using the nite volume method, the equation is integrated over a control volume. We use N discretization points in the r direction and M in the z direction so that after discretization the heat equation becomes:

dT i,j dt - (r i + r i+1 )(λ i,j + λ i+1,j )(T i+1,j -T i,j ) 4(ρC) e r i ∆r 2 + (r i + r i-1 )(λ i,j + λ i-1,j )(T i,j -T i-1,j ) 4(ρC) e r i ∆r 2 - (λ i,j + λ i,j+1 )(T i,j+1 -T i,j ) 2(ρC) e ∆z 2 + (λ i,j + λ i,j-1 )(T i,j -T i,j-1 ) 2(ρC) e ∆z 2 = 0 (4.35)
where ∆r = r i+1 -r i and ∆z = z j+1 -z j . (see gure 4.10). In our conguration, we used ∆r = ∆z = h. As in 1 D case, we use an automatic ODE solver based on a BDF scheme which possesses the property of sti stability to overcome the stiness of AHC method at phase change temperature.

Inverse problem (3D with axial symmetry)

In order to solve the parametric inverse problem consisting of nding the volumetric heat capacity (ρC) s , the conductivity λ s and the porosity φ of the solid matrix of the saturated soil, it is necessary to know the values of temperature T f gi,j at selected points (sensors) of the porous medium domain for times t f : T f gi,j = T g (r i , z j , t f ) where i = 1, 2, ..., N , j = 1, 2, ..., M and f = 1, 2, ..., F . F is the total number of time steps.

We use the least squares criterion to solve this inverse problem so we try to nd the soil parameters that minimize the cost function which is dened by:

S((ρC) s , φ, λ s ) = 1 2 T f i,j -T f gi,j 2 2 (4.37)
where T f i,j = T (r i , z j , t f ) are the temperatures being the solution of the direct problem for the assumed set of parameters at the point (r i , z j ), i = 1, 2, ..., M , j = 1, 2, ..., N for the time t f , f = 1, 2, ..., F and T f gi,j is the measured temperature at the same point (r i , z j ) for time t f . chapter 4

Method of resolution

To illustrate the method of resolution, we dene the following vectors:

T g =                                                               T 1 g1,1 ... T 1 g1,M T 1 g2,1 ... T 1 g2,M ... T 1 gN,1 ... T 1 gN,M T 2 g1,1 ... T 2 g1,M T 2 g2,1 ... T 2 g2,M ... T 2 gN,1 ... T 2 gN,M ... T F g1,1 ... T F g1,M T F g2,1 ... T F g2,M ... T F gN,1 1,(k) 1,1 ... T 1,(k) 1,M T 1,(k) 2,1 ... T 1,(k) 2,M ... T 1,(k) N,1 ... T 1,(k) N,M T 2,(k) 1,1 ... T 2,(k) 1,M T 2,(k) 2,1 ... T 2,(k) 2,M ... T 2,(k) N,1 ... T 2,(k) N,M ... T F.(k) 1,1 ... T F.(k) 1,M T F.(k) 2,1 ... T F.(k) 2,M ... T F.(k) N,1 ... T F.(k) N,M                                                                     p (k) =    (ρC) (k) s λ (k) s φ (k)    and r(p (k) ) = g(p (k) ) -T g
where r(p (k) ) is the residual vector at the iteration k and N δ = N × M × F . We adopted the use of the Levenberg Marquardt Algorithm [DS83] which is known to be locally convergent on almost all nonlinear least squares problems including large residual or very nonlinear problems [START_REF] Björck | Numerical methods for least squares problems[END_REF].

To minimize the cost function, we use the same method of resolution as in 1D where the sensitivity matrix, J(p k ) is dened by:

J(p (k) ) =                                                                     W 1,(k) 1,1 R 1,(k) 1,1 Z 1,(k) 1,1 ... ... ... W 1,(k) 1,M R 1,(k) 1,M Z 1,(k) 1,M W 1,(k) 2,1 R 1,(k) 2,1 Z 1,(k) 2,1 ... ... ... W 1,(k) 2,M R 1,(k) 2,M Z 1,(k) 2,M ... ... ... W 1,(k) N,1 R 1,(k) N,1 Z 1,(k) N,1 ... ... ... W 1,(k) N,M R 1,(k) N,M Z 1,(k) N,M W 2,(k) 1,1 R 2,(k) 1,1 Z 2,(k) 1,1 ... ... ... W 2,(k) 1,M R 2,(k) 1,M Z 2,(k) 1,M W 2,(k) 2,1 R 2,(k) 2,1 Z 2,(k) 2,1 ... ... ... W 2,(k) 2,M R 2,(k) 2,M Z 2,(k) 2,M ... ... ... W 2,(k) N,1 R 2,(k) N,1 Z 2,(k) N,1 ... ... ... W 2,(k) N,M R 2,(k) N,M Z 2,(k) N,M ... ... ... W F,(k) 1,1 R F,(k) 1,1 Z F,(k) 1,1 ... ... ... W F,(k) 1,M R F,(k) 1,M Z F,(k) 1,M W F,(k) 2,1 R F,(k) 2,1 Z F,(k) 2,1 ... ... ... W F,(k) 2,M R F,(k) 2,M Z F,(k) 2,M ... ... ... W F,(k) N,1 R F,(k) N,1 Z F,(k) N,1 ... ... ... W F,(k) N,M R F,(k) N,M      (4.38)
The sensitivity coecient J f (i,j),δ is dened by:

J f (i,j),δ = ∂T f (i,j) ∂p δ (4.39)
where W f,(k)

(i,j) = ∂T f (i,j) ∂(ρC) s | (ρC)s=(ρC) (k) s , Z f,(k) (i,j) = ∂T f (i,j) ∂φ | φ=φ (k) and R f,(k) (i,j) = ∂T f (i,j) ∂λ s | λs=λ (k) s .
In theory, LMA is based on a trust region approach which is: nd m (k) that minimizes chapter 4

r(p (k) ) + J (k) m (k)
2 subject to m (k) 2 ≤ ∆ where ∆ is called the step bound [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF]. On the other side, in Moré's implementation it is m (k) that minimizes

r(p (k) ) + J (k) m (k) 2 subject to Dm (k) 2 ≤ ∆ (4.40)
where D is a diagonal matrix which takes into account the scaling of the problem. (see section 4.3.2 for the detailed algorithm). There are two convergence tests used in the algorithm similar to those explained in section 4.2.5.

Governing Equations

In the following, we present the heat equation together with the three sensitivity equations resulting from the dierentiation of the heat diusion equation (4.34) with respect to the soil parameters p δ where δ = 1, 2, 3. (ρC) e (T )

∂ ∂p δ (ρC) e ∂T ∂t = div ∂ ∂p δ [λ e ∇T ]
∂U δ (x, t) ∂t + d(ρC) e (T ) dp δ ∂T (r, z, t) ∂t = div (λ e (T )∇U δ (r, z, t)) + div dλ e (T ) dp δ ∇T (r, z, t) (4.42)
where U δ = ∂T /∂p δ . The general sensitivity equation is accompanied with the following boundary and initial conditions:

t = 0 : U δ (r, z, 0) = U δ 0 = 0 in Ω U δ (r, z, t) = U D δ (r, z, t) on Γ D × (0, t end ] (Dirichlet) U δ (r, z, t).ν = U N δ (r, z, t) on Γ N × (0, t end ] (Neumann)
Similar to the 1D case, the heat diusion equation must be dierentiated with respect to the soil parameters (ρC) s , φ and λ s respectively. The derivatives of (ρC) e and λ e with respect to each soil parameter are as given in section 4.2.3.

The sensitivity equations in the 3D axisymmetric coordinate system are as follows:

Dierentiating with respect to (ρC) s :

∂W ∂t + 1 (ρC) e d(

Code verication

The code verication is based on the same plausible example used in the 1D case. We record the temperature at 5 sensors of our domain whose dimensions are l 1 = 30 cm along r axis and l 2 = 50 cm along z axis . The sensors are located at the depths z = 1 cm, z = 2 cm, z = 3 cm, z = 4 cm and z = 5 cm respectively for 6 hours every 10 seconds.

In tables 4.5 and 4.6 below, we used the same number of mesh cells in the forward problem (to create the synthetic data) and in the inverse problem. The target of our work is to perform a numerical simulation that is the closest possible to the real experimental case. For this reason, we generate the synthetic data using a very large number of mesh cells (around 150000 = 300 × 500) then we run the inverse problem with small number of mesh cells. Figure 4.11 represents the variation of the nal residue as function of the number of mesh cells. We can easily notice that the residue decrease as number of mesh cells increase. Table 4.7 shows that values of the parameters for dierent mesh sizes knowing that the synthetic data is generated by 150000 mesh cells.

Comparison between 1D and 3D axisymmetric inverse problems

As we see in table 4.6, the inverse problem succeeds to converge in 3D axisymmetric coordinate even if the initial guess is far from the exact solution which was not the case in 1D. In table 4.8, we present some examples using dierent initial guesses where we notice that for certain initial guesses the 3D axisymmetric inverse problem converges while the 1D inverse problem fails to converge using LMA knowing that we used the same number of mesh cells in both the forward (to generate the synthetic data) and

Inuence of noised measures on the convergence of the inverse problem

In reality, the experimental measurement of temperatures will include some errors.

To be close to reality, we perturbed the synthetic values of temperature by adding a Gaussian noise of null average and standard deviation equal to 5 • C. Tables 4.9 and 4.10 represent the convergent values of the heat capacity, the thermal conductivity and the porosity of the soil, estimated by taking into consideration the noised measures in 1D and 3D axisymmetric coordinate systems respectively.

We can see clearly that the inverse problem converges in both 1D and 3D axisymmetric coordinate system but toward dierent converging values which are not very far away from the values obtained without taking the noise into consideration. Figure 4.12 reveals that the residue decrease as function of the iteration number when noise is taken into consideration. chapter 4

Sensitivity Analysis

In our model, the parameters have dierent order of magnitude so it is better to use scaled parameters. Also, since our model parameters ((ρC ) s , λ s , φ) and model output (temperature) have dierent units, it may be advantageous to compute relative sensitivities dened as:

J = ∂T ∂p p T (4.46)
In this study, we present a method (discussed in [PEZ + 09]) that uses singular value decomposition of the sensitivity matrix J followed by QR factorization. The sensitivity matrix in our model (see equation (4.38)) need to be reformulated to match with equation (4.46). After that, the singular value decomposition of J (J = USV T ) is used to obtain a numerical rank for J. This numerical rank is then used to determine ρ parameters that can be identied given the model output T . To estimate the number of uncorrelated parameters we used an error estimate in our computation of the Jacobian as a lower bound on acceptable singular values. For example, in the study analyzed here we used ODE solver with an absolute error tolerance of 10 -6 , i.e., the error of the numerical model solution is of order 10 -6 and the error in the Jacobian matrix is approximately √ 10 -6 = 10 -3 . Consequently, singular values should not be smaller than 10 -3 . Once the number of identiable parameters has been determined, we nd the most dominant parameters by performing a QR decomposition with column pivoting on the most dominant right singular vectors. Below we summarize subset selection method as an algorithm.

Subset selection algorithm:

1-Given an initial parameter estimate, p 0 , compute the Jacobian, J(p 0 ) and the singular value decomposition J = USV T , where S is a diagonal matrix containing the singular values of J in decreasing order, and V is an orthogonal matrix of right singular vectors.

2-Determine , the numerical rank of J. This can be done by determining a smallest allowable singular value.

3-Partition the matrix of eigenvectors in the form

V = [V ρ V n-ρ ].
4-Determine a permutation matrix P by constructing a QR decomposition with column pivoting, for V T ρ . That is, determine P such that: V T ρ P = QR; where Q is an orthogonal matrix and the rst ρ columns of R form an upper triangular matrix with diagonal elements in decreasing order. 5-Use the matrix P to re-order the parameter vector p according to p = P T p. The rst ρ elements of p are identiable. They are ordered from most sensitive to less sensitive.

Applying the algorithm described above on our inverse problem in 3D axisymmetric coordinate system using synthetic data using J as the nal sensitivity matrix before convergence, the singular values of J are σ 1 = 4.86×10 3 , σ 2 = 4.24×10 2 and σ 3 = 26.79.

chapter 4

This means that the 3 parameters are identiable. Applying a QR decomposition with column pivoting for V T we nd the permutation matrix P

P =   0 0 1 1 0 0 0 1 0  
From all what preceded, we deduce that φ is the most identiable parameter then λ s and nally (ρC) s . Figure 4.13 represents the singular value ratio for the 3 parameters.

Conclusion

The idea of this chapter was the enhancement of the inverse problem presented in [START_REF] Muhieddine | Heat transfer modeling in saturated porous media and identication of the thermophysical properties of the soil by inverse problem[END_REF]. In our work, we estimate the thermophysical properties of the soil knowing the temperature prole at selected points of the domain (few sensors). In order to solve this inverse problem, we used the least square criterion where we try to minimize the error function between synthetic measures and calculated ones. The coupled system composed of the energy equation together with the three sensitivity boundary initial value problems resulting from dierentiating the basic energy equation with respect to the soil properties must be solved. To overcome the stiness of our problem due to the high nonlinearity of the AHC method and the problem of large residuals we used the Damped Gauss Newton and Levenberg-Marquardt methods.

Our approach has the advantage that all soil parameters including the volumetric heat capacity converges to the required exact solution which was not the case in any previous work since we adopted a variable-scaling technique and we removed any approximation which leads to an exact sensitivity matrix (Jacobian) and thus reducing the number of iterations till convergence. Moreover, we emphasized on the importance of the choice of ∆T where for a certain initial guess the inverse problem fails to converge.

We overcome this problem by chaining the inverse problems using dierent values of ∆T . Finally, a 3D-axisymmetric version of the inverse problem was presented since it is more adapted to the experiments performed in the laboratory. chapter 4

Chapter 5

Estimation of the thermophysical properties of the soil by Inverse problem: application to experimental data

This chapter basically uses the method explained in the previous chapter to estimate the thermophysical properties of the soil using experimental data. Here, we tried to apply the 3D axisymmetric inverse problem to real experimental data (water saturated case) performed by Szubert [START_REF] Szubert | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] and Bergonzi [START_REF] Bergonzi | Étude de transferts de chaleur dans un foyer préhistorique expérimental[END_REF] in a joint work. Unfortunately, the problem failed to converge due to the fact that the physical model is not close enough to reality due to many reasons. As a consequence, a more complex model is needed which will be discussed in the last chapter.

Moreover, we used the same numerical strategy but this time to estimate the value of thermal diusivity α of the soil in case where no phase change exists (i.e. in a dry porous medium) using experimental data. The experiments were performed by J. A. Cordero [START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF] in the archeology lab of Rennes CReAAH. We were obliged to reformulate the inverse problem and add the sensors' positions as unknowns. This latter forced some kinds of constraints to be present in order to obtain a unique solution. Comparing the results obtained with those using other methods like Laloy & Massard [START_REF] Laloy | Nouvelle méthod thermique d'étude des foyers préhistoriques[END_REF] and Laplace transform [START_REF] Carslaw | Conduction of heat in solids[END_REF], we deduce that our inverse problem algorithm works quiet well.

5.1 Inverse problem in 3D axisymmetric coordinate system using experimental data

Experimental hearth: materials and design

Replicating the evaporation of water in a saturated soil (due to phase change) is not an easy procedure. To perform an experimental hearth, the soil is replaced by a perfect porous medium.

Materials used:

A stainless-steel box of size 50 × 50 × 30 cm.

Thermocouples of type K (Chromel/Alumel).

An acquisition card.

Labview software (the number of thermocouples, measurements time interval and a stop watch that displays the total duration of measures and their numbers were specied).

An electric heating plate.

A voltage converter for feeding the plate (imposed temperature).

Length measurement tools.

Fontainebleau sand (quasi-mono-disperse of diameter 150 to 200 µm).

Procedure: A certain amount of sand is weighted and then 22% water by mass is added. Everything is mixed in a bucket then poured into the experimental box. The thermocouples are placed on the vertical axis which passes through the center of the heating plate (i.e they are plunged horizontally 25 cm) at dierent depths such that at Inverse problem in 3D axisymmetric coordinate system using experimental data 133 least one of the thermocouples is in contact with the bottom of the heating plate which is used to control its temperature. The initial temperature T 0 is recorded then the plate is heated till it attains 400 • C and the software is set to record the temperature every 30 seconds for 6 hours. When the experiment is over, the sand is removed carefully to locate the position of the thermocouples. The distance with respect to the two perpendicular walls of the box (x and y coordinates) is measured together with the height from the upper edge in order to deduce the depth of the thermocouples. In this experiment, thermocouples were found at the depths: 0 cm, 1.4 cm, 2.3 cm, 3.2 cm, 6.2 cm and 10.6 cm. Note that the Cartesian coordinates are transformed into 3D-axisymmetric coordinates (r-z).

Inverse problem

To solve the inverse problem consisting of nding the thermophysical properties of the solid matrix of the saturated soil: heat capacity (ρC) s , the conductivity λ s and the porosity φ during phase change and using real experimental data in 3D axisymmetric coordinate system, we use the same numerical strategy explained in section 4.5. We recall that the numerical method used in section 4.5 proved to be ecient and accurate using synthetic data but unfortunately it is not the case with real experimental data. In fact, the physical parameters attain negative values during the iterative procedure of the LMA algorithm and as a consequence the algorithm fails to converge. To overcome this obstacle, we used LMA with bound constraints over parameters to force the parameters values to stay in a certain physical domain [START_REF] Shan | A Levenberg-Marquardt Method For Large-Scale Bound-Constrained Nonlinear Least-Squares[END_REF]. In other words, given an arbitrary point x, the projection of x onto the feasible bounded region is dened as follows. The i th component of the projection of x is given by:

P (x, l, u) i =    l i if x i < l i x i if x i ∈ [l i , u i ] u i if x i > u i (5.1)
where l i and u i represent the lower and upper bound of x i respectively.

Results

The results obtained using the experiment performed by Szubert and Bergonzi is summarized in the table 5.1: The bound constraints over the parameters used in this example are: (ρC) s ∈ [1.0 × 10 5 , 1.0 × 10 7 ], λ s ∈ [0.1, 9.0] and φ ∈ [0.1, 0.9]. We notice that (ρC) s converges chapter 5

to its lower bound whereas the exact physical value is of order 10 6 , λ s converges to 0.11 whereas the exact physical value is of order 0.3. Moreover, changing the box constrained used to: (ρC) s ∈ [1.0 × 10 3 , 1.0 × 10 7 ], λ s ∈ [0.1, 9.0] and φ ∈ [0.1, 0.9] gives totally dierent results which are summarized in the table 5.2. Comparing tables 5.1 and 5.2, we notice that the results are not homogeneous and thus our algorithm failed to converge in case experimental data is used. This failure is due to the fact that the physical model used during phase change is not close enough to reality because:

The convection term is neglected (No vapor ow in the porous medium but it could be added and tested).

The gravity and capillary forces are neglected but in reality the gravity:

forces water to ow and ll the dry zone and thus transforming the shape of the liquid/gas interface from a basin to a smooth trough;

will cause Rayleigh Bénard instability in some zones, which causes a uid ow due to a dierence in densities.

The temperature of the plate is not uniform (see gure 5.2).

Radiation is neglected although strong heating is used. It is known that Silica is semi-transparent to certain wavelengths. Therefore, although conductive transfer is dominant there is also a radiative transfer.

A boiling model (temperature > 100 degree Celsius) is not used although there is certainly a mixed zone constituted of vapor bubbles (unstationary process) and a sensor inside this zone will remain at 100 degree Celsius. The boiling zone thickness could be deduced from the time duration of the plateau of temperature curves.

Stemmelen et al. [START_REF] Stemmelen | Ébullition en milieu poreux: analyse expérimentale de la stabilité du front zone liquide-zone diphasique[END_REF] study the boiling of water initially present in a saturated porous medium. The chosen conguration (heating from below) inevitably leads, in some cases, to instability because gravity and buoyancy forces are in opposite directions.

The bottom is waterproof, the vapor created must raise the liquid mass, and/or pass through it to exit. We then understood that the thickness of the two-phase zone is important.

In our case, and although the problem is very close, heating is carried out from above.

The steam created at the interface comes out naturally (open edge, while the bottom is closed) and it is conceivable that if boiling water exist, then the two-phase zone has a much smaller thickness than that in the cited study. As a consequence, a more complex model is needed in order to approximate the thermophysical properties of the soil during phase change but this is out of the scope of this work. On the other hand, in the next part of this chapter we will study in details another inverse problem to approximate the thermophysical properties in case of a dry porous medium where no phase change is present.

5.2 Identication of the thermophysical properties of the soil in a dry porous medium.

The soil is treated as a continuous medium subject to a purely diusive transfer. This last hypothesis is a priori questionable because the heating is intense (radiation is important) and the temperature gradients can be signicant (both vertical and horizontal components) convective transfer exists. The energy balance is simply reected as we have seen earlier by the heat equation:

(ρC) e ∂T ∂t = div (λ e ∇T ) (5.2)
Where T is the temperature, (ρC) e is the heat capacity per unit volume and λ e is the thermal conductivity of the medium.

chapter 5

If the medium is homogeneous and isotropic, and although the properties of air vary with the temperature, we neglect the variation of λ e with temperature to obtain:

∂T ∂t = α div ( ∇T ) (5.3)
Where α = λe (ρC)e is the thermal diusivity. We see that through all these assumptions (homogeneous eective medium, transfer by conduction only with conductivity independent of temperature T ), the transfer is completely dened by the geometry of the medium, its diusivity, boundary conditions and initial conditions. For a dry medium, which constitute of a solid phase of volumetric heat capacity (ρC) s and thermal conductivity λ s and uid phase (dry air) of volumetric heat capacity (ρC) f and thermal conductivity λ f , the eective volumetric heat capacity of the medium can be easily calculated by a simple linear formulation:

(ρC) e = φ(ρC) f + (1 -φ)(ρC) s (5.4)
Where φ represents the porosity of the medium. The volumetric heat capacity of the medium can be easily calculated using the previous equation. On the contrary, the expression of the eective thermal conductivity as function of that of solid and uid is unknown. Approximating the diusivity by an inverse problem similar to that used previously will be a good way to calculate λ e .

Experimental hearth: materials and design

To perform an experimental hearth, replace the dry soil by a perfect dry porous medium. The materials used are exactly the same described in the experiment done in the case of phase change but we used the small heating plate (radius = 5 cm) instead of the big one.

Procedure:

Pour a certain amount of dry sand into the experimental box. The thermocouples are placed at dierent depths such that at least one of the thermocouples is in contact with the bottom of the heating plate which is used to control its temperature. Record the initial temperature T 0 then we heat the plate till it attains 600 • C and we set the software to record the temperature every 30 seconds for 6 hours. When the experiment is over, we delve carefully to locate the position of the thermocouples. Measure the distance with respect to the two perpendicular walls of the box (x and y coordinates)

and the height from the upper edge in order to deduce the depth of the thermocouples.

In this study we explore two dierent experiments. Note that, porosity was estimated from the density of the sand (Fontainebleau sand) and the density of the solid alone.

The value of the latter (2540 kg/m 3 ) was determined from a measurement in water (measure of volume of solid for a given mass of sand). Porosity φ represents the volume of pores to that of the total volume.

Methods used to nd an approximate value of diusivity α

In order to validate the inverse problem method used to approximate the value of diusivity of a dry soil, we calculate the value of diusivity of the same set of data using two other dierent methods. These two (not fully numerical) methods are: Laloy & Massard and Laplace transform.

Laloy and Massard (Done by ARPHYMAT team)

In their paper [START_REF] Laloy | Nouvelle méthod thermique d'étude des foyers préhistoriques[END_REF], Laloy and Massard presented a way to estimate the duration of re. They assumed that the medium is homogeneous and isotropic and that physical characteristics such as conductivity, heat capacity and diusivity are independent of temperature and of spatial coordinates. In addition, they assumed that the dimension of re is larger than the height of the altered soil. These hypotheses allow a simple study of at re in a semi-innite domain and hence solving analytically the heat diusion equation in the case of a 1D geometry (function of the depth z).

The problem consists of nding T = f (z, t) which is the solution of equation ( 5.3) where div ( ∇T ) = ∂ 2 T ∂z 2 with α = constant. If T 0 is the temperature of re (at t = 0 and z = 0), T i is the initial temperature of the domain then:

T 0 -T (z, t) T 0 -T i = 2 √ π x 0 e -u 2 du where x = z 2 √ αt (5.5)
Knowing that:

e -2x 2 ≤ 1 - 4 π x 0 e -u 2 du 2 ≤ e -x 2
Laloy and Massard approximated the expression 1 - To determine α, according to Laloy and Massard, we simply need to know the variation of temperature in z-position at each instant t (T (z, t)) and then plot, for dierent values of z, z 2 as function of -ln 1 -T 0 -T (z,t) T 0 -T i

2

. If the assumptions are true, for each instant t, we should get a straight line passing through the origin, of slope ( 8 3 αt). Normally only one of these lines is enough to obtain α (if t is known). To improve accuracy and test hypotheses, we must determine, at dierent times, the slope a = This numerical calculation is done by E. Canot for the transfer of heat in a dry soil using 500 × 600 mesh cells. The total time of the simulation is 12 h.

Important Note:

Laloy and Massard method is valid in 1D only and is accurate if the temperature is measured at a depth close to the re (check gure 5.5 to see the dierence between the 1D and 3D-axisymmetric cases).

Results

First Experiment

Second Experiment Diusivity α in m 2 /s 2.5 × 10 -7 2.34 × 10 -7

Table 5.3: The values of diusivity (Fontainebleau Sand) obtained by Laloy and Massard method. Figure 5.4 represents the second experiment.

The two experiments are done using small heating plate of radius 5 cm and a box of dimensions 50 cm × 50 cm × 30 cm so the computational domain is in 3D axisymmetric coordinate system. Therefore, Laloy and Massard method will not be accurate due to the fact that we lost the 1D property and that the sensors are not all close to the re and to the axis respectively. As we noticed in gure 5.5, the error between 1D and 3D axisymmetric congurations is reduced close to the re and next to the axisymmetric chapter 5

axis whereas it increases away from the axis and as the depth increases. For all what proceeded, Laloy and Massard method is not accurate enough for an experiment done in 3D axisymmetric coordinate system.

5.2.2.2 Laplace Transform (Done by R. Delannay and É. Canot) In this part, we solve the heat diusion equation in 1D using Laplace transform in time u(x, t) → U (x, p) = L {u(x, t)} (p).

∂T ∂t = α ∂ 2 T ∂x 2
(5.7)

with T (x, 0) = 0 and T (0, t) = T 0 (t). Using Laplace transform, the heat equation (5.7) becomes:

p T (x, p) = α ∂ 2 T (x, p) ∂x 2

(5.8)

The solution of equation (5.8) is given by : T (x, p) = A(p) exp -p α x (see [START_REF] Carslaw | Conduction of heat in solids[END_REF]), where A(p) comes from the boundary excitation T 0 (t).

Methodology

Assume that we have two two sensors at two dierent depths z 1 and z 2 then we can say that T1 (p) and T2 (p) at z 1 and z 2 respectively are given by: T1 (p) = A(p) exp -p α z 1

(5.9)

T2 (p) = A(p) exp - p α z 2
(5.10) then dividing equations 5.10 and 5.9 we get:

T2 (p) T1 (p) = r(p) = exp - p α (z 2 -z 1 ) (5.11)
we end up with:

(∆z) 2 p = α [log (r(p))] 2 where ∆z = z 2 -z 1

(5.12)

Simply, to approximate α, we plot the graph of (∆z) 2 p as function of [log (r(p))] 2 and we calculate its slope.

Results

As we mentioned earlier, we need a couple of sensors only to approximate α. To ensure the validity of our method we use synthetic data where α = 1.0 × 10 -7 m 2 /s and we try to recover α using Laplace transform. The computational domain is in 3D axisymmetric coordinate system where r = z = 45 cm and the re radius is R = 15 cm. Figure 5.6

shows the variation of the Laplace of the temperature at two dierent sensors as function of p, gure 5.7 is a check for the validity of the ratio expressed in equation (5.11), gure 5.8 allows us to calculate the slope which eventually equal to the diusivity α. We can see clearly that the error between the exact and the calculated solution is 0.4%. Laplace transform method is valid for 1D conguration, so we need to check its precision in 3D axisymmetric coordinate system. We try to detect the variation of error obtained upon the choice of dierent sensors' positions. First, we study the eect of the sensor's depth on the accuracy for results. In gure 5.9, we took into consideration 8 sensors located on the axisymmetric axis at dierent depths such that the ratio of depth to that of re radius varies from 0.033 to 1.333 (0.033 ≤ z R ≤ 1.33). In gure 5.10, we notice that the error between the dierent calculated values and the exact one (1.0 × 10 -7 ) is between 0.4% and 0.44%.

In gure 5.11, we took into consideration 8 sensors that share the same property (r = z) but their positions vary from r = z = 4cm to r = z = 18cm, the values of α obtained using various sensors combinations are shown in gure 5.12 where we can see clearly that as the sensor's depth and distance from axis increase the error increases where it varies from 0.27% to 33.58%. Looking at the results obtained by synthetic data, we conclude that to obtain accurate result while using Laplace Transform method in 3D axisymmetric coordinate system, we must choose the sensors close to the axis to avoid large errors. Table 5.4 and table 5.5 represent dierent values of α obtained using dierent couples of sensors using real experimental data precisely experiments number 1 and 2. In the rst experiment, we choose 3 couples of sensors that are not far from the axis where 0.32 ≤ r R ≤ 0.46 and the same time not far from the re where 0.24 ≤ z R ≤ 0.76. For the second experiment, we also choose 3 couples of sensors that are not far from the axis where 0.28 ≤ r R ≤ 0.42 and the same time not far from the re where 0.3 ≤ z R ≤ 0.6. Diusivity α in m 2 /s 3.28 × 10 -7 2.04 × 10 -7 2.02 × 10 -7

Table 5.5: The values of diusivity (Fontainebleau Sand) obtained by Laplace transform method using the second experiment.

The average value of α obtained using the rst experiment is 2.63 × 10 -7 m 2 /s and that using the second experiment is 2.44 × 10 -7 .

Inverse Problem in 3D-axisymmetric coordinate system

Objective Similar to the inverse problem visited earlier, we need to estimate the thermal diusivity α of the soil by inverse problem knowing the history curves at selected few sensors of the domain. To do this, we use the least square criterion where we try to minimize the error function S which represents the dierence between the experimental temperature and the numerical temperature (see section 4.3.3):

S(α) = 1 2 M i=1 F j=1 T f i,num -T f i,exp 2
(5.13)
where T f i,num = T (x i , t f ) are the temperatures being the solution of the direct problem for the assumed parameter(α) at the point x i , i = 1, 2, ..., M for the time t f , f = 1, 2, ..., F and T f i,exp is the measured experimental temperature at the same point x i for time t f .

Forward Problem

The physical problem consists of heating the dry soil by a re (no phase change). To model this problem, we replace the soil by a perfect porous medium with constant and chapter 5 uniform properties heated from above by a constant temperature T c . In this case, we assume that the thermal conductivity λ e is independent of temperature and that the heat diusion equation simplify to:

∂T ∂t = α 1 r ∂ ∂r r ∂T ∂r + ∂ ∂z ∂T ∂z (5.14)
with the following initial and boundary conditions:

t = 0 : T (r, z, 0) = T 0 in Ω T (r, z, t) = T D (r, z, t) on Γ D × (0, t end ] (Dirichlet)
∇T.ν = 0 on Γ N × (0, t end ] (Neumann)

where T represents the temperature, T 0 is the initial temperature at t 0 = 0, T D is T c at the re and T 0 elsewhere; ν indicates the outward unit normal vector along the boundary of Ω. The computational domain is similar to that presented in gure 4.9.

Method of Resolution

The method of resolution is similar to that presented earlier in section 4.5.1 where

p (k) = (α (k)
) and the Jacobian(sensitivity) matrix ended up to be a vector:

      U 1,(k) α1,1 ... U 1,(k) α,M U 1,(k) α2,1 ... U 1,(k) α2,M ... U 1,(k) αN,1 ... U 1,(k) αN,M U 2,(k) α1,1 ... U 2,(k) α1,M U 2,(k) α2,1 ... U 2,(k) α2,M ... U 2,(k) αN,1 ... U 2,(k) αN,M ... U F,(k) α1,1 ... U F,(k) α1,M U F,(k) α2,1 ... U F,(k) α2,M ... U F,(k) αN,1 ... U F,(k) αN,M      (5.15) where U f,(k) α(i,j) = ∂T f (i,j) ∂α | α=α (k)
. We try to nd α by minimizing S(α) using LMA.

Governing Equations and Numerical Strategy

As in 4.5.2, we need to dierentiate the heat diusion equation in dry case (5.14) with respect to the unknown parameter α:

∂ ∂α ∂T ∂t = ∂ ∂α (α div [ ∇T ])
(5.16) chapter 5

which leads to the α sensitivity equation below:

∂U α (x, t) ∂t = div (∇T (x, t)) + α div (∇U α )

(5.17)

The initial and the boundary conditions for the sensitivity equation are similar to those in 4.5.2. The obtained system of coupled equations (heat diusion equation and the sensitivity equation with respect to α) forms a system of partial dierential equations. This system is solved using the same strategy explained in 4.1.1.

Results

We test our inverse problem using real experimental data because using synthetic data we are pretty sure that the inverse problem will provide accurate results. As we have seen in previous sections, using Laloy/Massard and Laplace methods, the value of diusivity is of order 10 -7 . Also, examining the literature, the value of α for Fontainebleau sand is always found to be of order 10 -7 . In [START_REF] Combarieu | Caractérisation mécanique d'un massif de sable compacté: Cohérence des essais réalisés[END_REF], it is reported that the density of Fontainebleau sand is 1480 kg/m 3 and in other references, we can nd that its volumetric heat capacity is 1.041 × 10 6 J/kgK and its thermal conductivity is 0.32 W/mK so we can simply deduce that 3 × 10 -7 m 2 /s. Moreover, Narsilio et al. [START_REF] Narsilio | Hydraulic and thermal conduction phenomena in soils at the particle-scale: Towards realistic fem simulations[END_REF] indicates that the thermal conductivity of Fontainebleau sand ranges between 0.28 and 0.42 W/mK and hence its diusivity will be of order 10 -7 . Comparing the results obtained in table 5.6 and what we mentioned now, we can see that the error is of order 10 2 which is huge.

We suspect that the error obtained is due to various experimental errors:

1. The thermocouples are exible and long (50 cm). During lling the box with sand, perhaps the thermocouple rod is curved down and upon removing the sand, to measure its position, it goes back to its straight shape.

2. The position of the thermocouple is located by measuring the distance with respect to two perpendicular walls of the box (coordinates x and y) and the height to the upper edge thereof, with a string (having previously identied, in the same manner, the height of the sand surface, we deduce the depth of the thermocouples).

The sensors' coordinates measured in cartesian coordinates are transformed by hand into 3D axisymmetric coordinates. Perhaps this process will encounter some errors.

3. Observing the measured initial temperatures, we notice that it is not uniform with depth.

4. The circular heating plate used in the experiment do not have a uniform temperature (see gure 5.2).

To treat the error resulted from measuring the positions of sensors we propose:

(i) Adding the sensors' positions as unknowns in addition to diusivity.

(ii) Use a curve t to treat the problem of having non-uniform initial temperature.

Estimation of α and sensors' positions

The analytical solution of the heat equation in 1D semi-innite domain [START_REF] Carslaw | Conduction of heat in solids[END_REF] is of the form:

T (z, t) = k erf z √ αt (5.18)
where T is the temperature, k is a constant, z is the position, α is the diusivity and t is the time. Based on equation 5.18, it is not possible to nd both the diusivity of a solid medium and the position of the sensors, because there are an innite number of solutions. This arise not only in 1D problems, but also in 2D or 3D congurations. (i) It is known that the analytical solution of the steady state heat equation in 3D

axisymmetric coordinate system is a family of ellipses.

(ii) Numerical tests showed that the solution of the transient heat equation in 3D axisymmetric coordinate system is very close to a family of ellipses.

As a result, we introduce a constraint by imposing the displacement normal to these isothermal lines. This writes:

sin(φ i ) δ r i = cos(φ i ) δ z i (5.19)
for all sensors i. In equation 5.19, φ i denotes the angle of the isothermal lines at sensor number i. δ r i (resp. δ z i ) is the (unknown) displacement of the sensor i in the r direction (resp. z direction).

Furthermore, we add new constraints about the unknown position of sensors, which specify that their mean displacement is zero, both in r and z direction. This comes from the assumption that the position errors obey a normal centered distribution law and it is generally true for a great number of sensors. This writes:

ri = r i + δ r i + shif t r zi = z i + δ z i + shif t z (5.21)
for all sensors i. shif t r and shif t z represent the bias and are also considered as unknowns, but global to our problem.

Results

As we mentioned earlier, we are going to apply our inverse problem to two experiments.

First Experiment (experiment number 13 in [START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF]): The thermocouples in this experiment are presented in gure 5.16, the re radius is R = 5 cm, the dimensions of the computational domain are r = 25 cm and z = 30 cm. We can see clearly that the sensors are close to the z-axis and their distance from the axisymmetric axis with respect to re radius is (0.18 ≤ r R ≤ 0.54) where as their depths with respect to the re radius is (0.18 ≤ z R ≤ 1.96). Figure 5.17 shows the initial temperature prole for the rst experiment: we can check that it is nearly uniform in space (the initial temperature varies from 24 to 26.5 • C). Using a numerical mesh of 300 × 500, we obtain a diusivity α = 1.73 × 10 -7 m 2 /s which is acceptable for Fontainebleau sand. The standard deviation for both δ r i and δ z i are respectively 0.8 mm and 5.5 mm, showing that the z position of the sensor is more dicult to be obtained accurately.

Besides, the shifts in r and z directions are respectively 1.4 mm and 9.8 mm; the last value appears to be somewhat too large to be attributed to an experimental protocol error. Therefore, the inverse problem should be investigated further on. However, we are condent of the numerical results, because the convergence curves during the iteration process are good (see gures 5.18 and 5.19). The mesh size varies from 30 × 50 to 300 × 500. The value of residue converges to 106.4 Second Experiment (experiment number 12 in [START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF]): The thermocouples in this experiment are presented in gure 5.21, the re radius is R = 5 cm, the dimensions of the computational domain are r = 25 cm and z = 30 cm. We can notice that the sensors are spread in the domain unlike the rst experiment where they were approximately close to the axisymmetric axis (0.28 ≤ r R ≤ 1.44) and (0.04 ≤ z R ≤ 1.9). Figure 5.22 shows that the initial temperature for the second experiment is not uniform, it varies from 25.5 to 32 • C. In this gure, we modeled the initial temperature by constant mean of temperatures (dotted line) and by an exponential law (continuous line). This kind of initial prole reveals that the soil was not in equilibrium and that the experiment occurred during a slow transient heating of the medium due to an increase in the room temperature. Nevertheless, the t of the temperature data to an exponential law works quite well.

Results: uniform initial temperature Using a 300 × 500 mesh, we used the inverse problem to nd the value of diusivity, δ r i , δ z i , shif t r and shif t z (i.e. the new sensors' positions). The values obtained using uniform (constant) and non-uniform (exponential t) initial temperature are shown in

Conclusion

The inverse problem algorithm failed to give good results using real experimental data during evaporation in a saturated porous medium. We render this failure to many reasons: the physical model used during phase change is not close enough to reality and the fact that the experiments are not well-controlled. Actually, the thermocouples positions are not exact, the temperature distribution on the heating plate is not uniform and some uncertainties concerning the exact initial state of the sand (moisture distribution, packing fraction, ...).

Concerning the inverse problem in dry case, we saw that the use of statistical and box constraints were unavoidable since the thermocouples position is dicult to be measured accurately due to the experimental setup. Moreover, each experiment has its own conditions and thus data should be treated accordingly.

Although the Fontainebleau sand has been used in both analyzed experiments, different values for the thermal diusivity were obtained (α = 1.73 × 10 -7 m 2 /s and α = 2.42 × 10 -7 m 2 /s). We think that this discrepancy comes from the fact that it is dicult to have well controlled experiments. However, it is worth noting that the magnitude order match quite well with our results obtained with other methods (Laloy & Massard, and the Laplace transform) and also with literature data. T < 100 • C then it is liquid water and otherwise it is water vapor). This method is well-adapted to industrial problems concerning melting/solidication of metals [START_REF] Hu | Mathematical modeling of solidication and melting: a review[END_REF][START_REF] Poirier | On numerical methods used in mathematical modeling of phase change in liquid metals[END_REF].

In case of phase change in a porous medium, the AHC method presented in [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] and [START_REF] Muhieddine | Coupling heat conduction and water-steam ow in a saturated porous medium[END_REF] considers the heating of saturated soil where the phase change occurs when liquid water turns into gas (i.e. authors studied evaporation only). Actually, archaeologists are interested in studying not only the heating phase of prehistoric hearths but also the cooling stage because the ancient hearths are made up of a sequence of heating and cooling stages.

In fact, the starting point of the models, discussed in [START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] and [START_REF] Muhieddine | Coupling heat conduction and water-steam ow in a saturated porous medium[END_REF], is a fully saturated porous medium which after heating becomes a medium formed of two separated phases one contains water only and the other water steam only. If after this heating stage, a cooling stage is started then we will numerically end up with a domain fully saturated with water which do not models properly the real situation. In reality, during the cooling stage of a hearth, air from surrounding enters the domain and diffuses inside the water vapor which condensed into liquid water so the old model which is based on the AHC method cannot be used to describe the cooling stage because of absence of air; this is quite dierent.

In this chapter, we present two physical models. The rst one aims to palliate the deciency of our current physical model (presented in previous chapters). The new constituent air, which can enter the soil from the atmosphere, is added to the model.

Beside that, in many situations, soils are not saturated but close to a dry granular medium with a few quantity of liquid water. This is the case of the experiments of dry sand used in chapter 5 and for which we suspect the presence of a small quantity of water.

chapter 6

During cold desert nights, water present in the atmosphere condenses transforming the dry sand into a humid granular medium. One of the most important application is the stability of sand dunes [LVeM + 13]. Our second model presents the case of a granular medium in pendular regime (small quantity of liquid water) where the results related to the calculation of the eective thermal conductivity for a wet granular medium (obtained in chapter 2) can be used. This second model is proposed without numerical results.

Phase change in a saturated porous medium in presence of air

The rst model that will be presented in this chapter takes into account a saturated porous medium in 1D coordinate system in presence of air which is a non-condensable gas. The medium now is formed of four components: the solid porous matrix, liquid water, water vapor and air. The model is therefore a separated phase model where the rst phase is composed of liquid water only while the second contains only gases (water vapor + air). We treat the physical problem by coupling the heat conduction in a saturated porous medium with gas ow taking in consideration the diusion of air inside water vapor.

The easiest manner to represent this model is to slightly modify the previous model (water saturated model) by restricting the presence of air in the gaseous part only and by considering it as water-insoluble. Our new model, as announced above, has two distinct zones: a dry one that contains gas only and another one that remains saturated (where temperature doesn't exceed 100 • C). The dried zone, in contact with the atmosphere, and which is the seat of diusion between air and water vapor, is found in many models that treat soil evaporation. But almost all of these models has an unsaturated zone, which is the seat of complex phenomena related to capillary forces and gravity.

The models that deal with the evaporation of the surface of the soil are numerous and so important that there is even a journal that deals just with this phenomenon Vadose Journal (Vadose is the Latin word for the supercial part of the soil). In reality, what happens in this gaseous region is coupled to both the top with the atmosphere (this latter may contain more or less moisture, be animated by a certain speed, etc.) and with the bottom with the phreatic nap that feeds the unsaturated zone with water by imbibition [START_REF] Gowing | Evaporative ux from a shallow watertable: The inuence of a vapour-liquid phase transition[END_REF].

A model close to ours is that of Gardner [START_REF] Gardner | Some steady-state solutions of the unsaturated moisture ow equation with application to evaporation from a water table[END_REF], which considers both zones, the rst is a deep saturated zone; the other is the supercial zone, where the gas ows freely but which may contain liquid water trapped in certain areas.

Another class of important applications is related to the drying of materials, whether in the eld of civil engineering where fresh concrete walls need to dry thoroughly [START_REF] Pel | Moisture transport in porous building materials[END_REF][START_REF] Daian | Condensation and isothermal water transfer in cement mortar. part i-pore size distribution, equilibrium water condensation and imbibition[END_REF][START_REF] Daian | Condensation and isothermal water transfer in cement mortar. part ii-transient condensation of water vapor[END_REF] or in the eld of wood industry, where the wood must be dried before use [START_REF] Perré | Multiscale aspects of heat and mass transfer during drying[END_REF]. Most of these studies, if not all, treat the wet porous medium as unsaturated.

Our model is simple but it constitutes a supplementary step towards a more realistic model (including drainage and imbibition).

Physical description and mathematical modeling

Consider a horizontal surface of a saturated porous medium in presence of air. To describe the evaporation of water in a soil and to simplify the model, we suppose that we have a separated phase model as shown in gure 6.1.

Water saturated zone Diffusion of air inside vapor

Gas (vapor + air)

Ground

Fire

Pure Conduction Water Vapor Flow

Liquid/Vapor interface Contrary to the assumptions used in the AHC method where the computational domain is considered as one region, here, we separate the liquid zone from the gaseous one.

Moreover, the interface position is calculated explicitly whereas in the AHC method the phase change interface position is embedded in the method.

Governing Equations

Many main variables can be chosen for the computation; in this part, we chose: the temperature (T ), the pressure (P ) and the molar density of water vapor (n w ).

1. The energy conservation in the two zones is expressed as:

(ρC) e ∂T ∂t + (ρC) f v f .∇T = div(λ e ∇T ).

(6.1) Using Darcy's law, which is expressed usually by:

v f = - K µ f ∇P (6.2)
the energy conservation becomes: chapter 6

(ρC) e ∂T ∂t - K(ρC) f µ f ∇P .∇T = div(λ e ∇T ). (6.3)
where (ρC) e is the eective heat capacity which is dened by the following equation:

(ρC) e = φρ f C f + (1 -φ)ρ s C s . (6.4)
φ is the porosity, K is the permeability, µ f is the viscosity of the uid, ρ f is the density of the uid, C f is the capacity of the uid, T is the temperature, v f is the ltration velocity of Darcy for the uid phase, P is the pressure of the uid, ρ s and C s are respectively the density and the capacity of the solid matrix and λ e is the eective conductivity calculated using Kunii and Smith model [MHB + 12], β = 0.895 and Φ = 0.058 (see equation 6.5).

λ e = λ f φ + (1 -φ) β λ f λs + Φ (6.5)
2. The mass conservation of the uid in the gaseous zone is given by:

∂(φρ f ) ∂t + div(ρ f v f ) = 0 (6.6)
Using Darcy's law, the mass conservation becomes:

∂(φρ f ) ∂t + div(- ρ f K µ f ∇P ) = 0 (6.7)
We assume that the gases obey the ideal gas law and thus:

P ρ f T = R M = β (6.8)
where R is the ideal gas constant and M is the molar mass of the gas mixture.

The mass conservation can be written now in the following form:

φ 1 βT ∂P ∂t - P βT 2 ∂T ∂t + div(- ρ f K µ f ∇P ) = 0 (6.9)
Notice that, since we are in the gaseous phase, the subscript f can be replaced by the subscript g that stands for the gas mixture.

3. The conservation of mass for water vapor in the gaseous phase is given by: ∂(φn w ) ∂t + div(n w v w ) = 0 (6.10)

where n w is the molar density of water vapor and v w is the average velocity of water vapor which can be written in the form:

v w = v g + v * w (6.11)
where v g is the barycentric (mass average) speed of the gas (air and water vapor) and v * w is the relative velocity of water vapor. This velocity decomposition is classical in diusion models in order to easily write Fick's law. Substituting 6.11 in 6.10, the mass conservation of water vapor becomes:

∂(φn w ) ∂t + div(n w v g + n w v * w ) = 0 (6.12)
The molar ux of water vapor J * w is given by: J * w = n w .v * w so using the linear Fick's law which is an approximation of the true diusion process, we obtain:

J * w ≈ -D w,a ∇n w (6.13)
the mass conservation of water vapor ends up as follows:

∂(φn w ) ∂t + div(n w v g -D w,a ∇n w ) = 0 (6.14)
where D w,a is the coecient of diusion between air and water vapor. Here, we use a temperature dependent law for D w,a taken from [START_REF] Nellis | Heat Transfer[END_REF]. Anyway, the diusion coecient in a porous medium is slower than that in a free space; for this reason equation 6.13 should include a tortuosity factor 0 < τ < 1 (see [START_REF] Daian | Condensation and isothermal water transfer in cement mortar. part i-pore size distribution, equilibrium water condensation and imbibition[END_REF]). In this work, this fact is neglected.

4. The equations at the interface ξ(t) are:

(i) the energy balance at the interface:

λ w ∂T ∂x ξ + -λ l ∂T ∂x ξ - = -φ ρ l L dξ dt (6.15)
where λ l is the conductivity of liquid water, λ w is the conductivity of water vapor and L is the latent heat of evaporation.

(ii) to write the mass balance equation at the interface, we rst write the mass balance equation for water only (liquid water + water vapor):

ρ l v l - dξ dt ξ - = ρ w v w - dξ dt ξ + (6.16)
but the liquid water in our model is assumed immobile (due to the boundary condition of pressure at the bottom of the domain, it will be written later), hence v l = 0 and the above equation reduces to:

ρ w v w = - dξ dt (ρ l -ρ w ) (6.17)
chapter 6

We know that ρ w << ρ l so the mass balance equation at the interface for water only is given by:

ρ w v w ≈ - dξ dt ρ l (6.18)
We actually have a mixture of two constituents water (liquid water + water vapor) and air. Air and water are immiscible so the velocity of air is equal to that of the interface i.e. v a = dξ dt . Taking into consideration that v l = 0 and ρ g v g = ρ w v w + ρ a v a , we can write:

ρ g v g = -ρ l dξ dt + ρ a dξ dt (6.19)
but ρ a << ρ l so we neglect the second term and end up with:

ρ g v g ≈ -ρ l dξ dt (6.20)
(iii) to track the interface position, we need an equation that relates its position with one of the main variables present in the system. Starting with equation 6.18 and taking into account equation 6.11 and Fick's law, we arrive at the following equation:

ρ w v w = ρ w v g + ρ w v * w -ρ l dξ dt = n w M w v g + n w M w v * w -ρ l dξ dt = -n w M w K µ g ∇P + M w J * w ρ l dξ dt = n w M w K µ g ∇P + M w D w,a ∇n w (6.21)
The system of equations to be solved is composed of all the preceding equations with the adequate initial and boundary conditions:

At t = 0 : T (x, 0) = T 0 (x) in Ω At x = 0 : T (x, t) = T c for t ∈ (0, t end ] (Dirichlet) At x = l : ∂T ∂x = 0 for t ∈ (0, t end ] (Neumann) At t = 0 : P (x, 0) = P 0 (x) in Ω At x = 0 : P (x, t) = P atm for t ∈ (0, t end ] (Dirichlet) At x = l : ∂P ∂x = 0 for t ∈ (0, t end ] (Neumann) At t = 0 : n w (x, 0) = n w,0 (x) in Ω At x = 0 : n w (x, t) = n w,atm for t ∈ (0, t end ] (Dirichlet) At x = ξ(t) : n w (x, t) = P sat (T ) RT for t ∈ (0, t end ] (Dirichlet) At t = 0 : ξ(t) = ξ 0
where Ω is a domain in R, T represents the temperature, T 0 and T c are the initial and re temperatures respectively, P represents the pressure, P 0 and P atm are the initial and atmospheric pressure respectively, P sat is the saturation pressure of water vapor, ξ is the interface position, n w represents the molar density of water vapor and n w,0 (x) is a function presenting the initial molar density of water vapor.

For stability reasons at the beginning of numerical integration, n w,0 (x) is dened as a linear function between n w,atm at x = 0 and n w,ξ at x = ξ where n w,atm is the molar density of water vapor in the atmosphere.

Finally, the system of partial dierential equations can be written in vector form as follows: where A and B are T dependent variables, whereas q corresponds to the interface tracking equation 6.21.

                                                 ∂T ∂t = f (t, x, T, P ) + energy source term              A ∂T ∂t + B ∂P ∂t = g(t,
chapter 6

Notice that:

the second equation (conservation of mass) has been previously declared as valid in the gaseous zone only, we extend the equation to the whole domain by dening a degenerated form ∂P ∂x = 0 in the liquid zone. This extension is coherent with the fact that liquid water is incompressible and with the boundary condition for the pressure at the bottom of the domain.

the interfacial balance equations 6.15 and 6.20 are embedded in the corresponding eld equations. The additional terms will be detailed below in the numerical part.

the equation presenting the mass conservation of water vapor doesn't have an interfacial mass ux term (6.18) because there exists a boundary condition for n w at the interface.

the rst three equations are PDEs while the last one is an ODE.

Numerical method and discretized equations

As in previous chapters, in order to solve the system of PDEs 6.22, we chose to use the method of lines which is a way of approximating PDEs by ODEs where space and time discretizations are considered separately. The spatial discretization is performed using the vertex-centered nite volume method which conserves the mass locally and preserves continuity of uxes. The detailed discretized equations are below:

1. At each time step, the 1D energy conservation 6.3 is numerically integrated over the control volume and the discretized equation is given in equation 6.23 for 2 ≤ i ≤ N -1 whereas at i = i ξ an additional source term is added for the interface mesh cell only.

(ρC) e,i dT i dt ∆x - K(ρC) f,i µ f,i ∆x (T i+ 1 2 -T i-1 2 )(P i+1 -P i-1 ) = λ e | i+ 1 2 T i+1 -T i ∆x -λ e | i-1 2 T i -T i-1 ∆x + S E (6.23) S E =    φρ l L (ρC) e,i ∆x dξ dt i = i ξ 0 2 ≤ i ≤ N -1 (6.24)
S E is an energy source term. The value of S E at i = i ξ is due to the latent heat released by phase change (equation 6.15).

The authors in [START_REF] Eckert | Analysis of Heat and Mass Transfer[END_REF] provided the following law for the specic heat capacity:

C f = C l ,
in the liquid phase

C g = x w C w + (1 -x w )C a ,
in the gaseous phase (6.25)

where C l , C g , C w and C a are the specic heat capacity of liquid water, gas mixture, water vapor and air respectively. x w is the mass fraction of water vapor. Unfortunately, more complex laws are needed for the calculation of λ g and µ g .

These laws can be found in literature but for simplicity and as a rst approach, we neglect these complex laws and we use a simple linear relation between the components of the gas mixture.

2. The 1D ow equation (mass conservation equation 6.9) is numerically integrated over the control volume in the gas region only i.e. 2 ≤ i ≤ i ξ . The discretized equation for 2 ≤ i < i ξ is given by equation 6.26 whereas an additional term (equation 6.20) is added to the ow equation at i = i ξ as seen in equation 6.27.

1 T i β dP i dt - P i T 2 i β dT i dt - (K i+1 + K i ) (ρ f,i+1 -ρ f,i ) 2φ (µ f,i+1 + µ f,i ) ∆x 2 [P i+1 -P i ] (K i + K i-1 ) (ρ f,i -ρ f,i-1 ) 2φ (µ f,i + µ f,i-1 ) ∆x 2 [P i -P i-1 ] = S M (6.26) S M = ρ l ∆x dξ dt i = i ξ 0 2 ≤ i ≤ i ξ -1 (6.27)
Note that in the liquid region, the ow equation reduces to ∂P ∂x = 0.

3. The conservation of mass of water vapor is numerically integrated over the control volume in the gas region only. The discretized equation for 2 ≤ i < ξ is given by: 4. The spatial discretization of the interface-position equation is given by:

dn w,i dt - (K i+1 + K i ) 2φ (µ f,i+1 + µ f,i ) ∆x 2 [(n w,i+1 + n w,i ) (P i+1 -P i )] + (K i + K i-1 ) 2φ (µ f,i + µ f,i-1 ) ∆x 2 [(n w,i + n w,i-1 ) (P i -P i-1 )] - D w,a i φ∆x 2 (n w,i+1 -2n w,i + n w,i-1 ) = 0
dξ dt - n w,i ξ M w K µ g,i ρ l ∆x P i ξ -1 -P i ξ - D w,a i ξ M w ρ l ∆x n w,i ξ -1 -n w,i ξ = 0 (6.30) chapter 6
By using the methods of lines, functions T , P and n w are transformed into 3N vectors corresponding to their values at each discretization point in addition to ξ which is an unknown. We end up having 3N + 1 unknowns and 3N + 1 equations: 6.23, 6.26, 6.28 and one one more equation 6.30. After the spatial discretization, we obtain a system of ordinary dierential equations (ODE) which can be written in the general form:

dY dt = F (Y ) (6.31)
where Y = [T, P, n w , ξ] T and F (Y ) is a vector of dimension 3N + 1. The ODE system is solved by a specic ODE solver based on a BDF scheme as explained in 4.1.1.

Numerical example

In this section, we provide a numerical example where we applied the previous system of discretized equations. We consider 1D surface where we assume that the interface between the liquid water and air is initially at few millimeters depth (5 mm) for a domain of length 20 cm. The initial position of the interface is not set to 0 cm because this would lead to a singular problem (having two dierent boundary conditions at ξ = x = 0). The re temperature is T c = 250 • C (see gure 6.2), the interface is initially at 5 mm and the rst region contains air only whereas the second one contains liquid water only. Heating the system over 4 hours a certain amount of liquid water turns out into vapor and the position of the interface becomes deeper approximately 2.33 cm. The dotted line in gure 6.3 shows the nal position of the interface. The rst region now contains both air and water vapor. In the gures below, we present the results obtained upon heating a porous medium when coupling the heat diusion with water steam ow in presence of air. represents the molar fractions of of water vapor and air after 4 hours where we can see clearly that the amount of water vapor increases as depth increase while that of air decrease. Figure 6.5 shows the temperature prole at t f inal where we can notice that the phase change occurs at x = 2.33 cm which is actually the interface position. Although we are representing a pressure driven ow, we can see that the pressure increases slightly from 1 to 1.0014 atm (gure 6.8). During the rst heating stage (1 hour), the interface position moved from 5 mm to 11.6 mm so a distance equal to 6.6 mm where as in the rst cooling stage (4 hours) it moved around 2.55 mm. In the second heating stage (1 hour) the change in the position chapter 6 of interface was 2.0 mm while in the last cooling stage which took 4 hours it moved 2.15 mm. This can be seen clearly in gure 6.10. The drawbacks of this model is that it is incapable to simulate a cooling phenomenon of a porous medium. In order to validate our model using the previous model (with AHC), we assume that the physical properties (density, thermal conductivity, capacity, ...) of air are exactly equal to those of water vapor and we place initially the interface very close to the top boundary at 1 mm depth. The gures below presents the temperature history as well as the temperature prole at the nal time using the two models. Comparing the two temperature histories, we notice that they are close in shape but for sure there are some dierences since we are not using the same method to deal with phase change in addition to the fact that the initial interface position is not the same. The temperature proles (gures 6.13 and 6.14) are quite similar except for the fact that phase change happens earlier in the current enhanced model and at a smaller temperature. Also, at the bottom of the domain, after 4 hours, the temperature in the current enhanced model is slightly greater than that in the old one. Looking at both the pressure proles and the pressure histories (gures 6.15, 6.16, 6.17 and 6.18). We can realize that the pressure dierence in the current model is twice that in the old model (it increases from 1 Patm to 1.0015 Patm in the current model whereas it started with 1 Patm in the old model and ended up with 1.00075 Patm).

A Future Physical Model Presenting Phase Change in a Wet Porous Medium

In the frame work of ARPHYMAT project, many experiments were performed using dry sand. In fact, due to the humidity of the experiment room (where the humidity in Brittany often rises above 60% even in summer days), we suspect that the sand don't stay 100% dry and a small quantity of liquid water exist between the sand grains in the form of liquid meniscus, i.e. liquid water in our granular porous medium is in pendular regime (see table 1.5 and chapter 2). To model the phase change phenomenon, which constitutes of the evaporation of liquid water (pendular regime) present in the granular porous medium, an appropriate physical model must be used.

In literature, a great number of dierent physical models have been presented to model the heat and mass transfer in an unsaturated porous medium. Among them, In all previous studies of heat and mass transfer in an unsaturated porous medium, the funicular regime was considered and the eective thermal conductivity λ e was taken as a constant. We mentioned that we are interested in pendular regime and we have shown in chapter 2 of this thesis that a small quantity of liquid water can increase enormously the value of thermal conductivity. For this reason, we present in this second part of this chapter a physical model that describes heat and mass transfer in a wet porous medium (pendular regime) with λ e dependent on water content.

In our model, the quantity of liquid water located around the contact points of the wet zone is an unknown. Following intense heating, evaporation takes place at the beginning of the wet zone after the limit line between the dry and wet zones. The transport of the gaseous species occurs under pressure and concentration gradients. Liquid water is assumed immobile at pore scale because it is blocked by the capillary forces between the grains whereas it can move at a macroscale. Thermodynamic equilibrium is assumed to exist between vapor and liquid in the wet region. Air and water vapor obey the ideal gas law because they remain at low pressure. The unknowns in our model are:

1. The temperature T (x, t).

2. The pressure P (x, t).

3. The liquid fraction ψ l (x, t) where ψ l = liquid volume void volume .

4. The partial pressure of water vapor p w (x, t).

5. The velocity of gas u g (x, t).

6. The total curvature of the liquid meniscus κ(x, t).

The main governing equations are:

1. Mass conservation for water (in its two phases).

The main unknowns in this equation are ψ l and p w .

Mass conservation for air.

The main unknown in this equation is P = p w + p a where p a is the partial pressure of air.

3. Energy conservation containing a source term which is the Latent heat released due to phase change (liquid/vapor).

The main unknown in this equation is T .

4. Darcy's law: u g = -K µg ∇P . It provides a relation between u g and P . The auxiliary governing equations are:

1. Ideal gas law: P ρg = R M T where ρ g is the density of gas mixture and M is the molar mass of that mixture.

It provides a relation between T , P and ρ. This is a micro-macro law for the thermal conductivity from pore scale to macro scale. The calculation of the value of λ e is described in chapter 2.

It is important to mention that there is no need to write a boundary condition at the limit-line between the dry and wet zones because the zone detection is based on the value of ψ l which is an unknown that will be calculated at each time step.

In 6.2, we present a complex physical model for the phase change in a wet granular porous medium where λ e is a key parameter which undergoes a hysteresis behavior as function of liquid content (as shown in chapter 2). The hysteresis behavior will add more challenge to the model which is still under study and is a subject of future work. what was their function?

Answering the previous questions contributes to the understanding of ancient human behaviors and the reason behind using the re.

The work of this thesis aimed to simulate the forced evaporation of water in a saturated porous medium via a sequence of alternating heating/cooling stages. Many paths had been explored under the framework of ARHYMAT project as mentioned in the introduction.

In order to compare various numerical methods that can be used to model the phase change phenomenon (heat diusion neglecting the convection term), we used the latent heat accumulation method and the apparent heat capacity method. Even though the AHC method was studied in M. Muhieddine thesis, a special attention is paid to the value of the phase change temperature interval ∆T and the way to obtain the most optimum value that insures accuracy of the solution with few uctuations in the temperature history.

An adaptive mesh technique was applied to reduce the above-mentioned uctuations and to improve accuracy. This adaptive mesh is of node-based type where nodes are added and removed for an initial mesh. Renement is done by a number of recursive subdivisions near the interface (mixed cell in LHA method and phase change temperature in AHC method). A comparison between the performance of the two methods was conducted in a 1D coordinate system where an analytical solution exists for the phase change in a homogeneous medium (ice/water).

To perform simulations in higher dimensional congurations, the AHC method is used because it can be implemented easier than LHA.

Moreover, a home-made node-based adaptive 2D mesh algorithm Zohour is presented. This algorithm, developed by É. Canot, is composed initially of a regular squared mesh, at each level of subdivision, only nodes are added or removed according

Conclusion

to some prescribed rules and verifying some criteria, one of them may be based on the Hessian matrix and the other on the phase change proximity. In this way, maximum information is conserved and thus few interpolations are needed which reduces the error.

Zohour is tested against a reference solution obtained using a regular ne rectangular mesh. The results obtained are very interesting and the algorithm showed that it is quite ecient, though some optimizations are needed to be fully competitive with other mesh tools.

Another target of this work was the estimation of the thermophysical properties of the soil by inverse problem. Based on the knowledge of the temperature curves at selected positions of the altered soil, the thermal conductivity, the porosity and the volumetric heat capacity of the saturated soil are simultaneously identied. To solve this problem, the least square criterion was used to minimize the error function between the synthetic temperature measures and the calculated ones. The least square problem was solved by the Damped Gauss-Newton and Levenberg-Marquardt algorithms accompanied with a parameter's scaling technique. This inverse problem was presented in a 1D and a 3D axisymmetric coordinate systems respectively. The eect of the choice of the value of ∆T on the convergence of the inverse problem was emphasized and a new technique to ensure convergence was used: chaining the inverse problem. A sensitivity analysis was conducted showing that the porosity is the most sensitive parameter followed by the thermal conductivity and nally by the volumetric heat capacity.

The application of our inverse problem to real experimental data has been performed

for both saturated and dry soils. For the saturated case, the inverse problem diverged which was rendered to many reasons: the physical model is not close to reality.

the temperature distribution of the heating plate is not uniform.

the thermocouples positions are not accurate enough.

However, concerning the dry case, the inverse problem provided good results after the application of some statistical and box constraints on the sensors' positions. The results obtained (i.e. the value of the thermal diusivity of the soil) matched up to a cer- On the other hand, the calculation of the eective thermal conductivity for a granular medium (in pendular regime) formed of three components: solid, liquid water and air was performed. A simple model of two solid spherical grains is considered where liquid water is attached (without gravity) to these grains in the form of a liquid meniscus.

First, the exact shape of the liquid meniscus is obtained by integrating a dierential algebraic system which is based on the fact that the liquid meniscus is in mechanical equilibrium and thus it has a constant total curvature. Second, the total heat ux is calculated by solving the steady state heat conduction equation. The eective thermal conductivity is then deduced using Fourier law, assuming that the grains' distribution forms a periodic pattern. A strong hysteresis behavior for the eective thermal conductivity has been revealed when changing the humidity of the granular medium. This

shows that the eective thermal conductivity is a function of the temperature and the relative humidity of the porous medium. This result will be used as an input parameter for the future model that treats the phase change phenomenon in a wet unsaturated porous medium.

It is important to mention that the discretization of all the systems of equations in this work is done using the method of lines where time and space discretizations are considered separately. The spatial discretization is done using a centered nite volume scheme and the discretization in time is done via an ODE or DAE solver that uses a BDF implicit scheme. The Jacobian is calculated by a Computer Algebra System (like maple).

As for the possible perspectives of this work, we mention:

generalizing the calculation of the eective thermal conductivity to a more realistic grains' distribution by using some statistical information concerning:

the size of the gap between close grains;

the direction of the line joining the center of the close grains.

A rst work based on this idea can be found in [START_REF] Canot | Reorganization of liquid menisci between solid grains during evaporation/condensation of water in a wet porous medium[END_REF].

applying the future model presented in section 6.2: phase change in a wet unsaturated porous medium. We are interested in the unsaturated medium in pendular regime because:

studying general unsaturated medium is too vast;

it is well-adapted to the ARPHYMAT experiments.

Other applications of this new model concerns: salinization (long time evaporation at the soil surface leads to accumulation of salts which in return modify the porosity of the medium), ssuring (strong heating of the soil surface leads to fractures due to the mechanical response of the medium), etc.

Conclusion

Less important:

taking radiation, gravity and capillary forces into account in the future models dealing with phase change in a porous medium.

studying the eect of the variation of the parameters related to the use of Zohour:

The Hessian coecient and the thickness of the phase change zone.

AcknowledgementI

  would like to thank CNRS-Lebanon and INRIA-Rennes for the funding of my project throughout the three years. Specically, I'm thankful for Dr Charles Tabet, the responsible of the fellowship program at CNRS-Lebanon. I'm also grateful to his excellency Dr Adnan Mansour who supported me at the beginning of my journey. Without his help, I could never acheived my goals. I would like to thank Mr Denis Maillet, Professor at Lorraine University and Mrs Amel Ben Abda, Professor at ENIT (Tunis) to accept to be reporters for this manuscript. I would like also to thank Mr Sean Mcnamara, Assistant Professor at University of Rennes 1 to accept to be an examinator of this work. A special thank you to Mr Renaud Delannay, Professor at University of Rennes 1, for his time, advices and fruitful long discussions. I would like to thank my PhD advisors: Édouard Canot, Mohamad Muhieddine and Nabil Nassif for supporting me during these last three years.

  Conduction dans un milieu poreux sec + Évaporation (saturé ou partiellement saturé) + flux de vapeur (1D, 2D, 3D axisymmetrique)  𝑴 𝟐 : Separation des phases( avec diffusion de Gaz) zone 1: gaz seulement(vapeur + air) zone 2: liquide seulement (eau) • Modèle futur:  𝑴 𝟑 : Modèle insaturé en régime pendulaire (diffusion de gaz) zone 1: gaz seulement(vapeur + air) zone 2: liquid (eau) + gaz (vapeur +air)

  Maillage adaptive (LHA, AHC || 1D, 2D).  𝑴 𝟐 : Suivi d'interface (1D ). • Problème Inverse:  𝑴 𝟏 (sans écoulement de vapeur): DGN et LMA + données expérimentales, contraintes de domaine, contraintes statistiques(position des capteurs).

Figure 1 :

 1 Figure 1: Propriétés physiques données par Bonacina (lignes pointillées) et les propriétés physiques régularisés (lignes continues). L'aire de la zone hachurée est égale à la chaleur latente de changement de phase.

  de chaleur avec la présence d'un terme convectif, l'équation d'écoulement du gaz, la diusion de la vapeur d'eau dans le gaz et une équation qui représente la position de l'interface. Il faut noter que les archéologues s'intéressent à ce problème, car les anciens foyers sont constitués d'une séquence d'étapes de chauage et de refroidissement. Des résultats numériques sont présentés et commentés.Dans la deuxième partie, nous présentons, mais sans l'implémenter, un modèle physique qui peut être utilisé pour traiter le changement de phase dans un milieu poreux insaturé avec une très faible quantité d'eau liquide (régime pendulaire). Dans ce modèle, la conductivité thermique eective doit être connue, ce qui requiert l'utilisation des résultats du chapitre deux. Deux zones existent (mais non matérialisées par une interface) : l'un contient le gaz uniquement, l'autre contient de l'eau liquide se trouvant comme des ponts entre les grains (régime pendulaire). Dans les deux zones, la vapeur peut diuser dans l'air. La condition d'équilibre de la pression de vapeur conduit à un problème de complémentarité.

ν

  Unit outward vector normal to the boundary of a domain ξ Phase change interface location in 1D research in ow, heat and mass transfer in porous media have received increased attention during the past years. This is due to the fact that a variety of industrial, agricultural and energy production processes involve heat and mass transfer in porous medium having single or multiple uid phases. In many such systems, these uids are subject to evaporation, condensation, melting, freezing and transport due to pressure, temperature or concentration gradients. Examples include the drying of porous solids and soils (whether natural or not), geothermal energy production, thermal enhanced oil recovery etc. Because phase change in porous media covers dierent scales, modeling should be initiated at the rst scale of interest, typically the microscopic scale, as a basis for successively upscaling to larger scales. At microscopic scale, problems involving phase change are characterized by the presence of a boundary separating the phases involved and that changes position with time. The determination of this interface position together with the interfacial geometry are of high importance in phase change problems at microscopic level. As a conclusion, at microscopic scale, phase change leads to a free boundary value problem. At macroscopic scale, the calculation of the eective physical parameters involved in heat and mass transfer is a key point that results from the upscaling technique after their calculation at microscopic scale. The upscaling technique depends basically on a statistical point of view. The phase change is hidden in the microscopic scale behavior.

  Conduction in dry porous medium + evaporation (saturated or partially saturated) + vapor flow (1D, 2D, 3D axisymmetric)  𝑴 𝟐 : Separated phases ( with gas diffusion) zone 1: gas only (vapor + air) zone 2: liquid only (water) • Future model:  𝑴 𝟑 : Unsaturated model in pendular regime (gas diffusion) zone 1: gas only (vapor + air) zone 2: liquid (water) + gas (vapor +air) Experiments: • Replication (Field). • In lab (controlled experiments) : Adaptive mesh (LHA, AHC || 1D, 2D).  𝑴 𝟐 : Interface tracking (1D only). • Inverse problem:  𝑴 𝟏 (without vapor flow): DGN and LMA + experimental data, box constraints, statistical constraints (sensors position).

  International Conference on Porous Media & Annual Meeting (Interpore), May 18 -21, 2015, Padova, Italy. The last chapter concerns two other physical models where one of them is the main perspective of the upcoming future work. First, we present a physical model in an initially saturated porous medium in a 1D domain in presence of air to treat alternating heating/cooling stages (with phase change). The coupled equations presented in this model are: the heat diusion equation with convection term, the ow equation, the diffusion of vapor inside air and an equation that represents the position of the interface.

chapter 1 Figure 1

 11 Figure 1.1: Transition from pore-scale (microscopic) to continuum (macroscopic) representation (from [Rot12])

  λ e = φλ f + (1 -φ)λ sSerial and parallel models (see gure 1.3), assume that the porous medium is made up of solid and uid layers which are perpendicular (series model) and parallel (parallel model) to the heat ux resulting in the lower and upper limit values of λ e respectively.

Figure 1 . 3 :

 13 Figure 1.3: Figure (a) represents the situation where an arithmetic mean can be used where (b) represents the situation where the harmonic mean is suitable.

Figure 1 . 4 :

 14 Figure 1.4: Dierent water content regimes in a granular medium (from [ILHE00]).

Figure 1 . 5 :

 15 Figure 1.5: The dierent liquid regimes in a granular medium as function of liquid fraction in pores.

Figure 1 . 5

 15 Figure 1.5 and table 1.6 show the dierent water regimes in granular media with their physical properties.

chapter 2 Figure 2

 22 Figure 2.1: Real granular medium in the pendular regime: dry air is located around the grains (in light grey) whereas liquid water menisci (in blue) are present as liquid bridges between them.

Figure 2 . 2 :

 22 Figure 2.2: Simplication of the grains' assembly (with a periodic pattern). The dashed rectangle is the computational domain used in this work.

Figure 2 . 3 :

 23 Figure 2.3: Sketch of the cylindrical computational domain (taking into account all the symmetries). Boundary conditions are of homogeneous Dirichlet type on top and bottom sides, and of zero Neumann type on the external vertical sides.

Figure 2 .

 2 Figure 2.4 is a sketch used to demonstrate the calculation of the volume of liquid water V . The volume of the hashed slice of height dz is πr(z) 2 dz so the volume of both the blue and the brown parts is V b =

Figure 2 . 4 :

 24 Figure 2.4: A scheme of the liquid meniscus (solid grain in brown and liquid water in blue).A is an arbitrary point dened as the beginning of the integration process.

Figure 2 . 5

 25 Figure 2.5 is a sketch to demonstrate the calculation of the liquid and solid interface contact angle θ. Based on the knowledge of the position (r B and z B ) of point B

Figure 2 . 5 :

 25 Figure 2.5: A zoom of the liquid-solid interface (solid grain in brown and liquid water in blue).

Figure 2 . 6 :

 26 Figure 2.6: This represents the family of the liquid menisci when the starting point of integration r 0 varies, in the contacting state. Each curve is the boundary of the half of a bridge linking the two grains. Note that, for dierent r 0 (especially small values), two menisci might correspond to the same liquid volume.

Figure 2

 2 Figure 2.7: This represents the family of the liquid menisci when the volume varies, in the non-contacting state. Each curve is the boundary of a sessile drop, and the same symmetric drop is on the top of the other grain (not represented).

chapter 2 Figure 2 . 9 :

 229 Figure 2.9: Example of contour curves of the temperature obtained after a numerical computation on a 500 × 500 mesh, in the case /R 0 = 0.1, θ = 30 deg. Due to Dirichlet boundary condition (bottom: T = 0 and top: T = 1), each contour line from bottom to top corresponds to the temperature: T = 0.1, 0.2, ..., 0.9 respectively.

Figure 2

 2 Figure 2.10: Thermal heat ux w.r.t. the liquid fraction, for the contacting state (blue) and the non-contacting state (red).

chapter 2 Figure 2

 22 Figure 2.11: This zoom of gure 2.10 shows with more evidence the hysteresis behavior when liquid fraction is increasing or decreasing. The jump in the heat ux occurs at the magenta dashed line; this jump is more pronounced in the increasing case, i. e. when liquid water increases. The arrows show the direction in the hysteresis cycle.

Figure 2 . 12 :

 212 Figure 2.12: Surface area w.r.t liquid volume (for small menisci), for the contacting state, shows that the same liquid volume corresponds to two dierent surface area.

Figure 2 . 13 :

 213 Figure 2.13: Thermal heat ux w.r.t. the liquid fraction, for the contacting state only, when θ = 30 deg. The curves show the inuence of the distance .

Figure 2 . 14 :

 214 Figure 2.14: This is a zoom of gure 2.13. It shows more clearly that only the cases corresponding to = 0 lead to a hysteretic behavior.

Figure 2 Figure 2

 22 Figure 2.15: Thermal heat ux w.r.t. the liquid fraction, for the contacting state only, when /R 0 = 0.1. The curves show the inuence of the contact angle θ, which is less important than those of the gap .

Chapter 3 Numerical

 3 Methods for Transfer of heat with phase change The special diculty posed by phase-change problems is the latent heat release on the moving interface whose position is usually an unknown function of time and needs to be determined as a part of the solution. In this chapter, we chose to handle the Stefan problem which is probably the simplest mathematical model of phenomenon of phase change. When a phase change takes place, a latent heat is either absorbed or released while the temperature of the material changing remains constant. The strong formulation of this model corresponds to a free boundary value problem involving a parabolic diusion equation for each phase and a transmission condition prescribed at the interface separating the phases.

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: Analytical solution of heat transfer with phase change (fusion). The gure represents the temperature prole at a xed time t.

  mK and λ s = 2.220 W/mK C l = 4.2 × 10 3 J/kgK and C s = 2 × 10 3 J/kgK ρ l = ρ s = 1000 kg/m 3 Part I: Phase change problem in 1D coordinate system 3.1 The Latent Heat Accumulation Method (LHA)

Figure 3 . 3 :

 33 Figure 3.3: The volumetric cells are equally spaced.

Figure 3 . 4 :

 34 Figure 3.4: Spatial temperature prole for number of mesh cells N = 640 at t f inal = 50 h. ∆t = 1.74 × 10 s. Analytical solution in red and numerical solution in blue using LHA method. The average L 2 relative error over the whole domain is 1.82 × 10 -2 %. In

Figure 3 . 5 :

 35 Figure 3.5: Temperature history at x = 5 cm. Number of mesh cells N = 640 and ∆t = 1.74 × 10 s. Analytical solution in red and numerical solution in blue using LHA method. The average L 2 relative error over the whole duration (t max = 50 h) is 4.11 × 10 -1 %. CPU time = 7.64 × 10 -2 s

chapter 3 Figure 3

 33 Figure 3.6: Physical properties given by Bonacina in doted lines and regularized physical properties in continuous lines. The area of the hashed zone is equal to the Latent heat of phase change.

Figure 3 . 7 :

 37 Figure 3.7: Quality factor versus the magnitude of the phase change temperature interval ∆T for dierent number of mesh cells N . Taken from [MCM09].

Figure 3 . 8 :

 38 Figure 3.8: ∆T optimum versus h. The blue dotted line (polynomial t of the data) shows that ∆T optimum ≈ kh.

4 Figure 3 4 Figure 3 .

 4343 Figure 3.9: Variation of temperature history at x = 5 cm. Analytical solution in magenta.

Figure 3 . 12 :

 312 Figure 3.11: Variation of the relative error in temperature prole and temperature history as function of ∆T .

Figure 3 . 13 :

 313 Figure 3.13: Spatial temperature prole at t f inal = 50 h. Analytical solution in red and numerical solution in blue using AHC method with an Euler explicit scheme. The average L 2 relative error (e x ) over the whole domain is 1.86 × 10 -2 %. In this gure we represent 1 8 of the whole domain.

Figure 3 . 14 :Figure 3 .

 3143 Figure 3.14: Temperature history at x = 5 cm. Analytical solution in red and numerical solution in blue using AHC method with an Euler explicit scheme. The average L 2 relative error (e t ) over the whole duration (t max = 50 h) is 3.68 × 10 -1 %. CPU time = 2.18 × 10 -1 s

Figure 3 .

 3 Figure 3.16: Temperature history at x = 5 cm. Analytical solution in red and numerical solution in blue using AHC method with a BDF implicit solver. The average L 2 relative error (e t ) over the whole duration (t max = 50 h) is 3.04 × 10 -1 %. CPU time = 3.39 × 10 -1 s

Figure 3 .

 3 Figure 3.17: A diagram representing a variety of mesh types with their strengths and drawbacks.

Figure 3 .

 3 Figure 3.18: An example of a non-conforming cell-based adaptive mesh. From http: //www.dealii.org/developer/doxygen/deal.II/step_12.html

Figure 3 .

 3 Figure 3.20: Our adaptive mesh. The dashed lines represent the boundaries of the cells while the continuous lines represent the nodes of the cells. Continuous black lines represent the nodes of the basic mesh cells. The blue, red, green and magenta lines represent respectively 1, 2, 3 and 4 level of subdivisions.

Figure 3 .Figure 3 .

 33 Figure 3.21: LHA: spatial temperature prole at t f inal = 50 h. Analytical solution in red and numerical solution in blue. N basic = 80. Number of subdivisions = 3. ∆x basic = 0.05 m. e x = 1.66 × 10 -2 %.

3 Figure 3 .

 33 Figure 3.23: The error in spatial temperature prole (e x ) decreases from 6.37 × 10 -2 % to 1.18×10 -3 % as the basic number of mesh cells (number of subdivisions = 2) increase from 80 to 1280 whereas it decreases from 1.66 × 10 -2 % to 7.16 × 10 -4 % for number of

3 Figure 3 .

 33 Figure3.24: The error in temperature history (e t ) at x = 5 cm decreases from 9.33 × 10 -1 % to 3.21 × 10 -2 % as the basic number of mesh cells (number of subdivisions = 2) increase from 80 to 1280 whereas it decreases from 4.05 × 10 -1 % to 1.41 × 10 -2 %

Figure 3 .

 3 Figure 3.25: In LHA: Variation of error in temperature prole (e x ) at t f inal = 50 h as function of CPU time in seconds for dierent basic meshes and subdivisions.

Figure 3 .

 3 Figure 3.26: In LHA: Variation of error in temperature history (e t ) at x = 5 cm as function of CPU time in seconds for dierent basic meshes and subdivisions.

Figure 3 .Figure 3 .Figure 3 .Figure 3 .

 3333 Figure 3.27: Spatial temperature prole at t f inal = 50 h. Analytical solution in red and numerical solution in blue (using Euler explicit method). N basic = 80. Number of subdivisions = 3. ∆x basic = 0.05 m. e x = 1.92 × 10 -2 %.

3 Figure 3 .

 33 Figure 3.31: Using Euler explicit scheme. e x decreases from 4.4×10 -2 % to 7.89×10 -3 % (level of subdivions = 2) as N basic increases from 80 to 1280.

3 Figure 3 . 3 Figure 3 .

 3333 Figure 3.32: Using Euler explicit scheme. e t at x = 5 cm decreases from 9.00 × 10 -1 % to 7.43 × 10 -2 % (level of subdivions = 2) as N basic increases from 80 to 1280.

3 Figure 3 .

 33 Figure 3.34: Using ODE solver. e t at x = 5 cm decreases from 7.17 × 10 -1 % to 4.26 × 10 -2 % (level of subdivions = 2) as N basic increases from 80 to 1280.

Figure 3 .

 3 Figure 3.35: Using AHC with Euler explicit method: Variation of error in temperature prole (e x ) at t f inal = 50 h as function of CPU time in seconds for dierent basic meshes and subdivisions. For N basic = 80 with number of subdivisions = 3, L 2 error is e x = 1.92 × 10 -2 % and CPU = 3.24 × 10 -2 s.

Figure 3 .Figure 3 .

 33 Figure 3.36: Using AHC with Euler explicit method: Variation of error in temperature history (e t ) at x = 5 cm as function of CPU time in seconds for dierent basic meshes and subdivisions. For N basic = 80 with number of subdivisions = 4), L 2 error is e t = 2.02 × 10 -1 % and CPU = 1.33 × 10 -1 s.

Figure 3 .

 3 Figure 3.38: Using AHC with an automatic ODE solver: Variation of error in temperature history (e t ) at x = 5 cm as function of CPU time in seconds for dierent basic meshes and subdivisions. For N basic = 80 with number of subdivisions = 4, L 2 error is e t = 1.85 × 10 -1 % and CPU = 1.49 × 10 -1 s.

Figure 3 .Figure 3 .

 33 Figure 3.39: Using both LHA and AHC methods: variation of error in temperature history (e t ) as function of CPU time for dierent N basic = 80, 160 and dierent subdivisions.

Figure 3 .

 3 Figure 3.41: Variation of error in temperature history as function of CPU time for the three methods (three levels of subdivision).

Figure 3 .Figure 3 .Figure 3 .

 333 Figure 3.42 represents our computational domain where we use Dirichlet and Neumann boundary conditions

  nodes are added or removed creating new cells. Each cell is dened by a shape, a position and subdivision level. 10 dierent shapes are counted for interior cells; they are listed in gure 3.45. The other cells (that rely on boundary nodes or corner nodes) may have even dierent forms: we speak of half-cell for boundary nodes and a quarter-cell for corner nodes. C. Beaucé worked on the optimization of the calculation of the Hessian during his Masters 1 training at INRIA [Bea15].

  node-basedThe numbers in black located inside the cells indicate the side length in terms of the distance from the central node.

Figure 3 .

 3 Figure 3.45: Dierent shapes of mesh cells obtained after subdivisions (due to the addition of nodes). All these gures can be rotated by 90 • , 180 • or 270 • . The same series of shapes inclined at 45 • also exist.

Figure 3 .

 3 Figure 3.46: Variation of L 2 error in temperature prole at nal time as function of CPU time for re-meshing performed at dierent time steps, using only two levels of subdivisions.

Figure 3 .

 3 Figure 3.47: Variation of L 2 error in temperature history at a xed position as function of CPU time for re-meshing performed at dierent time steps, using only two levels of subdivisions.

Figure 3 .Figure 3 .

 33 Figure 3.48: The mesh grid created by Zohour using N basic = 70 × 70 and a level of subdivision = 1. The total number of mesh cells is 6075.

Figure 3 .Figure 3 .

 33 Figure 3.50: Variation of L 2 error in temperature prole at nal time as function of number of basic mesh cells.

Figures 3 .

 3 Figures 3.50 and 3.51 shows the variation of e x and e t as function of basic number of mesh cells for dierent levels of subdivisions.

Figure 3 .Figure 3 .

 33 Figure 3.52: Temperature history at three dierent sensors. N = 70 × 70 and level of subdivisions = 1.

3 )

 3 Note that (ρC) f = ρ f C f where ρ f and C f are dened as in equations 4.8 and 4.6 respectively. (ρC) s = ρ s C s where the thermophysical properties of the solid matrix (ρ s = 1500 kg/m 3 , C s = 1300 J/kgK and λ s = 0.756J/kg). The eective conductivity in equation (4.3) is calculated using Kunii and Smith model [MHB+ 12], β = 0.895 chapter 4

. 5 )Figure 4 . 1 :

 541 Figure 4.1: Smoothed physical properties used in our model (AHC method) are in continuous lines and non-smoothed properties are in dashed lines. Smoothing is required to alleviate the singularity while dierentiating the physical parameters.

Figure 4 .

 4 Figure 4.2: 1 D control volume

( 4 .

 4 13): ∇S(p k ) = J T (p k )r(p k ) = 0(4.13) where J(p k ) i,j = ∂r i (p k ) ∂p k j , i = 1, 2, ..., N and j = 1, 2, 3. The sensitivity matrix, chapter 4

i

  for i = 1, ..., M and f = 1, ..., F 5-Calculate r (k) and the Sensitivity matrix J knowing that J = J.D chapter 4

Figure 4 .

 4 Figure 4.4 represents the variation of the nal residue as function of the number of mesh cells. We can easily notice that the residue decrease as number of mesh cells increase which assures the consistency of our method.

Figure 4 . 3 :

 43 Figure 4.3: Variation of the conductivity (in W/m.K) as function of iteration number for 120 mesh cells. The red line represents the exact value of λ s .

Figure 4 . 4 :

 44 Figure 4.4: Variation of residue (absolute value in Celsius) as function of number of mesh cells. The method is consistent (the error decreases as number of mesh cells increase). The order of convergence is O(h 1.116 ).

4. 3 . 1

 31 Introduction to LMALevenberg[START_REF] Levenberg | A method for the solution of certain nonlinear problems in least squares[END_REF] and Marquardt[START_REF] Marquardt | An algorithm for least squares estimation of nonlinear parameters[END_REF] proposed a very elegant algorithm for the numerical solution of equation (4.10). However, most implementations are either not robust, or do not have a solid theoretical justication. Moré[START_REF] Moré | The Levenberg-Marquardt algorithm: implementation and theory[END_REF] presented a robust and ecient implementation of a version of the Levenberg-Marquardt and show that it has strong convergence properties. In addition to robustness, the main features of this implementation are the proper use of implicitly scaled variables and the choice of the Levenberg-Marquardt parameter via a scheme due to Hebden[START_REF] Hebden | An algorithm for minimization using exact second derivatives[END_REF]. The implementation of LMA by Moré that is contained in Minpack[START_REF] Moré | User guide for MINPACK-1[END_REF] has proven to be very successful in practice. Several factors make LMA preferable to DGN: rst is that LMA possesses an embedded scaling technique, second it is well dened even when chapter 4

Figure 4 . 5 :

 45 Figure 4.5: Variation of the volumetric heat capacity (in J/m 3 K) as function of iteration number using LMA (Same number of mesh cells in both the forward and inverse problems). The red line represents the exact value of (ρC) s .

Figure 4 . 6 :

 46 Figure 4.6: Variation of the porosity as function of iteration number using LMA (Same number of mesh cells in both the forward and inverse problems). The red line represents the exact value of φ.

Figure 4 . 7 :

 47 Figure 4.7: Variation of the conductivity (inW/mK) as function of iteration number using LMA (Same number of mesh cells in both the forward and inverse problems). The red line represents the exact value of λ s .

Figure 4 .

 4 Figure 4.8: Cylindrical Coordinate System (In our case θ is assumed constant

chapter 4 numerical

 4 tests using Robin boundary condition have been performed showing slight dierences in temperature especially near the re.

Figure 4 . 9 :

 49 Figure 4.9: The computational domain in 3D axisymmetric coordinate system. The red part represent the re (re length = 15 cm). N represents Neumann boundary condition while D represents Dirichlet boundary condition

Figure 4 .

 4 Figure 4.10: 3D axisymmetric control volumes

( 4 .

 4 41)which leads to the general sensitivity equation below:

( 4 .

 4 45)These three sensitivity equations (4.43), (4.44) and (4.45) are completed with adequate initial and boundary conditions. W , R and Z are the unknowns of the sensitivity equations and T is the temperature. The obtained system of coupled equations (heat diusion equation + 3 sensitivity equations) is a nonlinear system of partial dierential equations. To solve this system, we use the same numerical strategy used to solve the 1D inverse problem.

(

  ρC) s (J/m 3 K) λ s (W/mK)

Figure 4 .

 4 Figure 4.11: Variation of residue as function of number of mesh cells.

Figure 4 . 12 :

 412 Figure 4.12: Variation of residue as function of iteration number.

Figure 4 . 13 := 5 .

 4135 Figure 4.13: Sensitivity Analysis: ratio of singular values σ k σ 1 . σ 1 = σ φ so σ (ρC)s σ 1 = 5.513 × 10 -3 and σ λs σ 1 = 8.733 × 10 -2 .

Figure 5 . 5 and

 55 1 represents the experimental setup performed by Szubert chapter Bergonzi in 2010.

Figure 5 . 1 :

 51 Figure 5.1: Experimental hearth. The sand is fully saturated with water.

(

  ρC) s (J/kg.K ) λ s (W/m.K)

(

  ρC) s (J/kg.K ) λ s (W/m.K) Physical properties of the soil obtained by inverse problem using real experimental data.

Figure 5 . 2 :

 52 Figure 5.2: The temperature of the heating plate is not uniform. In the experiment shown in this gure, it ranges between 250 and 320 degree Celsius.

Figure 5 .

 5 Figure 5.3 represents the experimental setup.

Figure 5 . 3 :

 53 Figure 5.3: Experimental hearth.

Figure 5 . 4 :

 54 Figure 5.4: Variation of slope as function of time (extracted from A. Cordero postdoctoral report [Cor13]) .

8 3 αtFigure 5 . 5 :

 355 Figure 5.5: The dierence in temperature prole for 1D and 3D-axisymmetric cases.

Figure 5 . 6 :Figure 5

 565 Figure 5.6: Laplace transform curves for two dierent sensors

Figure 5 . 9 :

 59 Figure 5.9: Laplace transform curves for dierent sensors at dierent depths. The information didn't reach yet the sensor: S 8 for p > 3.5 × 10 -3 (short time and far from re).

Figure 5 .Figure 5 .Figure 5 . 12 :

 55512 Figure 5.10: Diusivity Determination for sensors at dierent depths. The non-physical behavior of the curve representing sensors 7 and 8 comes from the deciency in the Laplace transform curve at sensor 8. (see gure 5.9)

Figure 5 .

 5 Figure 5.13: 1D heating process, showing the horizontal isothermal lines; there is a null sensitivity along the r-axis.

chapter 5 Figure 5 . 14 :

 5514 Figure 5.14: In a 3D-axisymmetric heating process, the isothermal lines are very close to ellipsoidal curves; there is a null sensitivity in a direction tangent to these curves.

n i=1 δ r i = 0 and n i=1 δ z i = 0 ( 5

 05 .20) for all sensors i.Lastly, we suspect that a bias (actually a shift, or an oset; not a statistical bias) may be present in the measures (see the second remark p. 148). In any statistical investigation, we can always attribute some of the variation in data to measurement error, part of which can result from the measurement instrument itself. But human mistakes, especially recording errors (e.g., misreading a dial, incorrectly writing a number, not observing an important event, misjudging a particular behavior), can also often contribute to the variability of the measurement and thus to the results of a study. In our experiments, the way the sensors' positions were measured might be aected by a bias (see gure 5.15).

Figure 5 .

 5 Figure 5.15: A bias in the measures might result from the way the sensors' position were calculated.

Figure 5 .

 5 Figure 5.16: Experimental positions of 10 sensors in blue. The red bar represents the re (unit is meter for both axes).

  Figure 5.17: Initial temperature prole of the rst experiment (in Celsius), showing that the temperature is nearly uniform.

Figure 5 .

 5 Figure 5.18: Variation of diusivity (in m 2 /s) as function of mesh size.The mesh size varies from 30 × 50 to 300 × 500. The value of diusivity converges to 1.73 × 10 -7 m 2 /s)

Figure 5 .

 5 Figure 5.20 represents the experimental (old) positions of sensors, the new numerical positions and the direction of displacement.

Figure 5 .

 5 Figure 5.20: The experimental positions of sensors are in blue and the new calculated positions are in red (unit is meter for both axes).

  chapter 5

Figure 5 .

 5 Figure 5.21: Experimental positions of 19 sensors in blue. The red bar represents the re (unit is meter for both axes).

Figures 5 .

 5 Figures 5.23 and 5.24 presents that variation of thermal diusivity as function of mesh size and the variation of residue as function of mesh size respectively for the case where we used a uniform initial temperature. Unfortunately, the curves are not well-convergent which incite us to use an exponential t instead.

Figure 5 .Figure 5 .

 55 Figure 5.23: Variation of diusivity (in m 2 /s) as function of mesh size using constant initial temperature. Mesh varies from 30 × 50 to 300 × 500

Figure 5 .

 5 Figure 5.26: Variation of residue (absolute value in Celsius) as function of mesh size using an exponential t for initial temperature.

Figure 5 .

 5 Figure 5.27: The experimental positions of sensors are in blue and the new calculated positions are in red

  Treat Phase Change in a Porous MediumThe previous work presented in this thesis and in M. Muhieddine thesis[START_REF] Muhieddine | Simulation numérique des structures de combustion préhistoriques[END_REF] presents phase change problem using AHC method with one uid constituent (water) whether liquid water or water vapor where phase selection depends on temperature (i.e. if

Figure 6 . 1 :

 61 Figure 6.1: Physical description of the domain.

  (6.28) At i = i ξ the equation reduces to the boundary condition:n w,i = P sat (T i ) RT i (6.29)which represents the thermodynamics equilibrium at the interface.

Figure 6 . 2 :

 62 Figure 6.2: The physical domain at initial time t 0 = 0 shows the initial constituents in the medium (air + liquid water) with the initial position of the interface in dotted line.

Figure 6 . 3 :

 63 Figure 6.3: The physical domain at the nal time t f = 4 h shows the existence of water vapor in addition to initially existing constituents with the nal position of the interface.

Figure 6. 4 ,

 4 represents the temperature history of three sensors, the rst sensor (in blue) attains phase change approximately at t = 26 mins and T = 95.5 • C.

  Figure 6.6, shows the evolution of the phase change interface position as function of time while gure 6.7

Figure 6 . 4 :

 64 Figure 6.4: The temperature history at three dierent depths 0.75 cm, 1.25 cm and 2 cm.

Figure 6 . 5 :

 65 Figure 6.5: The temperature prole at the nal time.A slope discontinuity is observed at the phase change as expected.

Figure 6 . 7 :

 67 Figure 6.7: The variation of the molar fraction of both air and water vapor as function of x at the nal time. The sum of the two quantities is always one.

Figure 6 . 8 :Figure 6 . 9 :

 6869 Figure6.8: The pressure prole at t f inal

  Figure 6.10: The change in the interface position during a sequence of respectively two alternating heating and cooling stages.

Figure 6 Figure 6 . 12 :

 6612 Figure 6.11: Temperature history for 3 different sensors using the old model.

  Looking at the rst sensor (blue), the phase change in the old model takes place at t = 38 min and T = 102 • C whereas it happens earlier in the current enhanced model at t = 36 min and at smaller temperature T = 98.5 • C.

Figure 6. 13 :

 13 Figure 6.13: Temperature prole at t f inal = 4 h using the old model.

Figure 6. 14 :

 14 Figure 6.14: Temperature prole at t f inal = 4 h using our current enhanced model.

  bottom: x = 2.00E-01 m.

Figure 6 .

 6 Figure 6.15: Pressure history at the bottom of the domain x = 20 cm using the old model.

Figure 6 .

 6 Figure 6.16: Pressure history at the bottom of the domain x = 20 cm using our current enhanced model.

Figure 6 .

 6 Figure6.17: Pressure prole at t f inal = 4 h using the old model.

Figure 6 .

 6 Figure 6.18: Pressure prole at t f inal = 4 h using our current enhanced model.

Udell [ Ude85 ]

 Ude85 presented a study of heat and mass transfer (one-dimensional steady state problem) in a porous medium saturated with liquid and vapor phases of a single component uid. Gravity, capillarity, multiphase ow and phase change eects are included. In a recent study of Bridge et al. [BBWW03], the authors present an extension of what was presented by Udell by adding an energy equation to the system. In these two studies the authors were interested in the steady state case. Another interesting study is that presented by Dayan and Gluekler [DG82] where a model for heat and mass transfer has been developed for a porous medium containing water (funicular regime), vapor and air. The physical model presented by Dayan is divided into two zones: dry and wet. A similar study is that of Min and Emmons [ME72] who studied the drying of a porous medium in funicular regime where two mechanisms of mass transfer were considered: the pressure driven convective ow governed by Darcy's law and the molecular diusion of the gaseous components.

  Figure 6.19 presents the physical domain in our new model where the dry zone is made up of both air and vapor chapter 6 whereas the wet zone contains liquid water and air. The value of the eective thermal conductivity is calculated depending on quantity of liquid water (liquid fraction) and on the temperature.

Figure 6 .

 6 Figure 6.19: The physical domain representing the unsaturated wet model. The limit line is not computed directly by the model, it will be deduced from the equilibrium condition at the interface.

Figure 6 .

 6 Figure 6.20: A diagram representing the expected results of the phase change problem showing the existence of an evaporation and condensation zones in the wet zone.

5 .

 5 Liquid/vapor equilibrium at interfaces, written at all positions:if ψ l > 0 then p w = k p sat (T ) otherwise, if ψ l = 0 then p w < k p sat (T )k is given by Kelvin's law and includes the total local curvature of each liquid meniscus κ.

  6. A relation between the liquid fraction ψ l and the total curvature of the liquid meniscus κ. This relation is a result of the use of some statistical information in the case of a random irregular grains' distribution (generalization of the results of chapter 2).

  2. A relation joining λ e with ψ l : λ e = f (ψ l ) taking into consideration an irregular grain distribution.

  this thesis was mainly the mathematical reproduction of the prehistoric res in order to answer the following questions: what was the duration of the re? what was their mode of operation?

  tain limit those obtained by other methods like Laloy & Massard and Laplace transform. Another part of this work simulate an interesting archaeological problem: Alternating heating/cooling stages of a porous medium. To simulate this problem, a separated phase model in a 1D saturated porous medium (in presence of air) is presented. The AHC method is not able to treat the cooling stage in a water saturated porous medium so a new model is required. The coupled equations used in this model are: the heat diusion equation with the convection term, the ow equation, the diusion of vapor inside air and an equation that represents the position of the interface. We used a front tracking method where the latent heat released due to phase change is embedded in the energy equation only at the phase change location.

  en 2014-2015. Le maillage initial est composé d'éléments carrés réguliers, à chaque niveau de subdivision, seulement des n÷uds sont ajoutés ou enlevés selon certaines règles prescrites et en vériant certains critères, dont l'un est basé sur la matrice du Hessien et l'autre sur la proximité avec

l'endroit où se passe le changement de phase. De cette façon, nous conservons un maximum d'informations, à savoir, le maillage initial reste xe et très peu d'interpolations sont nécessaires ce qui réduit l'erreur. Les mailles sont ensuite créées par un diagramme de Voronoi an de former de nouvelles formes de mailles à chaque nouveau niveau de subdivisions. Ce maillage adaptatif est appliqué à la méthode AHC pour résoudre le problème de fusion dans un système de coordonnées 2D.

Table 1 :

 1 Propriétés physiques du sol obtenu par le problème inverse.

		(ρC) s (J/kg.K ) λ s (W/m.K)	φ
	exacte	1.95 × 10 6	0.756	0.20
	Valeur initiale	2 × 10 6	0.8	0.18
	Calculé dans [Muh09]	2 × 10 6	0.7696	0.1979
	Calculé dans ce travail	1.9497 × 10 6	0.7559	0.2000

  Choice of ∆T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Estimation of α and sensors' positions . . . . . . . . . . . . . . .

	5.2.2	Methods used to nd an approximate value of diusivity α . . .
		5.2.2.1	Laloy and Massard (Done by ARPHYMAT team) . . .
		5.2.2.2	Laplace Transform (Done by R. Delannay and É. Canot)
	5.2.3	Inverse Problem in 3D-axisymmetric coordinate system	. . . . .
		5.2.3.1	Forward Problem . . . . . . . . . . . . . . . . . . . . . .
		5.2.3.2	Method of Resolution . . . . . . . . . . . . . . . . . . .
		5.2.3.3	Governing Equations and Numerical Strategy . . . . . .
		5.2.3.4	Results . . . . . . . . . . . . . . . . . . . . . . . . . . .
	5.2.4
		5.2.4.1	Numerical Strategy . . . . . . . . . . . . . . . . . . . .
		5.2.4.2	Solving the constrained inverse problem . . . . . . . . .
	3.2.1	Numerical Strategy . . . . . . . . . . . . . . . . . . . . . . . . . .
	3.2.2
	3.2.3	Results using AHC method with uniform mesh . . . . . . . . . .

3.3 Adaptive Mesh Renement . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 LHA with our adaptive mesh . . . . . . . . . . . . . . . . . . . . 3.3.1.1 Results and consequences . . . . . . . . . . . . . . . . . 3.3.2 Adaptive mesh applied to AHC method . . . . . . . . . . . . . . 5.2.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Other Approaches to Treat Phase Change in a Porous Medium 6.1 Phase change in a saturated porous medium in presence of air . . . . . . 6.1.1 Physical description and mathematical modeling . . . . . . . . . 6.1.2 Numerical method and discretized equations . . . . . . . . . . . . 6.1.3 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 A Future Physical Model Presenting Phase Change in a Wet Porous Medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 1: Examples of average porosities of some substances.

	Substance	Porosity
	Foam metal	0.98
	Silica grains	0.65
	Soil	0.43 -0.54
	Sand	0.37 -0.50
	Limestone	0.04 -0.10

Materials with low porosity are usually less permeable and typically have smaller pores, making it more dicult for gas or liquid to pass through them while materials with high porosity have large pores and are easily permeated.

Table 1

 1

	.2: Permeability of some common materials.
	Substance	Permeability (m 2 )
	Limestone	2 × 10 -15 -4.5 × 10 -14
	Sand	2 × 10 -11 -1.8 × 10 -10
	Silica Powder

Table 1

 1

	.4: Density of important substances.
	Substance	Density (kg/m 3 )
	Fresh water	1000
	Sea water	1025
	quartz	2650
	Sand, Silica	1538
	Air	1.2

Ice 934

Table 1 .

 1 6: Granular media with various amount of liquid[START_REF] Mitarai | Wet granular materials[END_REF].

	Content of pores	State	Physical description

Table 2 .

 2 1: Role of sphere distance and contact angle θ on the hysteresis

	/R 0 θ (degree) Hysteresis Length Flux jump
	0.05	10	0.0192	0.165
	0.05	30	0.002027	0.041026
	0.05	45	0.0008265	0.021695
	0.1	10	0.10684	0.285
	0.1	30	0.013434	0.0881
	0.1	45	0.005633	0.048

Figure 2.13 and 2.14 (resp. 2.15 and 2.16) shows the variations of the heat ux curves with respect to the distance (resp. w.r.t θ). They show clearly that the inuence of the distance is more important than that of the contact angle θ.

  and k is a constant chosen in a way to obtain good accuracy with fewer uctuations in the temperature prole. The value of ∆T that insures accuracy of the solution with few uctuations is called ∆T optimum .In order to nd ∆T optimum , the authors in[START_REF] Muhieddine | Various approaches for solving problems in heat conduction with phase change[END_REF] introduced a quality factor which combines accuracy and smoothness (see equation 3.20) and plotted it as function of ∆T for dierent number of mesh cells (see gure 3.7).

	Quality Factor = Error(T ) + 0.025 Error(	∂T ∂t	)

It seems reasonable to state that ∆T is proportional to h [MCM09] (∆T = kh) where h is the mesh chapter 3 size

  optimum 4 are the closest to the analytical solution but they suer from uctuations which are slightly less remarkable using ∆T optimum . Moreover, we can see that the red curve (∆T optimum × 4) is smooth with no uctuations but the numerical solution is far from the analytical one. The error in temperature history as well as in temperature prole as function of ∆T chosen are shown in gure 3.11 where we notice that the error is minimum in temperature history for ∆T = ∆T optimum whereas the error in temperature prole is minimum for ∆T =

	∆T optimum 2	. Having a look on the CPU time required to

solve the fusion problem using AHC method (see gure 3.12) with dierent ∆T values, we remark that the CPU time increases as the value of ∆T decreases. Depending on all what preceded, to have good accuracy in both temperature history and temperature prole with a moderate cost we will use ∆T = ∆T optimum in all what follows.

Table 4

 4

	.1: Physical properties of the soil obtained by inverse problem with or without
	scaling.
		(ρC) s (J/m 3 K) λ s (W/mK)	φ
	exact	1.95 × 10 6	0.756	0.20
	initial guess	2 × 10 6	0.8	0.18
	calculated without scaling calculated with scaling	2 × 10 6 1.9497 × 10 6	0.7696 0.7559	0.1979 0.2000

Table 4 .

 4 s , φ and λ s recede from the exact solution and the value of residue increases from 9.588 to 26.7198 as the value ∆T increases from 3 × ∆T optimum to 10 × ∆T optimum (see table4.2). The results in table 4.2 are obtained by using (ρC) s = 3 × 10 6 , λ s = 0.4 and φ = 0.12 as initial guesses. If the results obtained in the rst row (∆T = 3 × ∆T optimum ) are used as initial guesses then the inverse problem will converge to the exact values using ∆T optimum if we use the same mesh size in both the forward and inverse problems and 3: Physical properties of the soil obtained by inverse problem using ∆T = ∆T optimum using the calculated values in 4.2 as initial guesses. The residue is 7.928.

	Table 4.2: Values of soil parameters and residue obtained by varying become less accurate as ∆T kh increases.	∆T kh . Results
	X X X X X X X kh X X ∆T Soil parameter X X	(ρC) s (J/m 3 K) λ s (W/mK)	φ	residue
	1		Failed	Failed	Failed	Failed
	2		Failed	Failed	Failed	Failed
	3	1.149 × 10 6	0.5362	0.2319	12.146
	4	8.654 × 10 5	0.4748	0.2436	14.602
	5	6.42 × 10 5	0.4306	0.2526	17.122
	7	3.351 × 10 5	0.3770	0.2640	22.876
	10	6.159 × 10 4	0.3400	0.2724	33.845
	11		Failed	Failed	Failed	Failed
	exact		1.95 × 10 6	0.756		0.20
	initial guess	1.149 × 10 6	0.5362	0.2319
	calculated	1.939 × 10 6	0.733	0.2001
	4.3 Levenberg Marquardt Algorithm (1D)

toward acceptable values if we use huge mesh size to generate the synthetic data and 120 mesh cells in the inverse problem (check table 4.3). We call this special technique in solving the inverse problem using dierent values of ∆T over two or more steps to insure convergence: chaining of the inverse problem. ρC s (J/m 3 K) λ s (W/mK) φ

In section 4.2, we explained that our inverse problem can be viewed as a nonlinear least square minimization problem which is solved by the Damped Gauss Newton Algorithm.

In this section, we present a more robust algorithm to solve the nonlinear least square minimization problem known as Levenberg Marquardt Algorithm (LMA). LMA is the most widely used optimization algorithm for the solution of nonlinear least square problems. It outperforms simple gradient descent and other conjugate gradient methods in a wide variety of problems. It is a blend of original gradient descent and Damped Gauss Newton iteration.

Table 4 .

 4 4: Physical properties of the soil obtained by inverse problem using LMA.Scaling is used implicitly and approximation (4.26) is removed. Same number of mesh cells is used in both the forward and inverse problems (321 mesh cells).

		(ρC) s (J/m 3 K) λ s (W/mK)	φ
	exact	1.95 × 10 6	0.756	0.20
	initial guess	2.0 × 10 6	0.8	0.18
	calculated	1.957 × 10 6	0.758	0.1996
	In gures 4.5, 4.6 and 4.7 we notice that the calculated values are very close to

the exact ones due to the fact that we used same number of mesh cells in both direct and inverse problems. Moreover, we notice that convergence is obtained after about 10 iterations.

Table 4 .

 4 5: Physical properties of the soil obtained by inverse problem (5 iterations).

		(ρC) s (J/m 3 K) λ s (W/mK)	φ
	exact	1.95 × 10 6	0.756	0.20
	initial guess	2 × 10 6	0.8	0.18
	calculated	1.9499 × 10 6	0.75596	0.200003

Table 4 .

 4 6: Physical properties of the soil obtained by inverse problem (6 iterations).

Table 4 .

 4 7: Physical properties of the soil obtained by inverse problem using dierent number of mesh cells where the synthetic data are generated by 150,000 mesh cells.

	Mesh	(ρC) s (J/m 3 K) λ s (W/mK)	φ	Residue
	30 × 50	1.23 × 10 6	0.403	0.232	463.9
	60 × 100	1.97 × 10 6	0.651	0.198	173.8
	90 × 150	2.05 × 10 6	0.712	0.195	97.2
	120 × 200	2.045 × 10 6	0.735	0.195	62.9
	150 × 250	1.95 × 10 6	0.725	0.199	45.6
	inverse problems. We notice that the 3D version of the inverse problem is more stable
	than the 1D version.

Table 4 .

 4 8: Comparison between the results obtained by 3D axisymmetric and 1D inverse problems. We vary the initial guess to notice the advantage of the 3D conguration.

	exact	1.95 × 10 6	0.756	0.2
	initial guess 1	1.0 × 10 6	1.1	0.1
	calculated by 1D	Fails to converge using ∆T optimum
	calculated by 3D	1.949 × 10 6	0.7559	0.2
	initial guess 2	1.5 × 10 6	0.35	0.3
	calculated by 1D	Fails to converge using ∆T optimum
	calculated by 3D	1.949 × 10 6	0.7559	0.2

(ρC) s (J/m 3 K) λ s (W/mK) φ

Table 4 .

 4 9: Physical properties of the soil obtained by inverse problem taking into consideration the noised measures in 1D. The same number of mesh cells is used in both the forward (to generate synthetic data) and inverse problems.

		(ρC) s (J/m 3 K) λ s (W/mK)	φ
	exact	1.95 × 10 6	0.756	0.20
	initial guess	2 × 10 6	0.8	0.18
	calculated	2.002 × 10 6	0.790	0.188

Table 4 .

 4 10: Physical properties of the soil obtained by inverse problem taking into consideration the noised measures in 3D. The same number of mesh cells is used in both the forward (to generate synthetic data) and inverse problems.

	exact	1.95 × 10 6	0.756	0.20
	initial guess	2 × 10 6	0.8	0.18
	calculated	1.99 × 10 6	0.782	0.198

(ρC) s (J/m 3 K) λ s (W/mK) φ

Table 5 .

 5 1: Physical properties of the soil obtained by inverse problem using real exper-

	imental data.

  Diusivity α in m 2 /s 3.21 × 10 -7 1.92 × 10 -7 2.77 × 10 -7

	Sensors	1-2	2-3	7-8
	Table 5.4: The values of diusivity (Fontainebleau Sand) obtained by Laplace transform
	method using the rst experiment.
	Sensors	10-14	7-11	11-13

Table 5 .

 5 The two experiments studied are experiments number 13 and 12 performed by José Augustin Cordero in the archeology lab of Rennes CREAAH. Both experiments are performed in the same manner as described in 5.2.1 using same box, heating plate and Fontainebleau sand but dierent sensors positions. The value of porosity (measured in the laboratory)for both experiments is φ = 0.4. Using experiment 13 data[START_REF] Cordero | Transferts de chaleur dans un sol sec chaué par le dessus[END_REF], we obtained the results summarized in table 5.6. 10 -5 6.581 × 10 3 120 × 200 3.203 × 10 -5 6.584 × 10 3 300 × 500 3.060 × 10 -5 6.582 × 10 3 6: The values of diusivity (Fontainebleau Sand) obtained by Inverse Problem using dierent mesh sizes.

	Mesh	α (m 2 /s)	Residue
	30 × 50	3.909 ×

table 5

 5

	.7
		Uniform	Non-uniform
	α(m 2 /s)	2.1 × 10 -7	2.4 × 10 -7
	shif t r (mm)	-0.14	-6.64
	shif t z (mm)	4.73	1.87
	Standard deviation for δ r (mm)	6.58	1.97
	Standard deviation for δ z (mm)	5.44	4.44
	Residue	816.5	787.6

Table 5 .

 5 7: Results obtained using the second experiment.The new value of diusivity is α = 2.4 × 10 -7 m 2 /s which is acceptable, even if two dierent values for the diusivity have been obtained from experiments 12 and 13.

	Indeed, it is not uncommon to get some variations from one experiment to another,
	which may come from the ambient humidity of the room. In this second experiment,

we are condent of the numerical results because we obtain diusivity values close to the ones obtained by previously presented methods. The experimental (old) positions of sensors, the new numerical positions and the direction of displacement are presented in gure 5.27.

  x, T, P ) + mass ux at interface (gaseous zone)

	∂P ∂x	= 0 (liquid zone)
		∂n w ∂t	= h(t, x, T, P, n w ) only for 0 < x < ξ
			dξ dt	= q(t, P, n w )
			(6.22)

Acknowledgement

shared hard and good moments that we will never forget, you were next to me in all the ups and downs. Areej, simply, you are my unbiological sister. I love you... My family who were always the big support: my sisters, brother, brother in-law and my small Amar. Having you in my life is a bless.

Finally, Mom and Dad, thank you for your unconditional love and care. Believe me this work couldn't be without your support.

You are my true love! Résumé Détaillé

chapter 5

Reformulation of our inverse problem

The unknowns of our inverse problem are:

The diusivity α.

The displacement in both r-direction and z-direction for each sensor i (δ r i and δ z i ).

The bias in both r and z directions (shif t r and shif t z ).

The constraints applied to our inverse problem are: n i=1 δ r i = 0 and n i=1 δ z i = 0.

For each sensor, sin(φ i ) δ r i = cos(φ i ) δ z i .

Box constraints over the new sensors' position to ensure that they are in the physical domain under study.

Box constraint over the value of α to assure that it will attain a positive and reasonable value.

Numerical Strategy

The obtained system of coupled equations (heat diusion equation 5.14 + sensitivity equation 5.17) is a system of partial dierential equations. To solve this system, we use the same numerical strategy explained earlier (method of lines + Finite Volume method). After spatial discretization, the system of coupled equations can be written in the form of a system of rst order implicit ODEs:

where Y = [T, U α ]. The system in equation (5.22) can be solved by an ODE solver as in 4.1.1.

The ODE solver will provide us with the values of T and U α which are not enough to solve our inverse problem that have the sensors positions as unknowns. To approximate U r i = ∂T ∂r and U z i = ∂T ∂z at each sensor i, we use an interpolation of order 2 based on a bi-quadratic interpolation (9-point stencil) because a rst trial using interpolation of order one didn't gave accurate results.

Solving the constrained inverse problem

The Levenberg

Abstract

In this work, we are interested in studying heat and mass transfer in water saturated and unsaturated porous medium with a strong heating at the surface. Applications concerned are archaeology, agriculture and geothermal engineering.

The rst part of this work concerns the improvement of the AHC (Apparent Heat Capacity) method used in the numerical resolution of phase change problem in a homogeneous medium: the phase change temperature interval, over which the heat capacity varies, appears as a key parameter which must be chosen proportional to the mesh size. The third part presents a simple model to calculate the eective thermal conductivity of a granular medium which contains a small quantity of liquid water. The exact shape of the liquid menisci between the grains is calculated at equilibrium. The eective thermal conductivity experiences a hysteresis behavior with respect to the liquid volume. A future work that concerns a new unsaturated model, restricted to the pendular regime and detailed at the end of this thesis, should be able to use this result.