K Wolski Directeur

S Drapier

Genie Des

Maître Procedes F Gruy

Appliquees O Mathematiques

Maître-Assistant Roustant

O Informatique

Professeur Boissier

Professeur J C Pinoli

Genie Industriel

Nabil Absi

Maitre De

Professeur Ph Lalevée

Image, Vision, Signal CIS Vincent Cr

Tria Assia

Valdivieso François

DR Jean Paul

Wolski Krzy

Xie Xiaolan

Yugma Gallian

Romain Noël

J-C Afaf

Nico Fils

Sophie Joseph

Ozge Merci

Merci À Tetaz

Abel Foufouch

Pierre R Aura

Les Alix

Pauline & Millot Ainsi Que Afaf

Édouard

Cédric Sandy

Romain & Marjorie

Nicolas Durrande

Pauline Julou

Chaillou

Ily & Patrick

Audrey & Fouquet

Nicolas Alban

Pauline Valance

Sylvain & Mathilde

Fils, Clément Pierre G Laurène

Claire & Marie

Tété Joël

Spécialités doctorales Responsables : SCIENCES ET GENIE DES MATERIAUX

To determine the evolution in time of existing bodies, it is necessary to start from a given confi and perform predictive computations. To obtain information from within the matter without destroying it, is not an easy task. When degradation is not possible, a common way to acquire the initial configuration is the use of imaging systems. And due to the complexity of the behaviour of matter, it is also common to resort to numerical methods in order to simulate these temporal processes. A relatively recent numerical method called Lattice Boltzmann Method (LBM) tackles the evolution process at a mesoscopic scale. The combination of these two fields is originally investigated in this thesis, around the use of the LBM for numerical simulations of continuum medium aiming image-based diagnostics.

Foremost, an extension of the LBM to the morphological mathematics is suggested. It contributes to the development of a new concept: the coupling of image processing, mechanical and biological simulations on a single network. Then, the simulation of solids with LBM is addressed, using two different approaches. The first one provides an analytical Multiple Relaxation Times methods to generate arbitrary stress tensor and heat flux. The second introduces the divergence of the stress tensor into the Vlasov equation. Both approaches are numerically confronted with theoretical results in 1D and 2D and offer promising perspectives. 1 0 nous serions bien démunis. Je voudrais aussi remercier David P le courageux entrepreneur, et Miquel pour leurs conseils ; ainsi que Victor y Olfa los cálidos compañeros con los que puedo intentar hablar español.

J'en profite pour dire ici, merci à Hélène Navarro pour son soutien. Je remercie également Paolo de partager sa bonne humeur et David G de m'avoir « voler mon bureau ».

Il me faut maintenant dire merci à tous les doctorants avec qui j'ai partagé pas mal de souvenirs. Dans les plus anciens nous avons Armelle la footballeuse, Witold le Lyonnais, puis l'incontournable Fanette, Boris le Dude du dudisme, Pierre C le John Malkovich du ping et grand guide de Toronto, Cristina la référence en matière de tiramisu même si je n'ai jamais réussi à en manger, Thomas l'éternel compagnon du 19h12, et mes acolytes de l'ASEC Marc le master de TowerFall, Charles un des rares à ne pas être de l'EMSE, et Piou-Piou le danseur endiablé. Il me faut aussi remercier ceux qui sont dans une période difficile, la troisième année: Aymeric le Machopeur canard, Radha qui n'a pas peur de partir en grande voie improvisée, Mohsen, Yoann l'arracheur de bouton de chemise, Sophie la schtroumpf hibou bretonne, Fei l'expert du volume rendering. Vient ensuite le bureau des quatre fantastiques (je viens de l'inventer et je vous laisse vous battre pour savoir qui est qui) : Nicolas le papa éleveur de moutons du Yosemite, Joseph le Savoyard qui ne grimpe plus beaucoup (un comble), Joan la cabeza loca de las montañas chilenas, et Claudie qui mériterait de venir plus souvent aux pauses car sa conversation est toujours intéressante. Sans oublier la coloc des princesses, la Dráma queen Ozge qui m'a sauvé au Bul, et Cyriac la diva. Il me faut aussi exprimer ma gratitude la nouvelle relève Bastien le trailer sans barbe, Jules le mytho-breton nanto-parisien et Zeineb la discrète, mais effi relève en traitement d'image. Je n'oublie pas non plus les stagiaires plus ou moins récents, Rosaline Basseterre, Alexis le futur programmeur et Baptiste le nouveau responsable de la plante du bureau. Avant de continuer je voudrais en profiter pour exprimer ma gratitude toutes les personnes qui mon donner le gite et le couvert à Saint-Étienne au cours de ces

Résumé

Pour déterminer l'évolution temporelle de corps existants, il est nécessaire de partir d'une confi donnée et d'effectuer des calculs prédictifs. Obtenir des informations de l'intérieur de la matière sans la détruire, n'est pas une tâche facile. Lorsque la dégradation n'est pas possible, un moyen courant d'acquérir la confi de départ est l'utilisation de systèmes d'imageries. Et en raison de la complexité du comportement de la matière, il est également courant de recourir à des méthodes numériques pour simuler les processus temporels. Une méthode numérique relativement récente appelée méthode de Boltzmann sur réseau (LBM) permet d'aborder le processus d'évolution à une échelle mésoscopique. La combinaison de ces deux domaines est étudiée de manière originale dans cette thèse, autour de l'utilisation de la LBM pour les simulations numériques des milieux continus en vue de diagnostics à partir d'images.

Tout d'abord, une extension de la LBM à la morphologie mathématique est proposée. Elle concourt au développement d'un nouveau concept: le couplage du traitement d'images, des simulations mécaniques et biologiques sur un seul réseau. Ensuite, la simulation des solides avec la LBM est abordée, via deux approches différentes. La première fournit une méthode à multiples temps de relaxation analytique pour générer un tenseur de contrainte et un flux de chaleur arbitraires. La seconde introduit la divergence du tenseur de contrainte dans l'équation de Vlasov. Les deux approches sont confrontées numériquement à des résultats théoriques en 1D et 2D et offrent des perspectives prometteuses.

Dédicace

To my lost family members and friends.

Preface

This document is fulfi ment to the grade of Philosophiae Doctor in science. It relate s the research works done between May 2016 and May 2019. In order to help the readi ng, the symbols and notations adopted in this document are listed below. (D) volume measure m D dx (D-1) surface measure m D -1 q generalized momentum m.s -1 q θ heat flux W.m

Abbreviations

General introduction

Context & Generalities

A Danish proverb, mixed by physicist Niels Bohr and other authors, explains with simplicity and pleonasm the idea that "forecasting is always difficult, especially when it concerns the future". That is where one of the great science benefi stands, trying to understand surrounding things, allowing to foresee their future evolution. The present thesis manuscript follows this idea. In fact, in multiple cases it is useful, if not necessary, to anticipate the behaviour an object or a material is susceptible to adopt. In order to give answers, when the trajectory of an object like a planet is at stake, the shape and the internal actions do not really matter; and so, Isaac Newton's science resorts to an equations system and fundamental principles of mechanics. The latter enable to link an object movement and the forces it experiences. In simpler words: an action implies a reaction. But when the light is put on the intern actions and behaviours of the said object, it is necessary to use a finer mechanic that has interest in the continuum medium. This mechanic has to link forces or constraints to the deformation of a material containing no discontinuity.

Thus, to be able to describe the material, continuum medium mechanics needs to obtain constitutive laws linked to each individual material. These behaviour laws are the macroscopic reflect (visible by the human eye and quite generally about a millimetre or centimetre magnitude order) of microscopic interactions between elements in the materi al (molecular aggregate of a tens of nanometres magnitude order). If it is possible to see the se laws from an only macroscopic point of view, it is also very useful to tackle them from a microscopic point of view, since it purveys a finer information on the mechanical state of the material.

The complexity of the microscopic point of view stands in the number of elements or particles (as elementary objects of small size) to be taken into account. So, statistical physics and Ludwig Boltzmann approach allow to work with statistical distributions of elements rather than with a large number of elements. One then talks about mesoscopic point of view. This roughly simplifies the problem as it averages the interactions betwe e n the elements. Thereby, the Boltzmann equation allows to describe the evolution of a distribution of particles submitted to certain actions.

Nevertheless, the complexity of realistic situations requires quite often the nume ri cal tool to give an answer or a satisfying prediction. The terms numerical simulation is the n employed, in order to name the use of a computer program allowing multiple calculations based on a mathematical processes and models. When only the macroscopic continuous aspects are studied, the numerical processes like finite differences, volumes or elements are commonly used. On the other hand, numerical simulations interested in mesoscopic systems governed by the Boltzmann equation have benefi from using the LBM. The latter method is blooming nowadays.

Numerical simulations for direct problems (id est interested in predicting the final state of a system knowing its initial state and the actions undergoing) are distinguished into two main categories of uses. The first is for design: when designing a new vehicle, it is essential to ensure the good resistance of its mechanical parts composing. The second category is for diagnostic: when working on an old building, it must be verifi that the building can support them. In this second case, the accuracy of the results depends on the quality of the information provided by the sensors making it possible to measure the properties which will then be injected into the numerical simulations. So, to an equivalent model, if the first category has all the information necessary for computations during design (at least theoretically, since they are supposed to be desired), then the resulting variability of it lies in the precision of the machines shaping mechanical parts; the second has the maximum amount of information provided by the sensors and must compensate for the missing information.

The context of numerical simulations for biomechanical diagnostic is really representative of the second category. It is realistic to consider that the sensor system is limited to imaging systems because the direct sampling is avoided at maximum. So, the natural progression to realize numerical simulations for diagnostic is rather expensive in time, skil l s and money, because it requires successively:

1. to obtain images through a suitable acquisition device; 2. to process images so that only relevant information is kept for the rest of the procedure;

3. the use of classical continuous methods on meshes involve reconstructing a 3D geometric model based on the processed images;

4. build from the geometric model or images a mesh adapted to the computations that will be required;

5. build a mechano-numerical model from the mesh and properties observed or conjectured, allowing the numerical simulation of the studied problem.

Each of these steps involves a potential loss of information. Steps 2, 3 and 4 also requi re scientists to make important assumptions and decisions. A remarkable poi nt i s that in the vast majority, once the computations are made, the results available on the mesh are projected on the screen of the scientist microcomputer. Since this screen is nothing more than an image displayer, it means that the scientist has started from an image taken by an acquisition system to end with an image used for the interpretation and analysis of the results, to be able to predict. Moreover, the definition of the acquisition system is the maximum precision available without inputting a priori knowledge about the image anatomical structures. The choice of the method for biomechanical simulation and diagnostic from image s, is a question not to neglect. Each of these steps also requires time and resources: compute r equipment, electricity... Also, the reduction of one or more steps could bring a considerable gain to society. This would make the diagnostic more eff e and faster, either for a building that has suffered a strong earthquake or for a diagnostic of a patient requiring surgery.

A solution to remove step 4 is to avoid to work with standard numerical methods, but with meshless methods. Among the most famous, one can cite the Smoothed Particle Hydrodynamics or the Material Point Method which work on macroscopic particles. Their effi and performances have been proved and they allow nowadays to work on all types of continuum medium (solid, liquid and gas), even if they require more computational time. However, unlike mesh methods, they are not well adapted to perform image processing.

The LBM is also a "meshless" method, with the asset of allowing numerical simulation at a mesoscopic scale. This property is mainly due to the fact that, contrary to other meshless methods, it is based on statistical mechanics. In addition, since the beginning of the 2010s, the method has been applied to image processing with bolstering results. This means that with such a method it is possible to free oneself from step 4, and to carry out in step 2 with the same computational method. This simplifies the whole procedure. Furthermore, the LBM is a computational method known for its computational speed. This speed lies in its explicit and local character, which also makes it highly parallelizable.

Despite the important attractions of the LBM in terms of simulation for image-based diagnostic, it does not only have advantages. In particular, its mesoscopic nature makes the work on solid materials or high viscosity fluids difficult. Indeed, only few works on the capacity of the LBM have been conducted, and those works mention some locks. Thus, the use of such a method for the diagnostic from image is not something to be taken for granted because it still requires many developments.

Positioning the problem

Numerical simulations for image-based diagnostics is resource-intensive, but is a source of considerable progress for society, considering the increasing importance of images. Si nce , the LBM responds to many criteria that the numerical simulation for image-based diagnostic requires, the choice of relying on this method is done. This choice is motivate d by the ambition to bring missing elements to the method. So, the objectives of thi s the si s are to reduce the gap between the LBM and the image-based diagnostic. Indeed, se ve ral scientific locks persist to create such practical approach.

Although it is possible to work on images with the LBM, many image processing techniques have never been transposed nor discovered. This makes the LBM a perfectible and improvable tool to implement image processing algorithms. Moreover, the numerical simulation of fluid flows directly from an image represents an unprecedented experiment.

The LBM is known to be effi t when working with liquids and gazes. However, when it comes to solid matter, the LBM seems to encounter difficulties and has to be adapted. The desire to be able to study all types of continuum medium means overcoming the challenge of solids for LBM.

Also, this thesis suggests studying the use of the LBM for the simulation of continuum medium for image-based diagnostics. This subject is a scientific challenge for both the theory and the applications within its reach. With the objective of making a complete tool for the simulation of diagnostics from images with LBM, this thesis aims to obtai n a complete pipeline from the image to the simulation, using only the LBM as a computati onal tool.

In order to reach this objective, this manuscript is decomposed in several chapters. The first chapter recovers the context and State of the Art around the subject to exhibit the locks and opportunities which it offers. Since the subject aims to push forward the LBM, the second chapter reviews an implementation and the validation of our own numerical tools that serves as the foundation for the improvements to come. The third chapter explains how it is possible to use images directly as input of a LBM simulation. It also evokes the LBM-based image processing works developed for our purpose. Then, the fourth chapter aims to obtain an arbitrary internal stress with the LBM. This is, of course, motivated by the idea to mimic solid behaviour. To finish, the fifth chapter tackle the solid problem through a different angle: via the Vlasov equation. The manuscript ends with some conclusions and perspectives brought by these research works. The interested reader can also find in the first appendix, an additional work about the use of the LBM to address heterogeneous and multi-class traffic questions.

Abstract of the chapter

The LBM is an original method which presents some relevant aspects for image-based diagnostic. Since, the method is relatively younger than the conventional methods for continuum mechanics, it requires more developments to be fully applicable to image-based diagnostic. In particular, as it will be emphasized in this chapter, eff in image processing and computational solid dynamics have to be done with the LBM.

With a top-down approach, this chapter starts looking at the macroscopic scale. This allows to differentiate the Lagrangian from the Eulerian point of view. This difference is also felt in the numerical methods evoked after. As a link between the macroscopic and microscopic scale, the Boltzmann equation starts from the particle dynamics to end up with the Navier-Stokes-Fourier equations. The Boltzmann equation can be solved numerically with the LBM. As many numerical methods, the LBM implies multiple numerical aspects and developments briefly summed up here.

Résumé du chapitre

La LBM est une méthode originale qui présente quelques aspects pertinents pour le diagnostic à partir images. Étant donné que la méthode est relativement plus jeune que les méthodes conventionnelles de la mécanique des milieux continus, elle nécessite pl us de développements pour être pleinement applicable au diagnostic à partir d'images. En particulier, comme il sera souligné dans ce chapitre, des eff en traitement d'i mage s et en mécanique computationnelle des solide doivent être faits avec la LBM.

Ce chapitre adopte d'une approche partant du macroscopique vers le microscopique, Le premier permet de différencier le point de vue lagrangien de l'eulérien. Cette différence se ressent également dans les méthodes numériques évoquées par la suite. Puis, pour faire le lien entre l'échelle macroscopique et microscopique, l'équation Boltzmann part de la dynamique des particules pour arriver aux équations Navier-Stokes-Fourier. L'équation Boltzmann peut être résolue numériquement avec la LBM. Comme beaucoup de méthodes numériques, la LBM implique de multiples aspects numériques et développements qui sont brièvement résumés ici.

I.1 Introduction

To make this manuscript more comprehensible, the main foundations and concepts relate d to the subject are developed in this chapter. It is the opportunity to introduce the equati ons necessary to understand the ideas presented in this thesis.

In order to highlight the reasons and choices around the subject of the lattice Boltzmann method for numerical simulation of continuum medium aiming image-based diagnostics, it is necessary to look closer at equations and current progresses in the domain. This chapter comes back on the great physical principles of the macroscopic mechanics. This allows to understand the main numerical methods to solve this mechanics. Then, a brief review of the Boltzmann equation enables to understand the concepts and developments linking this mesoscopic equation with the macroscopic Navier-Stokes-Fourier equations. This link is also very important to evaluate the advantages of the numerical method based on this equation. So, from this equation the construction of the famous lattice Boltzmann method can be done. Once this method is introduced, numerical aspects and developments can be detailed, which leads naturally to unsolved questions.

I.2 Continuum mechanics

I.2.1 Macroscopic Lagrangian mechanics

To describe the dynamics of a macroscopic physical system, one must look at all the physical principles governing its evolution. At first look, these principles are:

• conservation of the mass,

• conservation of the quantity of the movement. The conservation of the mass, also called the Lavoisier-Lomonosov principle named after the scientists having demonstrated it [START_REF] De Lavoisier | Traité élémentaire de chimie[END_REF], stipulates that during a mechanical M dW or chemical transformation, the mass of the overall system does not vary with time. This result is translated in the following equation:

dmT dt = Sm = 0 (I.1)
where mT represents the system and Sm is the mass source term. The Fundamental Principe of the Dynamics (FPD) [START_REF] Newton | Principes mathématiques de la Philosophie naturelle[END_REF] of a body assimilated to its centre of mass reads:

dmT V dt = , Fi i=1 (I.2)
with Fi the forces applied to the system which the speed of its centre of mass is V . This equation expresses the idea that a variation of momentum (or mass times the acceleration) is equal to the sum of the external forces stressing the system. This FPD must be completed with the possibility for the system to start rotating. This leads to the equation of Euler-Newton:

dL = , dt i i=1 (I.3)
where L is the angular momentum and Mi is the torque associated to the i-th force. This means that the variation of angular momentum (or angular acceleration) is equal to the sum of the external torque experienced by the system.

With these two equations, one can describe the movement of planets. But it is harder to get information about the internal stress of a compressed spring.

In order to start looking at the internal solicitations of a system it is possible to exploit the two thermodynamics principles. The first is interested in the conservation of the total energy of the system, allowing conversion between energies eventually, and reads:

dEint = dQ + dW (I.4) dt dt dt
with Eint is the internal energy of the system, Q is the thermal energy of the system and dt is the power worked by external forces not used to generate movement. The second thermodynami cs principle introduced by Carnot [Car24], translates the irreversibility of the physical phenomenon, which leads to the following equation:

dS dt ≥ 0 (I.5)
where S is the entropy of the system, characterizing the disorder can only increase.

I.2.2 Eulerian mechanics

The previous principles and equations are expressed with a global and Lagrangian point of view, i.e. by following the material points in space. To obtain a local Eulerian formulation, it is possible to use the Leibniz-Reynolds transport theorem [START_REF] Reynolds | Papers on Mechanical and Physical Subjects: The Sub-Mechanics of the Universe[END_REF] which can be written:

d r dt V(t) A dx (D) = r V(t) ∂A dx (D) + ∂t f ∂V(t) ( v b • n A dx (D-1) (I.6)
where A is an extensive physical quantity which is integrated over the volume V, its-se l f depending on the time. The derivative ∂V represents the surface delimiting this volume, v b is the speed of the boundary surface points and n is the out pointing normal associated. By using, the Green-Ostrogradsky divergence theorem [Gre28], it is possible to transform the surface integral into a volume integral and so: f

∂V(t) (v ⊗ A) • n dx (D-1) = r V(t) ∇x • (v ⊗ A) dx (D) (I.7)
The application of the previous theorem allows to write the mass conservation in an Eulerian formulation through the form: r

r ∂ρ mT = r V(t) ρ dx (D) (I.8) r V(t) ∂t dx (D) + V(t) ∇x • (v ⊗ ρ) dx (D) = V(t)
Sm dx (D) (I.9)

where ρ is the volume mass density of the medium. As now, all the integration elements are equal, they can be removed. The expression of the mass conversation in its local Eulerian formulation, also called the continuity equation, reads:

∂ρ ∂t + ∇x • (v ⊗ ρ) = Sm = 0 (I.10)
The local Eulerian formulation of the momentum conservation from the FPD can be obtained in the same manner and is called the first undefined equation of (linear) motion of Cauchy [START_REF] Cauchy | Mémoire sur l'équilibre et le mouvement des liquides et des fluides élasitque[END_REF]:

, Fi = mT g (I.11) i ∂ρv ∂t + ∇x • (ρv ⊗ v) = ∇x • (σ) + ρg (I.12)
where σ represents the internal stress modelled by the Cauchy (second order) stress tensor.

Writing the equation of Euler-Newton gives rise to the second undefined equation of Cauchy (angular) motion:

, Mi = C (I.13) i ∇x ∧ (σ) = C (I.14)
with C the volumetric torque undergone by the material. This latter field is generally equal to zero in the absence of a magnetic field, which implies symmetry of the stress tensor of Cauchy. This assumption will be used throughout all this document and summarized by the equations

∇x ∧ (σ) = 0 ⇐⇒ σ T = σ (I.15a) ∂ρv ∂t + ∇x • (ρv ⊗ v) =∇x • (σ) + ρg (I.15b)
The scalar multiplication of the first equation of Cauchy with a virtual velocity (or real velocity) allows to obtain the local formulation of the virtual work principle [Ber25]:

(∂ t (ρv)) • v * + (∇x • (ρv ⊗ v)) • v * -(∇x • σ) • v * -ρg • v * = (∂ t (ρv)) • v * + (∇x • (ρv ⊗ v)) • v * -∇x • (σ • v * ) + σ : (∇x • v * ) -ρg • v * = P k -Pc + σ : (∇x • v * ) -Pg = 0 (I.

-

The integral formulation of the last equation allows the writing of the variational formulation of the problem. This latter weak form is used to build the Finite Element Method (FEM). In addition, this last equation identifies the term σ : (∇x • v * ) with the worked power as defined above (see eq. (I.4)). After injecting this equation into the first principl e of thermodynamics one can write the conservation of total energy:

( v 2 \ ∂ t ρ + Eint 2 + ∇x • (r v 2 l ρ + Eint 2 \ v -σ • v + q θ + r + ρg • v = 0 (I.17)
Finally, using the inequality of Clausius-Duhem [CD67], the writing of the second principle of thermodynamics in local Eulerian formulation gives:

r ∂ t (ρs) + ∇x • (ρsv) ≥ θ -∇x • ( q θ θ (I.18)
where r is the thermal radiation and q θ is the heat flux, s is the specifi entropy and θ is the temperature.

In summary, the evolution of the state of matter is translated by that of the macroscopic variables ρ, v, σ, θ, id est in 3 dimensions of 14 independent scalar variables or more generally in D-dimension of D 2 + D + 2 independent variables. However, the available equations are 1 + (D + D 2 -D ) + 1 + 1, that is:

• mass conservation (see eq. (I.10)),

• undefined equations of motion (see eq. (I.15)),

• first principle of thermodynamics (see eq. (I.17)),

• second principle of thermodynamics (see eq. (I.18)).

I.2.3 Closing relations

It is then necessary to have recourse to D 2 +D 1 supplementary equations to describe the evolution of matter, these are the constitutive laws describing the links that may exist between stress, density, velocity and temperature.

Classic examples of behavioural relationships are for a Newtonian viscous fluid, the stress is linked to the gradient of speed by the following Stokes linear stress constitutive equation

1 σ = ν 2 ( ∇xv + ∇xv T ) . (I.19)
When combined with the Fourier law for heat conduction, given by

q θ = κ∇x (θ) (I.20)
one would obtain the famous system:

∂ t (ρ) + ∇x • (ρv) = 0 (I.21a) ∂ t (ρv) + ∇x • (ρv ⊗ v + p -ν∇x v) = 0 (I.21b) ∂ t (Etot) + ∇x • (Etotv + p • v + κ∇x θ) = 0 (I.21c)
with ν the cinematic viscosity and κ the thermal conductivity of the fluid. The system is called the Navier-Stokes-Fourier equations [START_REF] Navier | Mémoires sur les lois du mouvement des fluides[END_REF].

Another classic example is homogeneous isotropic linear elastic solids, where the relation eq. (I.19) on the facing page is replaced by the Hooke law

σ = C : ε = GJ -5/3 ( Tr (B) \ B -3 I + K (J -1) I (I.22) dx ∂ t (u(X, t)) = v F = I + ∇X (u) = dX (I.23) B = FF T ≈ I + 2ε J = / det (B) (I.24) χ (x, t) = X = x -u(X, t) ∇x (χ) = dX = F - dx ∂ t (χ) + v∇x (χ) = 0 (I.25)
where G and K are respectively the bulk and shear modulus. This leads to the following system

∂ t (ρ) + ∇x • (ρv) = 0 (I.26a) ∂ t (ρv) + ∇x • (ρv ⊗ v + p -C : ε) = 0 (I.26b) ∂ t (Etot) + ∇x • (Etotv + σ • v + κ∇x θ) = 0 (I.26c)
These different behaviours are explained by the different interactions between the molecules that make up matter.

I.2.4 Macroscopic numerical methods for continuum mechanics

The differences of matter behaviour imply differences in the numerical methods employed to simulate them. Here, these methods are divided in two groups, those to simulate fluid and those to simulate solids. Both of these groups are subdivided in conventional methods on mesh and meshless methods.

I.2.4.1 Computational Fluid Dynamics

I.2.4.1.a Methods on mesh

It exists several numerical methods to solve the Partial Differential Equations (PDE) related to fluid dynamics. The conventional macroscopic methods to solve these equations require a mesh as a support for the computations. Sorted by order of complexity and precision these methods are:

• Finite Difference Method (FDM): This method uses a Taylor expansion to estimate the derivatives and solve numerically the PDE. It results in a great advantage of the method, its simplicity. However, when it comes to Computational Fluid Dynamics (CFD) some weaknesses appear. In a general case, it is not a conservative method, some errors are introduced in the mass, momentum and energy conservations. Also, this method has issues with complex geometries [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems[END_REF].

• Finite Volume Method (FVM): This method is intended to solve conservative equations. It integrates the conservation equation over a small volume and use the divergence theorem to evaluates all terms (except sources and sinks) as flux crossing the volume boundary surfaces. So, naturally this method has the advantages to be conservative. It is also easier to work with complex geometries. It is probably the most used method for CFD today. However, it is not straightforward to find an appropriate unstructured mesh compatible with the finite volume method, especially when higher orders are required [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method[END_REF].

1 2

• Finite Element Method (FEM): This method solves the PDE by using Galerkin methods of discretization on a weak form and the Lax-Milgram theorem proves the existence and uniqueness of the solution. The method has strong enough mathematical roots to work on unstructured mesh even with higher orders. It has the drawbacks not to be conservative, and be more complex than the two previous methods. The FEM may also suffer from instabilities which kept it a bit away from the CFD for years [START_REF] Zienkiewicz | -Clefs: Méthode de Boltzmann sur réseau, diagnostic à partir d'images, physique statistique des solides, equation de Vlasov, distribution normale multivariante asymétrique[END_REF].

I.2.4.1.b Meshless methods

By opposition with the conventional methods on mesh, the meshless methods do not require such support to perform their computations. This absence of mesh is a major advantage in numerical simulation, since many complexities arise from this mesh [START_REF] Li | Meshfree and Particle Methods and Their Applications[END_REF]. At a macroscopic scale, two main methods dominate this field, the Smoothed Particle Hydrodynamics (SPH) and the Material Point Method (MPM):

• Smoothed Particle Hydrodynamics (SPH): uses macroscopic particles and determines an area of influence around one particle. In this area of influence the interaction between the particles are summed to determine the movement of one particle. In a general case, this is called a Verlet algorithm. When the interaction forces between macroscopic particles are taken with a hydrodynamics nature, it leads to the SPH method. The SPH is used for many hydrodynamic applications. It conserves the mass and momentum quantities [START_REF] Violeau | Fluid Mechanics and the SPH Method: Theory and Applications[END_REF]. The disadvantages of this method are the boundary conditions and the weak mathematical proof of the simulated equations.

• Material Point Method (MPM): is also based on particles but also on shape functions as the FEM. The basic idea of this method is that these material particles (points) are following the movement of their vicinity. So, from the movement of the material points the evolution of the macroscopic quantities is extrapolated to a grid. Then, the equation of motion is solved on the grid, and the resulting terms are extrapolated back to the material points before resetting the grid. The boundary conditions are simpler with the MPM than with the SPH. However, the MPM is more expansive in computational time since it uses both particles and grid data. It also suffers from particular sensitivity to mesh refinement [START_REF] Zhang | The Material Point Method -A Continuum-Based Particle Method for Extreme Loading Cases[END_REF].

I.2.4.2 Computational Solid Dynamics

The same numerical methods are also used for computational solid dynamics, even if the rates of use are not the same. Moreover, although the sources of difficulty are fundamentally the same for solids and fluids, solids are described with greater ease in a Lagrangian setting.

I.2.4.2.a Methods on mesh

• Finite Difference Method (FDM): The finite difference is not a usual method to solve solid dynamics equations. This is mainly caused by its difficulties to deal with complex geometries that offers solids. Yet some works succeed to the simulate solids with non-linearity and contacts using the FDM [START_REF] Valkov | Eulerian Method for Fluid-Structure Interaction and Submerged Solid-Solid Contact Problems[END_REF].

2

• Finite Volume Method (FVM): This method is more famous for CFD problems. However, it is also able to solve various solid problems like biomechanical behaviours, large deformation, fractures [Agu+15; CD18]...

• FEM: Due to the lack of convection and potential turbulence, the solid dynamics is more often study with a Lagrangian point of view, which seems more natural with the FEM. The FEM is the most used method for solids. It is able to address many complex problems such as contacts, non-linear constitutive equations, fractures, large deformations [BW97; Wri08]...

I.2.4.2.b Meshless methods

The meshless methods to solve computational solid dynamics are the same as those used for the CFD, however other methods like the lattice-springs come in complement.

• SPH: Even if the SPH method is originally designed for CFD it is also able to be extended to solve solid dynamics. Examples of textile tension or penetration tests using the SPH can be found. Yet, the simulation of solids slow transformation compared to the sound velocity remains challenging [START_REF] Hoover | Smooth Particle Applied Mechanics: The State of the Art[END_REF].

• MPM: Thanks to its combination of material points and grid extrapolation, the MPM is capable of solving complex problems in solid dynamics such as failure, shock waves, penetration [START_REF] Zhang | The Material Point Method -A Continuum-Based Particle Method for Extreme Loading Cases[END_REF]...

• Lattice-spring: it is a method assimilating the solid matter a regular network of mass and spring. The masses are responsible for the inertia while the springs surrogate the internal stress [START_REF] Pazdniakou | Lattice Spring Models[END_REF]. This method is simple and efficient. It is capable of dealing with heterogeneous plasticity or crack propagation, for example [START_REF] Buxton | A Lattice Spring Model of Heterogeneous Materials with Plasticity[END_REF]. Yet, simulation on large deformations remains challenging.

In addition to the latter, there exist several other numerical methods that go deeper into the matter to simulate the particles themselves.

I.2.5 Particles numerical methods

The evolution of a matter point particle which is in interaction with a second fixed particle in the study reference frame is given by the continuity equation and the FPD. The fundamental forces of interaction between particles are of various types: attractive and repellent. At these scales, the only attractive force is the gravitational interaction and the repulsive forces are those of electromagnetic interactions. The combination of these two forces causes an equilibrium position when the kinetic energy of the particle is not too high. The Lennard-Jones model [START_REF] Lennard-Jones | On the Determination of Molecular Fields. II. From the Equation of State of a Gas[END_REF] can be mentioned as an illustration of these interactions. This model yields the equation:

V LJ = ε LJ ( r m \ 12 ra ( r m \ 6 - r a (I.27)
with VLJ the potential energy between atoms, εLJ is the depth of the energetic sink and ra is the distance between the two atoms. The integration of such interatomic forces can be incorporated in numerical microscopic particles methods. These methods simulate the matter directly from an atomic point of view. At that scale, only Newton laws and the expression of the interatomic forces are necessary. Thus, they are using a reduced number of different PDE in their model, but there are repeated as many times as the number of particles. The two main methods at this scale are the Molecular Dynamics and the Direct Simulation Monte Carlo:

• Molecular Dynamics (MD): is a method that integrates over time with a discrete scheme positions and velocities of a large number of particles. So, using the Newton laws and the expression of the interatomic forces, the new positions and velocities are updated. Since, the interatomic forces are more important in a certain range , the Verlet algorithm can be exploited. This kind of method is very detailed and keep track of all particles, but it is too detailed at macroscopic scale. As a conseque nce, it requires too many resources for practical results [START_REF] Frenkel | Understanding Molecular Simulation from Algorithms to Applications[END_REF].

• Direct Simulation Monte Carlo (DSMC): is a method using the probabilistic Monte Carlo simulation to solve collisions model. It means that the interaction model is reduced to a collision relation. So, the DSMC, rather than looking for a solution of the Boltzmann equation, directly simulates the particles and repeats the simulation a large number of times to perform a statistical average. These simulations offer description of physical eff sometime beyond the Boltzmann formulation but are very time expensive [START_REF] Bird | Approach to Translational Equilibrium in a Rigid Sphere Gas[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF].

I.3.1 Billiards game with N balls

It is possible to obtain the equations of Navier-Stokes-Fourier which describe the evolution of the density, velocity and temperature of a viscous fluid, from the fundamental principle of the dynamics of Newton, as several authors recalled it [START_REF] Gallagher | From Newton to Boltzmann: Hard Spheres and Short-Range Potentials[END_REF]. To do so, one needs the use of the bridge which is partially made with the equation of Boltzmann.

The description of a system with N particles whose statistical distributions fN (pi, qi, t) are given by the Liouville equation [START_REF] Liouville | Note sur la Théorie de la Variation des constantes arbitraires[END_REF] with p the generalized coordinates, q the generalized momentum and HN the Hamiltonian of the system which reflects the energy of the system. The classifi on of the system given by BBGKY [BBGKY35] allows to switch to the equation of Boltzmann. Roughly, under some assumptions the Hamiltonian is divided in two parts, one part to take into account the global action of a force field like the gravity, and another part for the interatomic interactions.

So the Boltzmann Equation (BE) established in 1872 [START_REF] Boltzmann | Further Studies on the Thermal Equilibrium of Gas Molecules[END_REF] is

∂ f ˇ ∂t + ξ • ∇x ( f ˇ) + g • ∇ ξ ( f ˇ) = Ω ( f ˇ, f ˇ) (I.29)
where g is the force field felt by the particles, like gravity for mass particles. Ω is the collision-interaction operator representing the interactions between particles. In general, the collision operator is expressed with integral terms (replacing the sums in the Liouville equation). Therefore, the Boltzmann equation is an integro-differential equation describing the evolution of an out-of-equilibrium thermodynamic system. The question of the existence and uniqueness of its solution is still opened. The unknown in this equation is the distribution f ˇ, which corresponds to the statistical distribution of the particles over the phase space.

In the previous equation, the quantity ξ is the microscopic velocity of particles. So, the term ξ • ∇xf ˇ is the convective term, corresponding to the transport of particles in the Eulerian field. f ˇ represents the number of particles N present in a small volume of the phase space dx dξ, these are particles close to the position x and having a microscopic velocity near ξ. This leads to the following equation:

N = f ˇ(x, ξ, t) dx dξ (I.30)
A common and practical change of variable is however done. Rather than having distribution of particles, it is convenient to have distribution of mass. So, this change of variable reads:

f = mp f ˇ (I.31)
In any case, it is possible to recover the physical quantities carried by these particles.

I.3.2 Statistical moments

q θ = 2 c 2 cf (x, ξ, t) dξ (I.33f)
It may also be useful to introduce c the difference between the microscopic velocity and the macroscopic velocity (the average velocity of all particles). This difference c is defi ne d by c = (ξv).

I.3.3 Macroscopic Passage formula

In order to recover the evolution of the macroscopic quantities, one should average the Boltzmann equation. To do so, it is useful to introduce the integral operator

I [•] as follows 1 r I [A] (x, t) = ρ A(ξ)f (x, ξ, t) dξ. (I.34)
So, the evolution of a quantity carried by the distributions is obtained by the multiplication of the Boltzmann equation (eq. (I.29)) by the desired quantity, then integrating the result over the velocity space:

r r [∂ t (f ) + ξ • ∇x (f ) + g • ∇ ξ (f )] A dξ = Ω (f, f ) A dξ (I.35)
Using the Green formula, with the fact that ∇x • (Aξ) = 0 (since none depends on space) and that f = 0 at the boundary of the velocity space, one gets the following macroscopic passage formula:

r ∂ t (I [A(ξ)]) + ∇x • (I [A(ξ)ξ]) -g • I [∇ ξ (A(ξ))] = A(ξ)Ω (f, f ) A dξ (I.36)

I.3.4 Collision invariant

The collision or interaction operator can be developed within the framework of kinetic theory, id est based on the assumption that particles interact very little with each other except during an elastic collision. Said more rigorously, the free mean path is much greater than the interaction distance. This is tantamount to assume that their kinetic energy is * high in front of those of cohesion. Also, this hypothesis naturally concerns gases. An usual writing [START_REF] Boltzmann | Further Studies on the Thermal Equilibrium of Gas Molecules[END_REF][START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF] of this operator is

Ω ( f, f ) = rrr ( f 1 f 1 -f f * ) B (ϑ, w) dϑ dw dξ * (I.37)
where the subscript * denotes the particles candidate to the collision, the superscript 1 is for the post collision, B is the cross section depending on the ϑ radius and solid angle w. This last collision operator possesses a number of invariant that verify the following equation:

r Ω (f, f ) ψ dξ = 0. (I.38)

By using the conservation of mass, momentum and kinetic energy during a collision between two particles, it can be shown [Bol72; Cer88] that ψ k ∈ 1; ξ; ξ 2 . Hence, the general 2 expression of an invariant is a linear combination of 1, ξ and ξ 2 .

I.3.5 The Boltzmann H -theorem

Another important quantity worth attention concerning the collision operator was demonstrated by Boltzmann [Bol72]. This quantity is not invariant but varies monotonically and giving rise to the H-theorem of Boltzmann. Thus the information entropy of the distribution density in the Gibbs-Shannon sense is introduced using the following formula:

r H(f ) = f (x, ξ, t) ln(f ) dξ (I.39)
Then Boltzmann H-theorem is analogous to the second principle of thermodynamics and is expressed as follows: dH dt ≤ 0 (I.40)

So, the information entropy can be linked to the thermodynamic entropy of a system [START_REF] Huang | Statistical Mechanics[END_REF]:

s = -kBH + constant (I.41)
where kB is the Boltzmann constant. The entropy can only grow until reaching its maximum when the system is in thermodynamic equilibrium. By using this theorem and the linear combination of terms, one can conclude that a general expression of the invariant is

ψ k = exp ( ξ 2 \ 1C1 + C2 • ξ + C3 2 (I.42)
where C1, C2, C3 are constants.

I.3.6 Equilibrium distribution

By combining macroscopic definitions with collisional invariant and Boltzmann's Htheorem, it is possible to prove the existence and uniqueness of a thermodynamic equilibrium [START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF], and determine the constants C1, C2, C3. This state is described by means of an equilibrium distribution density that cancels the collision operator, noted f (0) . This Chapter I. Context and State of the Art 1 8

2 distribution is the famous Maxwell-Boltzmann [START_REF] Maxwell | On the Motions and Collisions of Perfectly Elastic Spheres[END_REF] distribution and can be written:

f (0) (x, ξ, t) = ρ(x, t) (2πRθ(x, t)) D/2 ( exp -(ξ -v(x, t)) 2 \ (I.43) f ˇ(0) (ξ) = N dx ( m p 2πkBθ \ D/2 ( exp - 2Rθ(x, t) mp(ξ -v) 2 2kBθ \ (I.44)
where D is the dimension of the physical space in which the problem is living, R is the universal gas constant and θ is the absolute temperature. The fact that the temperature is taken in the absolute referential avoids a change of sign within the exponential, which is essential for the distribution to remain integrable.

I.3.7 Euler equations

Such as Cercignani recalls [START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF], the examination of the Boltzmann equation projected on the collision invariant for the equilibrium distribution of Maxwell-Boltzmann allows to obtain the conservation equations of Euler of a perfect compressible gas subjected to a force field. In other words, taking the eq. (I.36) on page 16 for A = ψ k ∈ 1; ξ; ξ 2 and f = f (0) , one recovers Euler's conservation equations of a perfect compressible gas, which are

∂ t (ρ) + ∇x • (ρv) = 0 (I.45a) ∂ t (ρv) + ∇x • (ρv ⊗ v + p) = ρg (I.45b) ∂ t (E k + E θ ) + ∇x • ((E k + E θ )v + p • v) = ρg • v (I.45c)
To close the system, the equation of state of a perfect gas must be used (pV = NRθ or p = Dθ/2).

I.3.8 BGK linearization

In an eff to facilitate further numerical implementation, it is useful to simplify the collision operator through a linearization. The most commonly used and known is the one of Bhatnagar, Gross and Krook (BGK) [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF], also called simple relaxation time. It is indeed a relaxation of the current distribution towards the thermodynamic equilibrium. Any system must tend to equilibrium state (as the Boltzmann H-theorem states). Through this simplification, the collision operator is then written:

Ω (f, f ) = -ω ( f -f (0) (I.46)
where ω is the relaxation frequency, a parameter controlling the rate of convergence towards equilibrium. This parameter is also related to the frequency of collisions between particles. So the Boltzmann-BGK equation is expressed by

∂ f ∂t + ξ • ∇x (f ) + g • ∇ ξ (f ) = -ω ( f -f (0) . (I.47)

I.3.9 Chapman-Enskog expansion

To move from Euler conservation equations to Navier-Stokes-Fourier equations, it is convenient to use the Chapman-Enskog expansion [CE17] (but it is not the only possible

t f (r) + x ξ path [CE18]
). By introducing, an asymptotic development in order of small perturbations for the distribution density [START_REF] Hilbert | Begründung der kinetischen Gastheorie[END_REF], it becomes:

∞ f = f (0) + E 1 f (1) + ... = , E r f (r) (I.48) r=0
with f (1) the first order of the out-of-equilibrium distribution. E is a small parameter compared to 1 and it is related to the size of the disturbances and the Knudsen number [START_REF] Knudsen | The Kinetic Theory of Gases, Some Modern Aspects[END_REF]. A similar development must be introduced for derivation operations, hence

∞ ∂t = E 1 ∂ (1) + E 2 ∂ (2) + ... = , E k ∂ (k) ∞ ∂x = E 1 ∂ (1) + ... = , E k ∂ (k) . (I.49) t t t k=1 x x k=1
Then, the linearized Boltzmann-BGK equation (see eq. (I.47)), becomes

∞ , E k+r ∂ (k) k=1, r=0 ∞ , k=1, r=0 E k+r ξ • ∇ (k) f (r) + ∞ , k=1, r=0 E k+r g • ∇ (k) ξ f (r) = -ω ( ∞ , r=1 \ E r f (r) (I.50)
which enables a resolution in the order of E. Thus, the zero order gives 0 = 0, which express the fact that the equilibrium distribution is a solution that negates the collision operator.

At first order, the out-of-equilibrium distribution can be expressed through the following equation

∂ (1) (1)
(1) (1) . (I.51)

t f (0) + ξ • ∇x f (0) + g • ∇ f (0) = -ωf

I.3.10 Study of the first order: the NSF equations

Since the equilibrium distribution function dependency in time and space is only through the variables ρ(x, t), v(x, t) and θ(x, t); the use of the chain rule is needed. Thus, by forgetting the derivatives superscript for the sake of simplicity:

∂tf (0) = ∂ρf (0) ∂tρ + ∂v f (0) ∂tv + ∂ θ f (0) ∂tθ (I.52a) ∂xf (0) = ∂ρf (0) ∂xρ + ∂v f (0) ∂xv + ∂ θ f (0) ∂xθ (I.52b)
The temporal derivatives can be obtained thanks to the Euler equations (see eq. (I.45)) this leads, after some algebra, to

∂tρ = -∇x • (ρv) (I.53a) 1 ∂tv = -ρ ∇ x • (p) + g (I.53b) 2 ∂tθ = -v • (∇xθ) -D θ(∇xv).
(I.53c) The derivatives of the equilibrium distribution function are easily obtained and give: Once the previous systems eq. (I.53) and eq. (I.54) are injected and the equations eq. (I.52) are injected in the expression of f (1) given by eq. (I.51), it becomes:

f (0) ∂ ρ f (0) = ∂ v f (0) = ρ f (0) (ξ -v) θ f (0) (ξ -v) (I.54a) (I.54b) ∂ ξ f (0) = - ∂ θ f (0) = f (0) θ ( ( ξ -v) 2 (I.
-ωf (1) f (0) r = ρ l -∇x • (ρv) + (ξ • ∇x)ρ f (0) (ξ v) + θ r 1 l -ρ ∇x(p) + g + (ξ • ∇x)v r + f (0) ( (ξ -v) 2 D \l r 2 l f (0) (ξ -v) 2θ 2 -2θ -v • (∇xθ) -D θ(∇xv) + (ξ • ∇x)θ -g • θ r l r( 2 \ ∇xv ( c 2 \ l = f (0) c c ⊗ c -D I + θ 2θ D + 2 -2 c • ∇xθ , θ (I.55)
where c = (ξ -v) is the microscopic variation of particles velocities around the mean value. Therefore reminding that consecutive centred moments of the equilibrium distribution function are r f (0) dξ =ρ

(I.56a) R r ξf (0) dξ =ρv (I.56b) R r c ⊗ cf (0) dξ =ρθI (I.56c) R r c ⊗ c ⊗ cf (0) dξ =0 (I.56d) R r c ⊗ c ⊗ c ⊗ c f (0) dξ =ρθ 2 P 3(I ⊗ I) (I.56e) R r c ⊗ c ⊗ c ⊗ c ⊗ c f (0) dξ =0 (I.56f) R r c ⊗ c ⊗ c ⊗ c ⊗ c ⊗ c f (0) dξ =ρθ 3 P 15(I ⊗ I ⊗ I), (I.56g) R
it is possible to determine the moments of the none-equilibrium distribution function f (1) , which gives

r f (1) dξ =0 (I.57a) R r ξf (1) dξ =0 (I.57b) R c ⊗ cf (1) dξ = - ρθ R r r 1 ( 2 ∇x v + ∇x v T ) - 2 (∇ l • v)I (I.57c) c 2 f (1) dξ =0 (I.57d) R r c 2 cf (1) dξ = ρθ R 2ω r l (D + 2)∇xθ (I.57e)
where I is the identity matrix. The three null moments are consistent with the invariant of the collision operator. One can finally obtain the Navier-Stokes-Fourier (NSF) equations by using the BE at first order:

∂tρ + ∇x • (ρv) = 0 (I.58a) ∂t (ρv) + ∇x • (ρv ⊗ v + p -ν∇xv) = 0 (I.58b) ∂t (E k + E θ ) + ∇x • ((E k + E θ ) v + p • v + κ∇xθ) = 0 (I.58c)
where ν is the viscosity of the fluid that can be expressed by ν = ρθ and κ = ρθ(D+2) , when the flow is incompressible.

I.4 Lattice Boltzmann Method

The

(ξ) = √ 2πs exp ( lξl \ -2sG (I.59) (-sG) k d k w s G k k H k (ξ) = w r dξ k = w ∇ ξ w (I.60) (f k | g k ) = 2 r f :gw dξ (I.61) lH k l = H k :H k w dξ = k!sG k (I.62) (Hi,α | H j,β ) = i!sG i δij δ (i+j) (I.63)
where δ (i+j) is a generalized Kronecker symbol which is equal to 1 if α = (α1, ..., αi) is a permutation of β = (β1, ..., βj ) and 0 otherwise. The Hermite basis is a complete orthogonal basis of L 2 (R D , w). The first vectors of the Hermite basis, taken in its probabilistic form (by opposition with physical form) are given by:

H0(ξ) = 1 (I.64a) H1(ξ) = ξ (I.64b) H2(ξ) = ξξ -sGI (I.64c) H3(ξ) = ξξξ -P 3 (II) (I.64d) H4(ξ) = ξξξξ -sGP 6 (ξξI) + sG2P 3 (II) (I.64e)
where Pn (•) is an operator yielding the sum of the n cyclic permutations. The projection of the density distribution over this basis reads:

∞ f (x, ξ, t) = w(ξ) , a k (x, t):H k k=0 1 sG k k! (I.65) r a k (x, t) = f H k dξ (I.66)
where the a k are the coordinates of the density distribution on the basis. Hence, the first coordinates of the equilibrium distribution associated to the first vectors are

a0 = ρ (I.67a) a1 = ρv (I.67b) a2 = ρvv + ρ(Rθ -sG)I (I.67c) a3 = ρvvv + P 3(ρ(Rθ -sG)vI) (I.67d) a4 = ρvvvv + P 3 ( ρ(Rθ -sG) 2 II ) + P 6 (ρ(Rθ -sG)vvI) (I.67e)

I.4.1.2 Projection of the passage formula

Just like with the macroscopic passage formula eq. (I.36), it is useful to project the BGK macroscopic passage equation on the Hermite basis, which leads to

∂t (a k ) + ∇x • a k+1 + P k (∇x • (a k-1 I)) -P k (g • (a k-1 I)) = -ω ( a k -a (0) (I.68)
One could notice that from the latter equation it is possible to adopt the same process as in the previous section section I.3.9, then by going until the k = 2 in the previous equation

G r i 1 n 1
(which implies to compute coordinates a k up to the third vector), one would also obtain the NSF equations (see eq. (I.58)). Thus, it is useful to eff tuate a truncation of the density distribution up to the third order over the Hermite basis, in order to obtains the foreseen results while simplifying the distribution, and bringing the structure of a Hilbertian basis.

I.4.1.3 Truncation order

In a general context, the writing of the density distribution truncated at N -th order reads:

N f N = w(ξ) , a k (x, t):H k k=0 1 sG k k! + O ( ξ N +1 ) (I.69)
The projection over the Hermite's basis procedure enable to isolate the interesting characteristics of the BE while limiting oneself to a basis of N vectors. So, N has to be chosen according to the problem at hand. Some authors remind that the choice of N to reconstruct the macroscopic equation is not straightforward. The order N must be at least equal to 2 for recovering the weakly compressible isothermal Navier-Stokes Equations (NSE), 3 for compressible isothermal NSE and 4 for compressible NSF. For example, the fact that the third order is not enough to recover the compressible NSF is due to an error in the reconstruction of the heat flux about 20% [START_REF] Grad | On the Kinetic Theory of Rarefied Gases[END_REF].

However, these vectors H k remain, for now, continuous functions of the speed ξ. Aiming to build a numerical discrete version of the BE, it is necessary to proceed to a discretization of the microscopic speed space.

I.4.2 Gauss-Hermite quadrature

I.4.2.1 Optimization

To discretize, it is recommended to resort to a quadrature method that preserve the values of integrals. It is natural to call on the Gauss-Hermite quadrature [START_REF] Greenwood | Zeros of the Hermite Polynomials and Weights for Gauss' Mechanical Quadrature Formula[END_REF]. This quadrature corresponds to the optimization of the parameters ξ i and wi such that:

min ξ i ,wi (1 +∞ 1 1 1 -∞ w(ξ)p k (ξ) dξ q - , wip k i=1 1\ (ξ ) 1 1 1 ∀ p k ∈ R D [X] (I.70)
with q the number of points ξ i used to evaluate the function and p k any polynomials of order k in dimension D. Obviously, the larger is q the better is the approximation of the function.

The Gauss-Christoffel quadrature formulae [START_REF] Christoffel | Über Die Gaußische Quadratur Und Eine Verallgemeinerung Derselben[END_REF] give for any polynomial family (p) up to order n: pq (ξ i ) = 0 ξ i are the roots of the q-th polynomial (I.71)

r b wi = a pq (ξ i ) (ξ -ξ i ) p 1 q (ξ) w(ξ) dξ (I.72)
where a and b are the bounds of the space where polynomials exist. For Hermite polynomials a = -b = +∞.

The Gauss quadratures have an algebraic precision m such that the expression eq. (I.70) is null for any polynomial of order up to m ≤ 2q -1. So the number q must be chose n in consistence with the desired level of precision of the chosen equations to simulate. As example in 1D, if one wants to simulate the only NSE for quasi-incompressible isothermal fluids, it requires to have k = 1 in the eq. (I.68). This means that N = 2 in eq. (I.69) on the preceding page (because of the k + 1 in eq. (I.68)). Since, one will have afterwards to recover (and preserve) second order moments, it implies a minimum algebraic order for the quadrature of m = 2 + 2 = 4. Therefore, the minimum number of discrete speed will have to be q ≥ m+1 = 2.5, so q will be chosen as 3.

A reasonable generalization of the previous example, is to take q ≥ ( m+1 ) D . After determining the values of ξ i and wi, the discretized density distribution is defined by: q wi , 1

f i (x, t) = w(ξ ) f N (x, ξ i , t) = wi k=0 a k (x, t):H k G k k! . (I.73)
And following the same projection eq. (I.68) and quadrature, the forcing term can be expressed on the Hermite basis by

q gi(x, t) = wi , a k-1 (x, t):H k • g k=1 1 sG k k! . (I.74)
Thus, to retrieve the fully compressible NSF equations, it is necessary to have an algebraic precision of 8 which implies using at least a discretization with 121 microscopi c speeds in 3D or 37 in 2D. Likewise, simulations of the quasi-incompressible isothermal NSE require an algebraic precision of 4 and therefore to use at least 19 discretized speed in 3D or 9 in 2D [START_REF] Malaspinas | Lattice Boltzmann Method for the Simulation of Viscoelastic Fluid Flows[END_REF].

I.4.2.2 Usual formulations

The number of discretized speeds and the dimension chosen are commonly called the discretization schema and noted Dnqm. The Figure I.3 illustrates for a D2q5 scheme the operation of the discretization over the velocity space.

ξ 3 f 3 f 4 ξ 0 f 0 ξ 4 f 2 f 1 ξ 1 Figure I.
ξ 4 ξ 2 ξ 0 ξ 1 ξ 3 -2 -1 0 1 2 -2 -1 0 1 2 (a) D1 q3 (b) D1 q5 Figure I.4 -Lattice schemes in 1D. ξ 4 ξ 3 ξ 0 ξ 5 ξ 2 (a) D2q5 ξ 6 ξ 1 (b) D2q7 (c) D2q9 Figure I.5 -Lattice schemes in 2D.
Nevertheless, another solution to find easily the roots of the Hermite polynomials that preserves rotational isotropy is to switch from Cartesian coordinates to polar coordinate s. With such operation, one could find other schemes like D2q7 (represented on Figure I.5b), which reminds the prior numerical method used in cellular automaton (see section I.5.1 on page 30).

In order to push forward the construction of the numerical methods based on the BGK equation, it is necessary to discretized the equilibrium distribution. A very classical formulation of the discretized equilibrium function used in the scheme of second order (i.e. usual Cartesian schemes D1q3 to D3q27) is

( f i (0) = wiρ 1 + ξ i • v s G 1 (ξ i • v) 2 2sG 2 v 2 \ -2sG 1 . (I.75)
And since, for these schemes the ξ i = √ sG (±1; ±1) with ξ0 = (0; 0), by choosing the appropriate standard deviation for the Gaussian measure, one would obtain

f i (0) = wiρ ( 1 + 3ξ i • v + 9(ξ i • v) 2 2 3v 2 \ -2 . (I.76)
and the forcing term at third order reads

wiρ g = ξ G 1 + wiρ 2sG 2 wiρ r (v • ξ)ξ s G 0 r v l -s G -1 ξ 2 ξ ( ξ (v • ξ)v o • ξ \l (I.77) + 6sG 3 (v • ξ) 2 ξ + (σ:ξ ⊗ ξ)ξ -s G -1 + P 3 s G -2 -s G -1 -s G -1 ξ 0 ξ 2 ξ 1 ξ 0 ξ 3 ξ 4 ξ 2 ξ 1 ξ ξ 3 7 ξ 0 ξ 6 ξ 2 ξ 4 ξ 8 ξ ξ 5 1 y x + Chapter I. Context and State of the Art 2 6 z ξ 10 ξ 9 3 ξ 18 z ξ 3 ξ ξ 8 ξ 5 ξ 13 ξ 4 ξ 0 ξ 1 x ξ 7 ξ 2 y ξ 14 ξ 10 ξ 17 ξ 5 ξ 13 ξ 4 ξ 0 ξ 12 ξ 9 ξ 1 x ξ 16 ξ y ξ 7 ξ ξ 6 ξ 12 ξ 11 ξ 11 ξ 6 ξ 15 (a) D3q15 ξ 21 ξ 10 ξ 25 ξ 17 ξ 5 ξ 20 ξ 18 z ξ 3 ξ 4 ξ 0 ξ 22 ξ 8 ξ 9 ξ 16 ξ 2 ξ 26 (b) D3q19 ξ 19 y ξ 7 ξ 13 ξ 11 ξ 24 ξ 12 ξ 6 ξ 1 x ξ 15 ξ 14 ξ 23 (c) D3q27 Figure I.6 -Lattice schemes in 3D. ξ Scheme Roots ξ i / √ 3 Number Length Weight wi D1q3 ξ0 = (0) 1 0 w0 = ξ1-2 = (±1) 2 1 w1 2 = - 3 D1q5 ξ0 = (0) 1 0 w0 = ξ1-2 = (±1) 2 1 w1 2 = - 6 ξ1-2 = (±2) 2 1 w 3 4 = - D2q5 ξ0 = (0, 0) 1 0 w0 = ξ1-4 = (±1, 0), (0, ±1) 4 1 w 1 4 = - D2q7 ξ0 = (0, 0) 1 0 w0 = ξ1 6 = (±1, 0), (± cos ( π ) , ± sin ( π ) ) - 3 3 6 1 w1 6 = - D2q9 ξ0 = (0, 0) 1 0 w0 = ξ1-4 = (±1, 0), (0, ±1) 4 1 w1 4 = - 9 ξ1-6 = (±1, ±1) 4 √ 2 w5 8 = - D3q15 ξ0 = (0, 0, 0) 1 0 w0 = ξ1-6 = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) 6 1 w1 6 = - 9 ξ7-14 = (±1, ±1, ±1) 8 √ 3 w 7 14 = 1 - 72 D3q19 ξ0 = (0, 0, 0) 1 0 w0 = ξ1-6 = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) 6 1 w 1 6 = - ξ7-18 = (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) 12 √ 2 w 7 18 = 1 - 36 D3q27 ξ0 = (0, 0, 0) 1 0 w0 = 8 27 ξ1-6 = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) 6 1 w 1 6 = - ξ7-18 = (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) 12 √ 2 w 7 18 = 1 - 54 ξ19-26 = (±1, ±1, ±1) 8 √ 3 w 19 26 = 1 - 216
Table I.1 -Values of the weight and roots of the usual Gauss-Hermite quadratures.

After all, the velocity space discretization of the BE reads

∂f i (x, t) + ξ • ∇ (f ) + g = Ω (f , f ) (I.78) ∂t i x i i i i
Unfortunately, the second-order Hermite series expansion is not sufficient to guarantee Galilean invariance [START_REF] Dellar | Lattice Boltzmann Algorithms without Cubic Defects in Galilean Invariance on Standard Lattices[END_REF], because of an error term about O ( v 3 ) in the macroscopic equations. It also limits its use to small Mach number flows.

I.4.3 Discretization in space and time

Before or after the velocity space discretization, one can start the discretization in space and time. These two dimension (space and time) have to be discretized together because of the convective nature of the BE. Even if, the usual discretizations are accomplished on a regular grid, some works developed a LBM on an unstructured mesh [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF][START_REF] Van Der Sman | Convection-Diffusion Lattice Boltzmann Scheme for Irregular Lattices[END_REF].

A first idea to do this discretization could be a first-order finite difference. So, one could rewrite the BE discretized over the velocity space eq. (I.78) as

∂f i (x, t) + ξ • ∇ D (f ) = f (x, t) = Ω (f , f ) -g (I.79) ∂t i x i Dt i i i i i Chapter I. Context and State of the Art 2 8 Dt 2 i i 2 ω 2 f � � f � eq
where D is the particles derivative. Therefore, the first-order discretization is called the Lattice Boltzmann Equation (LBE) and reads

f i ( x + ξ i ∆t, t + ∆t) = f i (x, t) + ∆t (Ω ( f i , f i ) -gi) + O ( ∆t 2 ) . (I.80)
However, a better discretization can be easily obtained [START_REF] He | A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit[END_REF]. Indeed, using a well chosen surrogate, one can use the LBE at second order. To do so, a direct integration of the eq. (I.79) on the preceding page with the method of characteristics yields

r ∆t f i ( x + ξ i ∆t, t + ∆t) -f i (x, t) = 0 Ωi(t + s) ds (I.81)
where Ωi(t) = Ω (f i , f i )gi. Then, using a trapezoidal rule on the right-hand side, one obtains

f i ( x + ξ i ∆t, t + ∆t) -f i (x, t) = ∆t ( Ω (t + ∆t) + Ω (t) ) + O ( ∆t 3 ) (I.82)
This approximation is now at the second order but implicit. Nevertheless, thanks to the appropriate change of variable

f � i (x, t) = f i - one can find the following explicit expression ∆t Ωi(t) (I.83) 2 f � i (x + ξ i ∆t, t + ∆t) -f � i (x, t) = ∆tΩi(t) + O ( ∆t 3 ) . (I.84)
Thus, one can recognize the LBE (eq. (I.80)) for the new variable f � i . In order to dive deeper, the BGK operator needs to be exploited. For the BGK operator the previous equation can be written

f � i (x, t) = f i + ∆tω (f i (x, t) -f i (x, t)) - ∆t gi(x, t) (I.85) 2 f � i (x + ξ i ∆t, t + ∆t) -f � i (x, t) = -∆tω ( f i (x, t) -f i eq (x, t) ) . (I.86)
where f i eq (x, t) = f i eq (x, t) -g i . Therefore, to work only with the new variable, the change of variable can be inverted, which yield

2 f i (x, t) + ∆tf i eq (x, t) f i (x, t) = ω � ω + ∆t . (I.87)
In the end, injecting eq. (I.87) into eq. (I.84) one can obtain the Lattice Boltzmann-BG K Equation (LBGKE)

( i (x + ξ i ∆t, t + ∆t) = f i (x, t) -∆tω i (x, t) -f i (x, t) + O ( ∆t 3 ) ( i (x, t) -∆tω f � i (x, t) -f i eq (x, t) ( -∆t 1 - ω \ gi(x, t) + O ( ∆t 3 ) (I.88)
where 1 = 1 + ∆t . From now, the simulation will be performed over f i rather than f i .

ω ω 2
This change of variable also has an impact on the macroscopic variable. With the use of the eq. (I.85), one can compute

q , i=0 q f � i = ρ (I.89) ∆t , ξ i f � i = ρv = ρv - i=0 2 gi (I.90) q ( , ξ i ξ i f � i = p (0) -1 + i=0 ∆tω \ 2 σ (1) - ∆t 2 (g ⊗ v + v ⊗ g) (I.91) = f � eq (
These summations justify the correctional terms in the macroscopic speed and the "latti ce viscosity" that have to be compensated when fitting the physical cinematic viscosity with the collision frequency ω. This change of variable is necessary to obtain a second order accuracy. For sake of simplicity, from now f � i → f i and ω → ω, which means that the reader should read f � i rather than f i and ω rather than ω.

I.4.4 Algorithm: Streaming and colliding

From a practical point of view, it is useful to observe that the eq. (I.88) can be separated in two different parts. Indeed, one part of the equation comes from the integration along the characteristics while the other part comes from the local collision operator. Because it is simpler from an algorithmic point of view, these two parts have been virtually distinguished into a streaming (or propagation) sub-step and a collision (or interaction) sub-step. This idea is represented on the The time evolution can be simulated by a succession of collision step and streaming step described by the following equations

f i c (x, t) = f i s (x, t) -∆tω (f i s (x, t) -f i eq (x, t)) -∆t 1 -ω 2 gi(x, t) (I.92) f i s ( x + ξ i ∆t, t + ∆t) = f i c ( x, t) (I.93)
Of course, since the equilibrium distribution depends on the macroscopic quantities, between streaming and the collision steps a computation of these macroscopic variables must be performed. This separation between the streaming and colliding step also has some numerical advantages. The local nature of the collision step makes it easy to parallelize, while the streaming requires exchange of information at the edges of the domain with the other computational units.

I.5 The numerical developments of the LBM

This section is a quick review of the numerical aspects related to the LBM. More precision will be furnished when necessary in the next chapters. A detailed State of the Art has recently been published [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF], the reader wishing to have more information on the subject might refer to it.

I.5.1 The origins from the Lattice-Gas Cellular Automata

The LBM finds its origins in cellular automatons and other games of life. Indeed, the first Lattice Gas Cellular Automata (LGCA) was born in 1973 with the model HPP from the names of its creators [START_REF] Hardy | Time Evolution of a Twodimensional Model System. I. Invariant States and Time Correlation Functions[END_REF]. This latter proposes a simulation of low-density gazes, represented by balls on a regular grid. To this is added a rule determining the speed and direction of the balls after their eventual meeting at the same point on the grid at the same time. Subsequently this model gives birth to the FHP I [FHP86] model, which is followed shortly by FHP II [START_REF] Frisch | Lattice Gas Hydrodynamics in Two and Three Dimensions[END_REF]. This latter offers the simulation of hydrodynamic phenomenon rather than gaseous ones. To do this, the authors rely on a hexagonal network such as D2q7, which increases the symmetry of the model.

Although all these models are extremely fast since they enable computer programming based on Boolean variables (operations are therefore performed directly on the bits), they all suffer from a major issue. They are subject to significant digital noise due to the discretization of particles on the network.

In the first stage, this difficulty was bypassed by repeating a large number of simulati ons and performing a statistical average of the results. This naturally counterbalanced the computational performance of the method.

I.5.2 The pioneers

Lattice Boltzmann methods then appeared in an attempt to work directly on statistical distributions rather than performing multiple simulations. The pioneers in this field obtained results without having a reliable demonstration of the equations and evolution obtained [MZ88; HSB89]. They replaced the collision management rules with an equivalent operator. In 1992, the use of a collision operator based on linearization BGK [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] was introduced for the first time.

In the aftermath, it has been proved that there exists a rigorous manner to pass from the Boltzmann equation to these methods. [START_REF] He | Theory of the Lattice Boltzmann Method: From the Boltzmann Equation to the Lattice Boltzmann Equation[END_REF]. Quickly followed refinements of the understanding and the underlying mechanisms [START_REF] Shan | Discretization of the Velocity Space in the Solution of the Boltzmann Equation[END_REF][START_REF] Chen | Lattice Boltzmann Method for Fluid Flows[END_REF]. Then this opened up a new branch of the numerical simulation for fluids.

I.5.3 Boundary conditions

As in all numerical methods for resolving PDE, boundary conditions play an important role in the LBM. The determination of boundary conditions in the LBM turns out to be exigent. The mesoscopic nature of the LBM induces boundary conditions that must also be expressed at this scale. There are several types of relevant boundary conditions for fluid simulations with the LBM, including conditions to impose velocities at the edges of the domains, and others to cope with the bounce on the solid walls or periodicity conditions [START_REF] Skordos | Initial and Boundary Conditions for the Lattice Boltzmann Method[END_REF].

To imposing a speed at the edge of the computation domain, it is commonly useful to have recourse to the approach suggested by Zou and He [START_REF] Zou | On Pressure and Velocity Flow Boundary Conditions for the Lattice Boltzmann BGK Model[END_REF]. On his side, Inamuro and his team [START_REF] Inamuro | A Non-slip Boundary Condition for Lattice Boltzmann Simulations[END_REF] proposed a method to manage bounce-back condition, offering quite good accuracy but lack of stability for high Reynolds numbers [START_REF] Latt | Straight Velocity Boundaries in the Lattice Boltzmann Method[END_REF]. About bounceback conditions, it is also useful to distinguish between those performed at discretization nodes commonly referred to as full-way bounce-back [START_REF] Ginzbourg | Boundary Flow Condition Analysis for the Three-Dimensional Lattice Boltzmann Model[END_REF] and those performed midway between two nodes called halfway bounce-back [START_REF] Ladd | Numerical Simulations of Particulate Suspensions via a Discretized Boltzmann Equation. Part 1. Theoretical Foundation[END_REF]. As Succi [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF] reminded, the so-called halfway bounce-back conditions offer a second order accuracy while the full-way conditions offer only a first order accuracy. In both cases the no-slipping bounce-back condition can be used through an operator defined as

Bounce-Back(f i ) = f i = f (-ξ i ) (I.94)
where i is the opposite direction of the i-th direction in the LBM scheme.

Based on the generalization of bounce-backs at different distances between nodes, some authors have focused on curvilinear walls [START_REF] Bouzidi | Momentum Transfer of a Boltzmann-Lattice Fluid with Boundaries[END_REF]. Afterwards, boundary conditions for fixed walls allowing bounce-backs at different distances from the nodes according to the curvature of the surfaces were created and called immersed boundaries [START_REF] Feng | The Immersed Boundary-Lattice Boltzmann Method for Solving Fluid-Particles Interaction Problems[END_REF].

More recently, refinements of these boundary conditions have appeared, allowing a further increase in accuracy [MCL11; Wis+17].

I.5.4 Forces

The notion of volume force is of great importance in the LBM. Many approaches exist to integrate forces into the LBM, some using simple time integration and the use of FPD like the one brought by Shan [START_REF] Shan | Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components[END_REF] or He [START_REF] He | Discrete Boltzmann Equation Model for Nonideal Gases[END_REF] or Kupershtokh [START_REF] Kupershtokh | On Equations of State in a Lattice Boltzmann Method[END_REF], while others use the same projection based on Hermite as suggested by Guo [START_REF] Guo | Discrete Lattice Effects on the Forcing Term in the Lattice Boltzmann Method[END_REF].

The latter approach is the one developed in the section I.2, as it seems more natural in the unfolding of concepts. On the other hand, comparisons between different approaches [SS11; Son+13], show that the Guo approach is one of the most eff e, notably because it allows the incorporation of greater spatial and temporal variations.

I.5.5 Multiphase-Multicomponent

Very quickly after its inception, the LBM has been able to address complex problems such as multiphase flows. Over the years, several approaches have emerged to address this problematic. The bests known are colour-based models, Shan-Chen methods and free energy methods.

The oldest approaches are those based on colours. In 1988, Rothman and Keller [RK88] proposes a model of LGCA for immiscible fluids. Then in 1991 Gunstensen [Gun+91] suggests a similar approach for the LBM. These approaches are based on the location of each colour distribution in order to compute the gradient to perform a redistribution over time. Although effi t, this method is quite time-consuming to compute at each time step.

The Shan-Chen [SC93] approach is probably the most famous and widespread of the multiphase approaches. This uses a Green function to create a pseudo-potential between phase populations. This greatly simplifies the use of the method. Only one parameter can be used to adjust the size of the interface, which has contributed to the success of this approach. However, it has not achieved to overcome all the difficulties such as the large differences of density or viscosity between fluids and the thermodynamics spurious aspects.

The approach initiated by the Swift team [START_REF] Swift | Lattice Boltzmann Simulations of Liquid-Gas and Binary Fluid Systems[END_REF], called the free energy approach, allows to improve the differences in density or viscosity and offers in addition a physical origin of the fluid separation. This method is thermodynami cally validated. Its origin is based on the free energy of the mixture from the physical models such as Peng-Robinson [START_REF] Peng | A New Two-Constant Equation of State[END_REF] or Carnahan-Starling models [START_REF] Carnahan | Equation of State for Nonattracting Rigid Spheres[END_REF]. This energy is then derived as a local volume force that is exerted on each phase of the fluid and makes possible to separate them as well as to estimate the surface tension between them. Despite its advantages, some limitations 2 have been encountered with this approach, including the emergence of spurious speed at the interface.

On this latter point, Pooley et al. [START_REF] Pooley | Eliminating Spurious Velocities in the Free-Energy Lattice Boltzmann Method[END_REF] suggests an improvement by adding specific terms in the equilibrium distribution, that allow to reduce this spurious flows. Research in this field are active; the interested reader could refer among many to Béchereau's PhD thesis [START_REF] Bechereau | Élaboration de Méthodes Lattice Boltzmann Pour Les Écoulements Bifluides à Ratio de Densité Arbitraire[END_REF].

I.5.6 Multiple Relaxation Time

One of the main concerns with the LBM-BGK is that it is not possible, as it stands, to simulate thermal fluid flows at arbitrary Prandtl number with only one distribution. Indeed, the thermal conductivity coeffi t and the cinematic viscosity coeffi t are both driven by the same relaxation time.

The team of d'Humière and Lallemand [LL00; dHum+02], and the work of Ginzburg [Gin05; GVdH08], develop a generalization of the collision operator BGK, in order to have a better representation of their models. Soon, it has been noticed that this generalization is able to mitigate the inconvenience of the fixed Prandtl number. The BGK is then called simple relaxation operator, for the distinguo to come. This new collision operator enables, through a matrix computation, to introduce several relaxation times, and which so confers on it the name of Multiple Relaxation Time (MRT) operator. It introduces as much relaxation times as there are discrete speeds chosen over the lattice. Some of these times can have a physical meaning that allows to have variable Prandtl numbers; while others allow the improvement of stability and numerical accuracy as Dellar [START_REF] Dellar | Incompressible Limits of Lattice Boltzmann Equations Using Multiple Relaxation Times[END_REF] emphasizes.

A complete comparison between the single and multiple time relaxation operators is proposed in 2014 [START_REF] Aslan | Investigation of the Lattice Boltzmann SRT and MRT Stability for Lid Driven Cavity Flow[END_REF].

I.5.7 Entropic

Although the LBM is based on the BE and the Maxwell-Boltzmann distribution, it does not necessarily respect the Boltzmann H-theorem. To restore this principle, a variation of the collision operator is suggested [START_REF] Ansumali | Single Relaxation Time Model for Entropic Lattice Boltzmann Methods[END_REF]. This also increases method stability for high Reynolds numbers flows. This version is called entropic LBM, and is based on the idea that entropy can only increase during the physical phenomenon, which is characterized by a convergence from f to f (0) .

The Figure I.8 on the facing page illustrates the concept of the entropic LBM. The contour lines are the entropy levels of the distribution function. By definition, the maximum point of entropy corresponds to the equilibrium distribution function. The black arrow represents the distribution modification during the collision. Starting from the streamed distribution, f i s has a given entropy and after the collision the distribution f i c can not have less entropy without violated the Boltzmann's H-theorem. Thus, the solid line delimits the watershed of the authorized possible distribution, while the dashed lines are the unreachable distribution. Of course, if the ω = 1, the streamed distribution f i s will directing reach f i c = f i (0) . Many contributions have been brought to this branch of the LBM. Nevertheless, among the remarkable contributions made by the original team of this entropy LBM led by Karlin, they extended it to 3D problems [START_REF] Chikatamarla | Entropic Lattice Boltzmann Models for Hydrodynamics in Three Dimensions[END_REF], provided compatible boundary conditions [START_REF] Chikatamarla | Entropic Lattice Boltzmann Method for Turbulent Flow Simulations: Boundary Conditions[END_REF], improved Galilean invariance [START_REF] Chikatamarla | Entropy and Galilean Invariance of Lattice Boltzmann Theories[END_REF] and even addressed multi-phase prob-

I.5.8 LBM for solids

The simulation of solids with the LBM remains a task not straightforward. Several authors try different approaches. Despite some results, few works exist about this specifi field.

The oldest but probably the most furnished works are those done by Chopard's team. They replace the discrete probabilities of particles by discrete forces. Then by expressing these forces they are able to simulate complex cases like crack propagation and contact dynamics [CL99; Mar02; MC03]. However, the link with the BE is lost by this approach. Others authors work around these difficulties by coupling the LBM for fluids and the lattice spring method for the solids. Thus, they are able the obtain fluid-structure interaction like in the sphere breathing example [START_REF] Buxton | Newtonian Fluid Meets an Elastic Solid: Coupling Lattice Boltzmann and Lattice-Spring Models[END_REF][START_REF] Wu | Lattice-Boltzmann Lattice-Spring Simulations of Influence of Deformable Blockages on Blood Fluids in an Elastic Vessel[END_REF] Others approaches are experimented wave propagation in elastic solids. Some approaches are based on the knowledge of the elastic properties and the prior static solution [START_REF] Wang | Elastic Property of Multiphase Composites with Random Microstructures[END_REF]. Some other approaches are more based on the likelihood of the wave equation at macroscopic scale [START_REF] Frantziskonis | Lattice Boltzmann Method for Multimode Wave Propagation in Viscoelastic Media and in Elastic Solids[END_REF]. Some attempts also introduced the fluid-structure interaction through the pressure interface [START_REF] Velasco | A Lattice-Boltzmann Method for the Interaction between Mechanical Waves and Solid Mobile Bodies[END_REF]. And more recently, some work simulated the propagation of P-wave through the add torque force in the LBM [START_REF] Murthy | Lattice Boltzmann Method for Wave Propagation in Elastic Solids[END_REF].

I.6 Image-based diagnostic

I.6.1 Some notions about image processing

Image processing is a very broad field, as it covers lots of different scientific areas. Among al l the topics composing this field, denoising, classifi segmentation and morphological image processing constitute the main areas that are used in the context of medical image processing, and that can be addressed with the LBM. In this section is presented the idea lying behind these three topics and some techniques that are representative of the concepts. Literally, denoising signifies the action of removing noise. In other terms, it is linked to image regularization. The goal of this process is to improve regularity without losing information contained in the image. To this aim, different methods can be employed: filtering, diffusion,... A comprehensive review of denoising algorithms is given in [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]. One can cite one of these methods, anisotropic diffusion, as it is an application that can be directly addressed with the LBM. This point will be detailed in the next section I.6.1.5. Anisotropic diffusion [START_REF] Perona | Anisotropic Diffusion[END_REF] is basically a diffusion algorithm using edges of the image to control the diffusion velocities. More precisely, the algorithm increases diffusion velocity when the gradient is low and vice versa.

I.6.1.2 Classification

Classifi is used in different types of problems of image processing. More spe ci fi , the goal of classifi is to make a decision about the belonging to a part of an image to a specifi class. Classifi algorithms are more willingly associated with machine learning and Pattern recognition as they are often based on a priori knowledge in the form of databases. However, in several cases classifi can be done in an "unsupervised" way, without any a priori knowledge. A review focusing on medical images classi fi cation methods can be found in [START_REF] Smitha | A Review of Medical Image Classifi cation Techniques[END_REF].

I.6.1.3 Segmentation

The segmentation principle aims at regrouping parts of the images that have common characteristics and disconnect the other parts. The idea is to "make sense" from an image by selecting parts of interest. A review on basic image segmentation techniques can be found in [START_REF] Pal | A Review on Image Segmentation Techniques[END_REF]. Among all the different image segmentation techniques, one can cite a simple one that is widely used, the Otsu method. The Otsu method [Ots79] is based on the separation of an image into several classes using its histogram. The algorithm searches for the minimum intra-class variance, which is defined as a weighted sum of the variance s of each class. Thus, this method is global and does not take into account the spatial aspect of the image.

I.6.1.4 Morphological image processing

Morphological mathematics [START_REF] Serra | Mathematical Morphology and Its Applications to Image Processing[END_REF] uses basic principles of the set theory to describe morphological operators. It is mainly based on two operators: dilation, which is equivalent to a Minkowski addition, and erosion, which is the dual operation. From these two operators, it is possible to defi a lot of more complex ones such as the opening, closing, skeletonization... In fact, morphological image processing can be viewed as two different things. The first one is a mathematical theory which is really comprehensive and can be useful to define new algorithms but also establish proofs of existing ones. The second is a set of basic tools that allows to manipulate spatial information easily in the image.

I.6.1.5 Image processing with the lattice Boltzmann method

Due to the ability of the LBM to solve PDE of a diffusive-convective nature such as the BE, some authors have had the idea of reducing this characteristic to only diffusive nature that is useful in image processing. Thus, the pioneers of the field [JLS99] used anisotropic diffusion for image denoising. Afterwards, Zhao [START_REF] Zhao | Lattice Boltzmann Based PDE Solver on the GPU[END_REF] then Chen [START_REF] Chen | An Anisotropic Diffusion Model for Medical Image Smoothing by Using the Lattice Boltzmann Method[END_REF] proposed an improvement by introducing a bouncing propagation probability, which allows it to do image editing in particular and image segmentation.

The introduction of concepts such as total variation [START_REF] Chang | A Lattice Boltzmann Method for Image Denoising[END_REF] are used for contour detection and filtering [START_REF] Chen | Lattice Boltzmann Method for Filtering and Contour Detection of the Natural Images[END_REF].

Some authors artfully used the BE with special definitions of forces to implement a dithering of images [START_REF] Hagenburg | A Lattice Boltzmann Model for Rotationally Invariant Dithering[END_REF] The optical flow is also introduced, based on the Vlasov equation, i.e. a collisionless Boltzmann equation [START_REF] Ding | Optical Flow Analysis Based on Lattice Boltzmann Method and Lower Order Approximation with Relaxation Factors[END_REF]. An example of the use of the latter concept has also been provided in the field of crowd motion analysis [START_REF] Xue | Abnormal Prediction of Dense Crowd Videos by a Purpose-Driven Lattice Boltzmann Model[END_REF]. Very recently, the BE has been used to improve the stability of motion analysis on the hand [START_REF] Alahyane | A Lattice Boltzmann Method Applied to the Fluid Image Registration[END_REF].

I.6.2 Image-based diagnostic with macroscopic methods

Image-based diagnostic is a task requiring many specialists, and it is more often separate d in different domains. Some recent work uses many sensing systems such as tracking se nsors and pressure plate to calibrate and validate their FEM model of feet [START_REF] Akrami | Subject-Specifi Finite Element Modelling of the Human Foot Complex during Walking: Sensitivity Analysis of Material Properties, Boundary and Loading Conditions[END_REF].

Some authors started from X-ray and CT-scan images to diagnostic the resistance of bones to solicitations [Ulr+98; Pis+04]. Others used CT-scans to create a patient specifi model and then perform FEM to predict the deployment of stent-graft inside patient aneurysms [START_REF] Perrin | Patient-Specifi Numerical Simulation of Stent-Graft Deployment: Validation on Three Clinical Cases[END_REF].

Knowing the eff and time saved by meshless methods, some works develop methods close to the FEM but meshless. An example called the Total Lagrangian initially developed to help surgical preparation from medical images [START_REF] Horton | A Meshless Total Lagrangian Explicit Dynamics Algorithm for Surgical Simulation[END_REF][START_REF] Miller | Beyond Finite Elements: A Comprehensive, Patient-Specifi Neurosurgical Simulation Utilizing a Meshless Method[END_REF] offers acceptable numerical errors. The same team also used images of the tibia to create a discretized geometry and mesh but to perform computations on this bone with FEM [START_REF] Wittek | From Finite Element Meshes to Clouds of Points: A Review of Methods for Generation of Computational Biomechanics Models for Patient-Specifi Applications[END_REF].

I.7 Emanating questions

Image processing tools have been experimented with the LBM. However, image-based diagnostic requires more exploration in this direction. The Mathematical Morphology (MM) is a powerful theory useful in both image processing and numerical simulation. So, one can wonder if it is possible to create a bridge between the MM and LBM. Also, what are the image processing tools that the LBM can incorporate to simplify the images? Assuming that, the answer to the previous questions are reachable, is it possible to create a complete pipeline that starts from the images to numerical simulations using only the LBM technic.

Diagnostics may concern all kinds of matter. If the LBM is famous for its effi with fluids, a lot of eff still have to be done about solids. Also, is it possible to reach, with the LBM, any stress tensor as a solid equilibrium state in first approach? But afterwards, is it possible to simulate the solid dynamics from the LBM?

These opened questions can be gathered around a single subject: the lattice Boltzmann method for numerical simulation of continuum medium aiming image-based diagnostics. Indeed, this subject widely embedded the previous questions and the eventual links between the fluids and the images. 6

I.8 Conclusion of the chapter

This chapter drew the main principles and equations of the macroscopic mechanics. The main numerical methods to deal with the macroscopic mechanics were reviewed. Then, the capacities of the Boltzmann equation to solve indirectly the macroscopic NSF equations was underlined, such as the path to the LBM. From the current developments around the LBM, it appears that even if the LBM is proved to solve numerically the NSF equations from a mesoscopic scale; it has some challenges to overcome before to be able to simulate solid as do the macroscopic methods. It also appears that the LBM is an efficient method for image processing. But some extensions are required to build a complete tool for image-based diagnostic with only the LBM. Before to be able to perform the developments necessary to investigate these wondering, it is compulsory to possess a reliable computational basis. This is the purpose of the next chapter.

Abstract of the chapter

Further developments are required to answer the questions raised by the subject of a pipeline between the images processing and numerical simulations, so a reliable LBM tool s is required. Since, no existing LBM software solution is able to deal with image processi ng, the creation of this pipeline implies the implementation of such kind of solution. This chapter describes roughly the implementation done.

The implementations done are divided in three groups. The first one is about the LBM kernel itself while the second is about the environment wrapped around it and the last deals with the graphical interface. Then, a benchmark of usual test cases to certify the implementation done, is browsed. This benchmark is made of a hydrostatic pressure profile, a Poiseuille flow under gravity and Zou-He boundary conditions and Von Kármán vortex street. The results of the tests allow to validate the implemented tool and go furthe r with the developments.

Résumé du chapitre

Des développements supplémentaires sont nécessaires pour répondre aux questions soulevées par le sujet d'un pipeline entre le traitement des images et les simulations numériques, d'où la nécessité d'outils fiables en LBM. Étant donné qu'aucun logiciel LBM existant n'est capable de gérer le traitement des images, la création de ce pipeline implique l'implémentation d'une telle solution. Ce chapitre décrit grossièrement la mise en oeuvre eff Les implémentations eff sont divisées en trois groupes. Le premier concerne le noyau LBM lui-même tandis que le second concerne l'environnement qui l'entoure et le dernier concerne l'interface graphique. Ensuite, une base de comparaison consititué de cas tests habituels est parcouru pour valider l'implémentation eff Cette base de comparaison est constituée d'un profil de pression hydrostatique, d'un écoulement de Poiseuille entrainé par gravité ou par des conditions aux limites de Zou-He et de l'allée de tourbillons de Von Kármán. Les résultats des tests permettent de valider l'outil implémenté et d'aller plus loin dans les développements.

II.1 Introduction

The previous chapter states on the idea that the LBM for image-based diagnostic is something that have to be created. Aiming this construction, some developments will have to be done. So, before starting the modifications of a LBM, it is necessary to validate this initial code to lay on a solid ground. This is the purpose of this chapter.

Before to list the main works done about the implementation, the specifi imposed by the global aims are explained. Then, some technical choices and aspects are discussed leading to a global structure and capacities of the proposed tool. And to ensure the reliability of the tool, a validation around a benchmark is carried.

II.2 Specifications and paradigms

In order to create new tools to respond to the problematic, it is necessary to question about the implementation of a version of the LBM algorithm. Indeed, none of the existing LBM solutions is able, for now, to perform both image processing and CFD. The main advantage to implement a new version of the LBM algorithm is to be able to master in detail the whole underlying mechanism. But the main drawback is the necessity to implement from scratch all the basics of the existing method. The use of established codes like Palabos or OpenLBM would make save considerable development time to the detriment of understanding and pedagogy.

At the end, the choice of a new implementation has been made. So, to frame the implementations, the following requirements are chosen:

• cross-platform implementation that allows to work on Windows and Unix (Linux-like and MacOS) indifferently,

• fast and effi t computational tool,

• generic tool that can be adapted to work with 3D images,

• free and open if possible.

• be able to work with or without a graphical user interface (computational clusters usually do not have graphical devices).

Some of these requirements drive technical choices, while others ask some supplementary work and precautions. The proposed solution will then have to be tested in several benchmarks test cases in order to prove its effi . Before diving deeper in the technical aspects of the implementation, it is useful to remind the nature of the LBM core. A usual LBM pseudo-code is presented in Algorithm II.1. The LBM is traditionally exploited through an explicit time iteration loop. This loop is preceded by some initialization. One iteration of the latter loop is divided in 5 operations that are defined in the previous chapter: the computation of the equilibrium, the collision, the streaming, the application of the boundary conditions and the computation of the macroscopic variables.

II.3 Structure and capacities

From a first look, the implementations can be divided in 3 groups. The first group is the computational core itself that runs the LBM simulations. The second is the environment that manages all the files and operations. The last is the graphical user interface that he l ps the use of the two previous groups.

II.3.1 Computational code

II.3.1.1 Choices: programming language and compilers

The LBM, like most of the numerical simulation tools, is dealing with matrix operations to solve the discretized equations; so, the most suitable programming language is the FORmula TRANslator (Fortran). This choice is also in the very heart of the other requirements for the desired tool: cross-platform, fast, free and open. Indeed, the Fortran language is the fastest one when it comes to matrix operations, but another choice could be well-known C-type languages.

For every compiled programming language like Fortran or C, several compilers exist, each with its specifi Three possible compilers have been studied:

• the Intel compiler which provides faster optimization but is not free nor open;

• the PGI, which is free and not open, provides effi t accelerators but still contained bugs in 2018 (specially about pointers);

• the GNU compiler is free, open compatible with the accelerators but a bit less effi cient.

The Fortran implementation of the computational code has been fulfi in order to be compatible with these three compilers. Thus, the code is sure to be accepted by the standards and the needs of everyone.

II.3.1.2 Libraries

The implementation of this computational core required the writing of a lot of libraries usually missing in Fortran. This is mainly due to the decrease in the Fortran use over the last years. One illustrative example is the libraries to read and write a standard file that stores field data in 3D. A very common choice is to use of Visualization ToolKit (VTK) format. Unfortunately, there is no standard libraries for this format. This leads to a complete re-implementation of the following libraries:

• PorPre: The implementation of a library called Portable Precision (PorPre) insures the portability of the code through all the platform and compilers, as recommended by many authors.

• FoST: A library implemented is the Fortran String Transformation (FoST) ease the manipulation of String in Fortran which is certainly one of its biggest weak points. Furthermore, this library also allows to create arrays of a chain of characters with different lengths.

• Logger: To be able to communicate with the users even when multi-processing is enabled, it is necessary to have a thread-safe protocol. This leads to the creation of the library Logger, which allows to sort different levels of information and printing messages to prompt and log files at the same time even in pure thread-safe procedures through a saved buffer.

• VTK_IO: To be able to visualize the results in 3D in an effi t way, the VTK solution is classical in then numerical simulation field. Indeed, this standard gives the possibility to read and write tensor data fields in parallel with all mesh types. Plus, this standard provides binary and American Standard Code for Information Interchange (ASCII) format of data in the so-called "legacy" and eXtensible Markup Language (XML) types. Thus, a re-implementation providing all these features wi th a complete type structure was needed.

• FoXML: The VTK and the input for commands needed a complete implementation of an XML for input and output. This allows to have a well-structured input file.

Furthermore, the reading of images was necessary, which leads to the creation of a library called Fortran Image Library (FIL). This library was partially based on a wellimplemented existing library for Graphics Interchange Format (GIF) files.

II.3.1.3 Optimization

Since one of the objectives is the performance of the computational core, some fundamental optimization was made. The first simple optimization realized was to avoid doing loops, which is time consuming. When it comes to generic LBM algorithms, loops are nested in many different actions at least because it is necessary to range the discretized velocities (one example is during the computation of macroscopic variables from the Dnqm scheme). The first idea would be to use an if-conditional test and put a sequential procedure to compute macroscopic variables. Nevertheless, this idea implies a loss of performance due to the evaluation of the conditional test, if they are in the time loop. Thus, another solution that maintains flexibility and performance has been adopted: a procedure pointer is correctly confi earlier, when the Dnqm scheme is given. This railway switch is operated by the module called Option_Selector. The same kind of method is used in I/O and Logger libraries to save some computational time.

The only part that cannot be done completely sequentially is the application of the boundary conditions and followed by the appropriate macroscopic variables. Indeed, the se steps depend on each simulation with its boundary conditions and equations to solve. To do such procedures with performance and flexibility, the choice of creating a list of procedure s associated to a cloud of points on which the procedure will be applied was done. Thi s l i st of procedures can then be executed in the correct order.

All these libraries and modules are organized around the structure represented on the Figure II.1.

In order to improve further the performance, the Fortran programming language offers the possibility to create procedures called "pure". This kind of declaration for procedures insures that it has no side eff and it is thread-safe. These declarations give to the compiler the opportunity to proceed a stronger and a deeper optimization. All the procedures in the LBM time loop were declared and implemented as pure. This kind of procedure offers another advantage, since they are thread-safe, they are naturally ready and compatible with parallel computing. The choice of a parallel programming of heterogeneous CPU/GPU was made. This is provided by the Open ACCelerators (OpenACC) library

II.3.2 Environment

The computational tool is wrapped in an environment in order to cope with the compilers, the management of files, plotting of results or interact with other numerical tools. This environment is written entirely in Python, which provides cross-platform compatibility and a good interoperability with other languages or numerical tools. Moreover, it is intrinsically cross-platform.

The Python environment is made of multiple bricks:

• compilation tool: the compilation brick for the Fortran kernel which accept any compiler and any option of compilation and can be adapted easily to any code that needs to be compiled.

• execution tool: Another brick is made for the execution of the computational tool, which is operated in a temporary folder with all the extra files needed in order to protect the architecture in case of a crash.

• batch: A batch script allows the sequentiality of the actions (compilation, execution, plotting...). The main advantage of the batch is based on its confi file which makes it possible to work without a graphical user interface: clusters or servers do not have graphical access commonly. • database: One major risk in the implementation of a computational tool is the regression. This means that between two versions of the codes some functionality or precision is lost. To handle this issue, a database containing the results of the previous simulation of the benchmarks is stored and a methodical comparison is made at every new version of the code. Thus, the evolution of the whole software is tracked and the insurance of no regression is provided.

• installer: The last environment part lead off is an installer that has been prototyped to work with Python. This installer has been tested on Linux, more work is nece ssary to finalize it and create a Windows version.

As can be seen, eff have been made to ensure that this environment provides a solid foundation for future developments, code enhancements and distribution.

II.3.3 Graphical User Interface

Since it is much more convenient to work with a Graphical User Interface (GUI), it was more comfortable to lend one to the software. In addition, as the final goal of the software is to work with images, the most direct way of manipulating images is using a GUI. Such a GUI is implemented using the Qt graphical technology and the Python core which is embedded in the PyQt5 library. Thanks to this library, a powerful, ergonomic, following standards and cross-platform interface could be developed. Through some dedicated software that Qt incorporates, it is a flexible and easily maintainable tool.

The GUI can edit the XML input files for the computational tool providing a link with the computational core. But the communication with the entire environment is made through the batch file. This implies that the GUI is generalizable to many numerical simulator if an associate confi file for the batch is furnished. Moreover, to not interrupt the interactions between the user and the interface duri ng the execution of the batch, a multi-thread implementation is made. Another feature of the GUI is its link with the database. Indeed, it is possible to launch all the benchmark test cases from the GUI, and plot information from the database associated to the benchmark, to graphically track the evolution of the tool. The

II.4 Benchmark to validate the initial code

To verify the correct behaviour of the implemented code, this one is tested on multiple test cases. These test cases allow to validate different known situations and components of the code. This validation has two objectives, certify an appropriate basis for further development and insure the non-regressivity of the code to avoid the propagation of errors.

Environment Group

GUI database

The number of test cases, the level of detail or numerical precision are reduced, since this is not the originality nor the purpose of the thesis work. So, the benchmark developed here looks at the hydrostatic pressure, the Poiseuille flow under gravity, the Poiseuille flow with Zou-He boundary condition and the Von Kármán vortex street. In this first approach, only the BGK operator is used. Moreover, no 3D test is included here, since no 3D simulations are relevant yet, in the next chapters.

II.4.1 Hydrostatic pressure

To evaluate the capacity of the implemented LBM algorithm to compute a static solution in hydrodynamics, it is usual to simulate the case of the hydrostatic pressure of a column of water in a gravity field. Under the assumption of static solution of a weakly compressible fluid at rest (the speed flow is null), the Navier-Stokes equation becomes chosen lead to a variation of pressure along the column of water: ∆p = 3.820.10 -2 Pa.

According to the parameters, the simulation length is 118 cells plus two for the ful l way bounce-back condition. This yields the non-dimensional result

∆p * = ∆p.∆x -2 .∆t 2 = 3.508.10 -3 (II.3)
Starting from a resting initial state (without forces), after 15 000 iterations and using the Guo formulation of forces, the algorithm returns a variation of pressure of 3.4185.10 -3 . This means a relative error of 2.54%.

The Figure II.5 represents the field of pressure obtained along the column of water. This field is varying linearly with the altitude x, and the iso-pressure lines are horizontal, according to the theory.

II.4.2 Poiseuille flow under gravity

A short path between the hydrostatics and hydrodynamics is to reuse the gravitational forces, in order to create a stationary Poiseuille flow. A 2D Newtonian viscous fluid weakly compressible in a cylinder of diameter H with periodic boundaries reaches a stationary state and the Navier-Stokes equations are reduced to ρ∇x • (2ν∇x (v)) = ρg.

(II.4)

Using the boundary conditions stipulating that the speed of the viscous fluid is null at the interface with the pipe, the stationary flow is described by the equations:

g v(y) = 2ν y (y -H) (II.5) g ( H \ 2 R e = 2ν

II.4.3 Poiseuille flow with Zou-He conditions

Imposing speed at the border of the simulation is a necessary boundary condition to test. Probably the most famous is the one suggested by Zou and He [START_REF] Zou | On Pressure and Velocity Flow Boundary Conditions for the Lattice Boltzmann BGK Model[END_REF]. This boundary is LBM Theory

•10 -2 x v * x (y) v *
expressed in D2q9 at the top edge by the following equations

1 ρ = f 0 + f 2 + f 4 + 2 (f 6 + f 3 + f 7) 1 -v (II.8a) 2 f 1 = f 3 + 3 ρv (II.8b) 1 1 f 5 = f 7 + 2 (f 4 -f 2) + 6 ρv (II.8c) 1 1 f 8 = f 6 -2 (f 4 -f 2) + 6 ρv (II.8d)
In order to verify the implementation of this condition, it is possible to simulate again Poiseuille flow. For this simulation, the parabolic speed profile is imposed at the inlet of the pipe, starting from the rest initial condition (v = 0), waiting for the steady flow to establish itself, and the speed profile is measured at the outlet. This outlet speed profile should be a parabolic Poiseuille profile equal to the one given at the inlet due to the steadiness of the flow. The maximum speed resulting in the outlet is 1.4757.10 -1 , this means a relative error of 1.62%. This relative error is larger than the previous test case. A likely explanation is the fact that the periodic boundary condition was not used in the simulation.

II.4.4 Von Kármán vortex street

Since the temporal dynamic analytical solutions in fluid dynamics are not easy to obtain, one classical test case is the Von Kármán vortex street. This allows to validate the dynamic behaviour of the LBM algorithm. Also, in this situation, the dimensionless Reynolds number is given by Re = v max L (II.9) ν where L is the characteristic length of the cylindrical obstacle in the flow path. And the experiments establish different kind of flow in function of the Reynolds number:

• if Re < 1 the flow is in a laminar regime, the streamlines flow around the obstacle without any perturbation;

• if 1 < Re < 80 the flow is in a transitive regime, two symmetrical recirculation rotating in opposite sense appears behind the obstacle;

• if 80 < Re < 10 000 the flow is in a oscillating regime, two vortices are successively detaching from the recirculation behind the obstacle.

From a practical point of view, the Zou-He condition is applied to the left and right borders of the physical domain, while the top and bottom border use the full way bounceback condition. It has to be noticed that a zero speed gradient condition at the outlet would be a better condition, but not considered here for simplicity. The full-way bounceback condition is also used on the cells constituting the hard cylinder in the pipe. The pipe is constructed from 250 cells in the x-direction and 75 cells in the y-direction. The cylinder diameter is 15 cells. Table II.3 -Parameters taken for the Von Kármán vortex street different regimes.

II.5 Conclusion of the chapter

This chapter emphasis the main choices and feature of the LBM implementation done. The Fortran implementations are compatible with the 3 principal compilers and with the library OpenACC. Some optimization allows to keep 70% of the simplest sequential codes computational speed. The LBM computational tool is wrapped in a Python environment which allows to manage notably a database which verifi the non-regressivity of the computational tool on all platforms. In addition, the whole can be controlled by a batch script or a GUI. All these tools are rather generic and free. Thus, the initial specifi are well respected.

The implemented solution is tested with a benchmark of different situations such as hydrostatic pressure profile, a Poiseuille flow under gravity and Zou-He boundary conditions and Von Kármán vortex street. Through these series of benchmark tests, the implemented solution presents acceptable errors to be the foundation of further modifications and investigations related to image-based diagnostic.

III.1 Introduction

The previous chapter validated the initial implementation of the LBM tool. So, the coming improvements are based a solid ground. The discussions in this chapter are mainly around the interactions and interconnections between the LBM and image processing. Indeed, the will to build a pipeline from images to image-based diagnostics with the LBM requires to introduce unprecedented concepts. The construction of the pipeline is done here, but its construction is divided in successive steps. These steps are the simulation of CFD directly from images, and the processing of images with the LBM. This last step is enhanced here by the MM tool which has never been translated to the LBM. This enhancement opens the gate to further development detailed in the aftermath.

III.2 The image support as an LBM network

III.2.1 Density as a mathematical image of the "image"

To apprehend the link between images and the LBM it is necessary to come back on certain points. A grey-level image I in dimension D is a function that goes from R D into R (usually in [0, 1]). It is natural, in two dimensions, for discrete images to be considered as a square structured network, in which each pixel is surrounded by its eight neighbours. Furthermore, the LBM is constructed on a network commonly squared constituting a lattice (hence its name). In cases such as those encountered in isothermal and weakly compressible CFD, this lattice (let us say in D2q9) is also built on the connection of one point with its eight closer neighbours. Although, the interpretation of images in terms of density of particles f , is not simple at first glance; at a fixed time t, the density ρ(x, t) becomes on a spatial function going from R D into R and sharing the structure of the network (if discretized). Thus, the notion of density and image coincide in that sense. Then, the density of particles f can be interpreted as the information exchanged between two pixels. These notions can be easily generalized, without losing the link between LBM and images, in this way an image is the density of particles at a given time.

III.2.2 LBM simulation from images

Now that the link between the images and the LBM is exhibited, some simple but original simulation can be performed. The previous section highlights the link between the density and the image. Therefore, using the image as an initial state of the matter, and using the structured network provided by the image lattice, it is easy to construct LBM simulations. An illustrative example from medical domain is shown on the Figure III.1. This example is given without real biological meaning since the input parameters chosen were too coarse. It shows an image of an aortic aneurysm from a Computed Tomography-Scan, where the area of interest is segmented and its border is associated to bounce-back conditions. Then, inlet and outlet impose a constant velocity flow.

Moreover, the usage of the Boltzmann equation allows the resolution of some PDE used in image processing. Several differential equations are useful for denoising, inpainting, dithering, segmentation or for the Geodesic Active Contours (GAC).

III.3 The usage of the LBM for image processing

III.3.1 Anisotropic diffusion equation

LBM anisotropic diffusion for image processing has been implemented for image pre+processing such as smoothing and filtering operations [CYQ08; NCR15] and segmentation [START_REF] Wang | Segmentation of Giant Cerebral Aneurysms Using a Multilevel Object Detection Scheme Based on Lattice Boltzmann Method[END_REF][START_REF] Wang | Multilevel Segmentation of Intracranial Aneurysms in CT Angiography Images: The Segmentation of Intracranial Aneurysms[END_REF].

In segmentation and denoising of images, the anisotropic diffusion-reaction equation can be used. In such case, it is necessary to resort to the following discretized version of the BE

f i (x + ξ i ∆t, t + ∆t) = Gi (x) f i (x, t) + ω f i (0) (x, t) -f i (x, t) + [1 -Gi (x + ξ i ∆t)] f i (x + ξ i ∆t, t) + ω f i (0) (x + ξ i ∆t, t) -f i (x + ξ i ∆t, t) + a∆t [T -ρ (x, t)] (III.1)
where a is a parameter setting the threshold reaction, Gi is the probability of propagation (or the probability of bounce-back) and T is the desired level of threshold for the segmentation. The f i are still the distribution over the lattice and ξ the discretized velocities. By choosing the propagation probability parameter, one can set the ratio of "particles" transported by the f i to the next cell. As Chen [CYQ08] explains, the Gi can be illustrated by Figure III.2 but remains an artefact. This has the eff of improving the numerical stability of the model.

Moreover, the computation of the probability of propagation is accomplished via the following expression:

Gi (x) = r 1 + | G ρ (x, t) -G ρ (x + ξ i ∆t, t)| l -2 10 (III.2)
where G is the discrete matrix gradient operator. So, the convolution of G with ρ corresponds to the evaluation of the gradient of ρ. Therefore, the probability of propagation is related to the variation of the gradient of ρ.

To explain the genesis of this method, the Chapman-Enskog and Taylor expansions (or the projection on the Hermite basis) are applied to the LBGKE eq. (I.88) on page 28 and the equation of diffusion is recovered:

∂ρ ∆x 2 ( 1 = \ -0.5 ∇ ρ.
(III.3) ∂t 3∆t ω Then with a variable relaxation frequency ω, the anisotropic diffusion equation can be reconstructed.

III.3.2 Example of anisotropic diffusion

In order to illustrate the capacity of anisotropic diffusion to segmented images, by choosing the set of parameters presented in Table III

-1 -1 -1 G -1 8 -1 -1 -1 -1 a 0.5
Table III.1 -Parameters taken for the anisotropic diffusion reaction for segmentation.

III.3.3 Amelioration through volume rendering and level set

More recently a LBM segmentation tool has been elaborated to capture a GAC model. This new method is applied with success to the segmentation of aneurysms, and relevant scientific results are published [START_REF] Chen | Segmentation of the Thrombus of Giant Intracranial Aneurysms from CT Angiography Scans with Lattice Boltzmann Method[END_REF]. Another effi t approach for image segmentation is a combination of the anisotropic diffusion with threshold driving a volume rendering algorithm coupled with an Intensity Gradient Method (IGM). Such an approach has been recently explored by Ge [Ge+17] , and offers promising results.

III.4 Mathematical morphology and the LBM

The MM introduced by Matheron in 1967 [START_REF] Matheron | Elements pour une theorie des milieux poreux[END_REF] is a mathematical theory essential in image processing, still to this day. This theory is based on two elementary operations called erosion and dilation. From these two basic operators, many other more complex operators can be built, allowing their turn applications such as segmentation, denoising, skeletonization, gradient calculation, etc. [START_REF] Najman | Mathematical Morphology: From Theory to Applications[END_REF].

III.4.1 Reminding on the mathematical morphology

A grey-level image I in D dimensions, is mathematically defined as a function of R D in R. Furthermore, a greyscale structuring element S f is a function that goes from R D into R and such as the set

{ x ∈ R D | S f (x) /= {∞} ∪ {-∞} 1 is bounded.
Then the morphological dilation of a grey-level image I by a functional structuring element S f , is denoted (I⊕ S f )

y i i i f i i i z f
and is defined by the relation

(I⊕ S f ) (x) = sup (I (y) + S f (x -y)) .
(III.4) y Otherwise, a morphological erosion of an grey-level image I by a functional structuring element S f , is denoted (I8 S f ) and is defined such as

(I8 S f ) (x) = inf (I (y) -S f (y -x)) (III.5)

III.4.2 Equivalence between MM and LBM -LB3M

Proposition 1. One increment of time in the LBM at the equilibrium using the following definitions of the first moment and the equilibrium function:

ρ (x, t) = max {f i (x, t)} (III.6) f i (0) (x, t) = ρ (x, t) + S f (ξ i ) ; (III.7)
is equivalent to a morphological grey-level dilation of the image ρ by the structuring element S f .

Proof. In order to construct the equivalency of the greyscale morphological dilation operator with the LBM, it is necessary to adapt some operations in the different step of the LBM. Thus, by considering the density ρ no more as the sum over the probabilities f i , but as the maximum over these probabilities, this leads to the writing of the equation eq. (III.6). It is also necessary to integrate the evaluation of the image through the structuring element, which is introduced via the equation eq. (III.7). These simple modifications, without any changes in the collision and streaming steps, allow to obtain a new density after one time step; this latter density can be expressed in terms of the different steps constituting the LBM at equilibrium (i.e. for ω = 1). Thus, by rolling back successively the equations composing the LBM, one can obtain

ρ (x, t + ∆t) = max {f i s (x, t + ∆t)} = max {f i c (x -ξ i ∆t, t)} = max ( -ω f i (0) -i l (x -ξ ∆t, t) (III.8) i = max f (0) (x -ξ ∆t, t) i i i = max {ρ (x -ξ i ∆t, t) + S f (ξ i )}
Moreover, by reusing the definition of the morphological dilation and by applying the change of variable z = xy, it becomes possible to write this definition in its discrete form

(I⊕ S f ) (x) = max (I (x -z) + S f (z)) (III.9)
This last expression shows the equivalence between the suggested method and the greyscale morphological dilation.

In the same manner, it is possible to define the greyscale morphological erosion through the LBM. i i Proposition 2. One increment of time with the LBM at the equilibrium using the following definitions of the moments and the equilibrium function:

ρ (x, t) = min {f i (x, t)} (III.10) f i (0) (x, t) =ρ (x, t) -S f (-ξ i ) (III.11)
is equivalent to a morphological grey-level erosion of the image ρ by the structuring element S f .

Proof. By proceeding like in the previous demonstration, one can roll back the consti tuti ng LBM equations. Thus, one obtains under the same hypothesis than previously (ω = 1), the expression of the density after a time step

ρ (x, t + ∆t) = min {ρ (x -ξ i ∆t, t) -S f (-ξ i )} . (III.12)
This last expression also coincides with the discrete expression of the greyscale morphological erosion after the change of variable z = xy.

Therefore, it is possible to define rigorously the two operations funding the mathematical morphology through the LBM formalism. One can also notice the element S f (-ξ i ) is nothing more than bounce-back operator applied to S f (ξ i ).

Whole these set provides the elementary operation constituting the link between the LBM and the MM, which leads to the Lattice Boltzmann Method and Mathematical Morphology (LB3M) . Those equations are summed up in the following algorithm Algorithm III.1.

III. 1. 1 function LBM_DilationErosion()():

Data : Lattice and initial conditions Result : Density field Since, the MM lays down on these two operations, this is one way to prove that the LBM can provide MM operations. But other methods are imaginable such as integrating the structuring element in an external force. Indeed, this affects the collision operation (specially with Guo [GZS02] formulation of forces). So, the same result can be obtained. In the same manner, the use of the continuous version in time and space for the Boltzmann equation is thinkable.

f i ← Collide(f i , f i ( 0
Another remark can be mentioned. In comparison with the anisotropic diffusion, which is known to conserve the global grey-level of the image; the MM does not. Indeed, the anisotropic diffusion equation resulting from the BE conserves the density (or grey-level) through the Chapman-Enskog expansion (see eqs. (I.57a) and (I.58a) when v = 0). But, this property is not true with the MM, since the purpose of the dilation or erosion is to gain or lose some grey-levels through the structuring element. Thus, the MM does not conserve the density through a controlled and desired process.

III.4.3 Simple illustrations

To study the effi of the suggested method, a comparison between the LBM approach and the SciPy library was conducted. This comparison confi the perfect correspondence between the methods, for both dilation and erosion.

To illustrate these results, 30 dilations with a structuring element of size 3 × 3 were performed on Lena's image and the output is represented on the In accordance with the result stated by the proven equivalence, there is no difference between the classical dilation obtained from the library SciPy and that obtained with the LB3M. This is confi by a zero quadratic error. In order to illustrate the results of the comparison around erosion, 30 morphological erosions with a structuring element of size 3 × 3 were also performed on the image of Le na (see Figure III.5).

Once again, the erosion obtained with the LB3M and the one resulting from the library SciPy are identical, which is characterized by a zero quadratic error.

III.4.4 Region growth in grey-level (without biological interactions)

After these introductory examples on Lena's image and in order to show the use of mathematical morphology related to modelling with the LBM, an example of bone growth modelling is suggested. Indeed, a major advantage of the LB3M is that it can easily couple mathematical morphology and numerical simulation such as fluid dynamics.

For this purpose, a bone growth study is carried out. This occurs when some cells, transported by a biological fluid, succeed to anchor themselves to a scaffold. This fixing operates according to certain physical parameters such as the shear stress at the solid-fluid boundary, not taken into account in this illustration. Once the maximum cell density is reached, a new extracellular matrix layer is created on the previous one, creating a region growth [START_REF] Bidan | A Three-Dimensional Model for Tissue Deposition on Complex Surfaces[END_REF]. This growth is simulated using the following morphological equation

ρ(x, t + ∆t) = ρ(x, t) + α ((ρ⊕ S f ) (x) -ρ) . (III.14)
where α is a selected setting parameter equal to 0.5. To ensure that growth does not occur until the maximum density is reached, it is possible to use the parameter ω, previously assumed to be equal to 1. Indeed, this parameter allows a weighting of the morphological operation. So here it is suggested to use the following function as a choice (illustrative, as there are other solutions) for the parameter it is not dimensioned; however, it is quite possible to link biological parameters to LBM parameters to simulate progressive densification.

( 1 if I = 1.0 ω(I) = 0 otherwise (III.

III.5 Coupling concept and applications

III.5.1 A new concept

In order to understand the new coming concept it is useful to sum up the ideas encountered in the previous sections. In the section III.2 an image is used as a support to perform a CFD simulation with the LBM, in section III.3 and also in section III.4 some image processing is done with the LBM, then in the section III.4.4 biological growth is simulated. Therefore, inspired by [START_REF] Noël | Lattice Boltzmann Method for Modelling of Biological Phenomena[END_REF] a new concept emanates: the use of the same LBM network to perform at the same time the image processing coupled with the CFD coupled with the mechanobiology. This concept is schematized on Then each of these distributions can evolve separately or coupled at mesoscopic or macroscopic level.

In order to demonstrate the feasibility and efficient of the suggested concept, its application to a giant onionskin cerebral aneurysm is proposed.

III.5.2 Concept application to a giant onionskin cerebral aneurysm

To prove the effi iency of the suggested method, more precision about the technic used for this application is required. This is the purpose of this section 1 . The details about the image processing is evoked here, then the CFD model and the biological mechanisms.

A cerebral aneurysm is a vascular disorder characterized by a dilation of the vessel wall itself caused by vessel walls weakening. A cerebral aneurysm healing is conditioned by the formation of a clot or thrombus. The modelling of this biomechanical phenomenon named thrombosis is complex since the hemodynamic must be coupled to a biomechanical model in order to understand the thrombus formation based on patient biological factors

The segmentation of the medical images is obtained using a discrete (pixelic) levelset method performed with the LBM already published [START_REF] Chen | Segmentation of the Thrombus of Giant Intracranial Aneurysms from CT Angiography Scans with Lattice Boltzmann Method[END_REF]. Thanks to advance s in medical image techniques and reconstruction tools, patient-specifi geometries of ane urysms with parent blood vessels, are provided with the LBM (see Figure III.8).

The CFD foundations and general aspects are recalled in the first chapter of the manuscript. The hemodynamics field has been studied using CFD with LBM formulation and several studies demonstrated the LBM ability for providing a detailed analysis of the blood flow for patient-specific cases [SPC09; Ceb+05; Zha15] (see Figure III.9).

The biological mechanism of thrombus formation inside aneurysms is complex and deals with different nature of biological parameters with their variability related to patients [START_REF] Ouared | Thrombosis Modeling in Intracranial Aneurysms: A Lattice Boltzmann Numerical Algorithm[END_REF].

1 In this section the images and resources used, are from the European Thrombus project (http://www.thrombus-vph.eu) Thrombus with onionskin structure have been observed in some giant aneurysms from patients. They have the capacity to initiate the thrombus formation when the WSS i nsi de the cavity of the aneurysm is low enough to provide the right environment for the buildi ng of a blood clot [START_REF] Sagar | Dynamic Modeling of the Human Coagulation Cascade Using Reduced Order Effective Kinetic Models[END_REF]. As a result, the thrombus formation can be triggered or not, by changing the magnitude of the velocity profile [START_REF] Zhang | Hemodynamic Investigation and Thrombosis Modeling of Intracranial Aneurysms[END_REF].

In consequence, a thrombosis model has been implemented with the LBM taking in consideration the successive steps given by eq. (III.16) reflecting the biological schematic on In the above model for example, the thrombin particles cannot be generated by a thrombus node until this node becomes a new fresh wall node after a certain time to control this process [START_REF] Malaspinas | A Spatio-Temporal Model for Spontaneous Thrombus Formation in Cerebral Aneurysms[END_REF]. Since this kind of thrombosis are growing with an onionskin structure [START_REF] Zhang | Hemodynamic Investigation and Thrombosis Modeling of Intracranial Aneurysms[END_REF], when no thrombin particle stays inside the vessel, the proce ss of constructing the next layer is initiated. This is far from being usual and very different from well-studied cases [Tos+13; Oua+08; Mal+16].

III.6 Conclusion

This chapter illustrates the ability of the LBM to perform simulation directly from images which is unusual and original. Then, it reminds that the LBM can also run simulations for image processing, and goes deeper by demonstrating the equivalence between the LBM and the MM elementary operations. This proof leads to the LB3M which is then used for image processing and biological growth.

The fact of solving image processing, CFD and mechanobiology on the same support is original, but above all it has the strong advantage to allow crossing information at any time. This chapter provides the opportunity to imagine the modelling of physico-chemical and biological mechanisms by implementing the Lattice Boltzmann method on a single lattice solving different particular mechanisms.

The effi of these ideas is illustrated with the thrombus formation, but this method is given without any specifi focusing. The challenge for simulating the thrombus formation remains in resolving the biomechanical interactions between different components of the blood taking in consideration the characteristics of the blood flow in the parent blood vessel and in the cavity of the cerebral aneurysm. The presented method gives good results validated by the theoretical frame of CFD and image processing with LB3M. Plus, the mechanobiology results were validated by the medical staff of the Thrombus2 project. The implementation of LBM is a real opportunity to propose an unprecedented and effi t numerical simulation of an onion-skin structured thrombus formation by using the same network/lattice from image processing applied to DICOM medical images, to the modelling of the thrombosis mechanism coupled to the hemodynamics of the blood; leading to a novel investigation tool for medical scenarii without a priori models, based on experimental data.

Despite the promising results brought by the cerebral aneurysm, one non-negligible mechanism is not taken into account: the deformation of the blood vessel with the pressure. This deformation implies capacity to do simulations of hyperelastic solids with the LBM. These simulations are far from being straightforward and give rise to the next two chapters.

Chapter

Abstract of the chapter

The aim of the chapter is to simulate solids with the LBM. This implies to obtain generic stress tensor according to the constitutive equations chosen, which is more complex than the stress tensor given by the BGK operator. So, the MRT approach is regarded, and an analytical writing of the MRT is constructed. This leads to the GOEDF approach. Therefore, after giving a formal definition of the GOEDF, the study of these functions is necessary before it can be considered as a surrogate of the MRT. Moreover, properties of the quasi-equilibrium functions and their links with GOEDF are discussed. Then, the evaluation of the GOEDF to simulate the solids state is numerically investigated.

Résumé du chapitre

Le but de ce chapitre est de simuler des solides avec la LBM. Cela implique d'obtenir un tenseur de contrainte générique en fonction des lois de comportements choisies, ce qui est plus complexe que le tenseur de contrainte fourni par l'opérateur BGK. Ainsi, l'approche MRT est envisagée, et une écriture analytique de la MRT est construite. Cela mène à l'approche GOEDF. Par conséquent, après avoir donné une définition formelle d'une GOEDF, l'étude de ces fonctions est nécessaire avant de pouvoir être considérée comme un substitut de la MRT. De plus, les propriétés des fonctions de quasi-équilibre et leurs liens avec GOEDF sont discutées. Ensuite, l'évaluation du GOEDF pour simuler l'état des solides est examinée numériquement.

IV.1 Introduction

As it has been evoked, a major difficulty with the LBM is its capacity to cope with solids. These difficulties emanate from the theory. Indeed, in the solids the interactions between the atoms are strong compared to their mutual inertia and so cannot be neglected as it is the case in the kinetic theory. This leads to a great complexifi of the collisioninteraction operator. This complexity drives an unsolved problem, the rigorous path from the Boltzmann equation and its collision-interaction operator to the macroscopic solids equations.

One solution could be to consider that there is no collision but only forces between particles, and include these forces in the external forces term. This solution will be studied in the next chapter.

To avoid this complexity, one can try to guess the resulting form of the BGK linearization of the collision-interaction operator. This is equivalent to guess the form of the equilibrium distribution. The equilibrium is partly based on the molecular chaos. So, the molecular chaos with no interatomic forces is kept in this chapter (which is a very discussable hypothesis but helps); and some mathematical constrains are added to the equilibrium in order to obtain the correct equations of motion for solids. These constrains remind the idea of MRT, by looking at the moments rather than the distributions. Thus, the problem can be set the other way around: how to create an analytical MRT ope rator? The answer to the question might allow the impose on the distributions of an internal stress corresponding a solid stress. This is the philosophy adopted in this chapter. So, the corresponding procedure is summed up on With a view to construct an analytical MRT operator in a first attempt, and to mimic the solid state in a second attempt, it is necessary to get new settable parameters. Here, the manner adopted to accomplish these objectives is to generalize the "equilibrium" distribution. This generalization is then linked to the MRT method. A potential opening on the quasi-equilibrium method will be evoked before the final section. The final section will finally study the application of this distribution to simulate solids behaviour.

IV.2 Study of the generalized out-equilibrium distribution

IV.2.1 Choice of the distribution

To build an analytical MRT it is necessary to look at moments. From a generic point of view, f is a distribution that has many different moments; and once integrated the Boltzmann equation it reveals these moments. Thus, if firsts moments of the distribution f are given by

ρ = ρv = σ = D r f dξ (IV.1) R r ξf dξ (IV.2) R r (ξ -v) ⊗ (ξ -v) f dξ (IV.3) R r 2 E θ = kBθ = 2 (ξ -v) R r f dξ (IV.4) q θ = (ξ -v) 2 (ξ -v) f dξ (IV.5)
∂t + ∇x • (ρv ⊗ v + σ) = ρg (IV.6b) ∂E k + E θ + ∂t x • ((E k + E θ )v + σ • v + q θ ) = ρg • v (IV.6c)
which is the wanted equations for the solids if the stress σ is describing a solid behaviour. Therefore, in a general case, it is needed to have a distribution which possesses at least of 4 free parameters (ρ, v, σ, q θ , since θ = Tr(σ)). So, let us introduce the following definition:

Definition 1. a Generalized Out-Equilibrium Distribution Function (GOEDF) at the order n, noted f � , as a distribution such that the n first moments are free parameters and able to reproduce the Gaussian distribution.

Since the space of distribution is infinite-dimensional, there is an infinite number of distributions respecting the first criteria of the definition. But, since the Maxwell-Boltzmann distribution (which is a Gaussian distribution) is the unique one that describes the perfect gases, it seems reasonable to look for the distribution that becomes the Maxwell-Boltzmann distribution when the stress tensor is a constant multiplied by the identity and the heat flux is zero. Thus, a generalization of the normal Gaussian distribution represents a viable option.

As suggested, f � at third order is exhibited. Such a distribution, not uniquely defi can be expressed in function of a multivariate skew-normal distribution introduced by Azzalini [START_REF] Azzalini | The Skew-Normal Distribution and Related Multivariate Families*[END_REF]. A skew-normal distribution expresses as:

S N (ξ) = ρ ( 1 + erf ( α T λ -1 (ξ 2 -µ) ( 1 exp (ξ µ) T Λ -1 (ξ \ µ) (IV.7) (2π) D/2 (det Λ) 1/2 -2 - -
where Λ is a definite positive matrix, λ is the diagonal standard deviation matrix of Λ given by λij = δij / Λij , and α is a vector that drives the skewness of the distribution.

So for now, the conjecture that a f � at third order is proportional to a multivariate skew-normal distribution is made.

f � (ξ) ∝ S N (ξ) (IV.8)
The multivariate skew-normal distribution can be expressed in terms of physical variable as the next section shows it; it leads to the generalized out-equilibrium distribution function. This work is necessary to verify that the distribution is adequate with the desired equation system. This verification requires to compute at least the third order moment. These moments are also necessary to compute a discretized version of the distribution that can be used in a LBM algorithm, which is done in the section IV.2.3. The skew multivariate normal distribution yields the wanted moments as the following discussion explain it.

√ √

IV.2.2 Study of the moments yield by a generalized out-equilibrium distribution function

The multivariate skew normal distribution has been introduced by Azzalini [Azz05], as a generalization of the normal distribution. To study moments, several methods exist. However, the same author suggested that its moments can be studied thanks to its Moment-Generating Function (MGF) which reads

MG(t) = ρ exp ( 1 µ • t + 2 t • Λ • t \ ( 1 + erf ( δ • λ -1 \\ √ 2 t (IV.9)
where δ = λ -1 Λλ -1 α 1+αΛα ment leads to :

and

Φ(x) = ( 1 + erf ( x .
Therefore, the consecutive order mo-

2 M 0(0) = MG(t = 0) = ρ (IV.10) ( M 1(0) = ∇t M 0(0) = ρ µ + 2 \ λδ π (IV.11) M2(0) = ∇tM1(0) = ρ ( µ ⊗ µ + ( 2 \ P 2 (λδ ⊗ µ) + Λ (IV.12) M3(0) = ∇tM2(0) = ρ µ ⊗ µ ⊗ µ + P 3 (Λ ⊗ µ) 2 \ (IV.13)
+ π P 3 (µ ⊗ µ ⊗ λδ + Λ ⊗ λδ) -λδ ⊗ λδ ⊗ λδ with P ermut[n]
• is the sum of permutation index containing n terms.

It appears clear that the parameters ρ, µ, Λ, δ which depends f � are independent, so are the generated moments. A little algebra allows to recover the following centred mome nts for the distribution f � :

r f � dξ = ρ (IV.14) r ( ξf � dξ = ρ µ + r ( 2
\ λδ π 2 = ρv (IV.15) \ c ⊗ cf � dξ = ρ Λ -π λδ ⊗ λδ = σ (IV.16) r 2 ( 4 \ c ⊗ c ⊗ cf � dξ = ρ π π -1 λδ ⊗ λδ ⊗ λδ = Q (IV.
λδ = ( 2 \ -1/3 ( 4 \ -2/3 q -1 (IV.18) π π ( 2 \ 1/6 ( 4 lql 1/3 \ -2/3 q µ =v -π π -1 lq l 1/3 (IV.19) ( 2 \ 1/6 ( 4 \ -4/3 q Λ =σ - α T λ -1 = √ π δ T λΛ -1 π -1 lq l 2/3 (IV.20) (IV.21) 1 -δ T λΛ -1 λδ
it is possible to define the GOEDF at third order that satisfies the desired moments and reads:

f � (x, ξ, t) = f � ρ,v,σ,q θ (x, ξ, t) ρ = (2π) D/2 ( det σ -2 1/6 ( 4 -1 ) -4/3 q l D/2 π π q 2/3 ( \ T -1 1 ( 2 \ 1/6 ( 4 \ -2/3 q 2 1/6 ( 4 \ -4/3 q exp -2 ( ξ -v -π 2 1/6 ( 4 π -1 \ -2/3 q σ lq 1/3 π \\ π -1 lql 2/3 ξ -v -π ( ( 1
π -1 δ T λΛ -1 ( lql 1/3 2 1/6 ( 4 \ -2/3 q \\� × 1 + erf √ 2 √ 1 -δ T λΛ -1 λδ ξ -v -π π -1 lq l 1/3 (IV.22)
One can verify thanks to the equations eqs. (IV.14) to (IV.17) that the Boltzmann equation is still yielding the eq. (IV.6) when the collision-interaction operator is nullified, but for any σ and q θ settable With regards about how would be integrated this generalized out-equilibrium distribution function in the algorithm, it needs to be discretized to be used effi tly. This is the purpose of the next section.

IV.2.3 Discretization

Since the solution presented is different from the equilibrium distribution function, to obtain a numerical method, one only have to project the GOEDF on a Hermite basis to apply a standard lattice Boltzmann method.

So, by reminding that from eq. (I.73) on page 24 a discretized distribution function fi, in its general form on the Hermite basis truncated at the N -th order is

f (x, t) = w i i w(ξi) N , n=0 1 n!c 2 a (n) (x, t)H (n) (ξi) (IV.23)
with c 2 the celerity of the sound of the lattice. The expression of the coordinates a (n) of the GOEDF can be obtained easily from the moment yield by f � (see eqs. (IV.14) to (IV.17))

f � (0) s s ξ s s s √ f � (0) (q ) 2 +
Thus, the generalized out-equilibrium distribution function from the equation (IV.7) truncated at the second order moment is

( i (x, t) = wiρ 1 + v • ξi 2 s (v • ξi) 2 + 2c 4 + σ : ξi ⊗ ξi 2c 4 2 i -2c 2 v 2 -2c 2 Tr(σ) D \ -2c 2 + 2 (IV.24)
where common choice is cs = 1

(3)

in the lattice units.

And the generalized out-equilibrium distribution function truncated at the third order moment is

r i (x, t) =wiρ 1 + ξ i • v c s 2 + ξ i ⊗ ξ i :σ 2cs 4 (ξ i • v) 2 2cs 4 ξ i 2 -2cs 2 Tr(σ) -2cs 2 v 2 D -2cs 2 2 + ξ i ⊗ ξ i ⊗ ξ i :Q + 6cs 6 P 3 (ξ i • vξ i ⊗ ξ i :σ) 6cs 6 (ξ i • v) 3 6cs 6 P 3 ( (ξ i ) 2 ξ i • v ) - 6cs 4 P 3 (ξ i • q θ ) - 6c s 4 - P 3 (ξ i • σ • v) 6cs 4 - P 3 ( (v) 2 ξ i • v ) 6cs 4 DP 3 (ξ i • v) 6cs 2 (IV.25)
where Q = q θ ⊗q θ ⊗q θ .

θ Now that the generalized out-equilibrium distribution function is ready to be used in a LBM algorithm, one tries to emulate the MRT method through the generalized outequilibrium distribution function.

IV.3 Emulation of MRT

IV.3.1 MRT theory

From the Chapman-Enskog expansion, it is clear that the BGK approximation restrictively links the values of the thermal conductivity to the viscous coeffi t. This is equivalent to say the Prandtl number is not variable. The most spread technique is to resort to two distributions: one for the fluid and one for the temperature. But another numerical solution that is only based on one distribution (with some drawbacks) is the use of a generalization of the BGK operator. This latter is called MRT [START_REF] Humieres | Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions[END_REF]; it uses matrix formulation to introduce a different relaxation per moment. Indeed, the idea behind the MRT is to relax to the equilibrium in the moment space through the different times for each moment.

The discretized distribution f i can be stored in a vector of length q, the number of speed. And for example, the momentum along the axis x in D2q9 is given by

8 ρvx = , ξ i x f i = f 1 -f 3 + f 5 + f 6 -f 7 -f 8 = [0, 1, 0, -1, 1, 1, -1, -1] • f. (IV.26) i=0
But more generally, for a moment m k , the sum can be turned into a scalar product between f and a vector linked to m k , noted MM k . The m k can be computed from the equation

q m k = , MM k i f i m = MM • f (IV.27) i=0 c + + + Chapter IV. Generalized Out-Equilibrium Distribution 76 (0)
� Hence, the matrix MM is a distribution-moment transformation matrix that turns the q elements composing the discretized distribution f i (from the distribution space) into the q moments m (from the moment space). Then, assuming that MM is invertible, the BGK collision operation can be written

f (x + ξ i ∆t, t + ∆t) -f (x, t) = -ω ( f (x, t) -f (0) (x, t) ∆t = -MM -1 MM ω ( f -f (0) ∆t = -MM -1 ω ( MM f -MM f (0) ∆t = -MM -1 ωI ( m -m (0) ∆t = -MM -1 SM ( m -m (0) ∆t (IV.28)
where SM = diag(ω, ..., ω) is the diagonal matrix containing, the relaxation time, and m (0) is the vector build from the moments yield by the equilibrium distribution function. The MRT technique is the release SM from the value Iω to a relaxation time per moment through diag(ω0, ..., ωq-1). Indeed, the fact that all the elements in the diagonal matrix SM (which is now called multiple relaxation time matrix) allow to relax to the equilibrium at a different time for each moment.

To sum up, the idea underneath the MRT is to relax moments rather than the distributions directly. One obvious counterpart of this theory is that its limitation to numerical approach, since the matrix is constructed on the discretized distributions.

IV.3.2 Analytical MRT via the GOED function

A similar idea could be used, but in analytic formulation this would imply to the use infinite-dimensional matrix to obtain moments, which is not very convenient. Therefore, a workaround is suggested in the following developments.

As in the MRT approach, a distribution which yields free moments is needed. For the usual NSF equations, the moments up to the third order are required.

Proposition 3. The introduction of the GOEDF at third order in the collision-interaction operator allows to release analytically the Prandtl number, as the MRT does numerically.

Proof. Once the generalized out-equilibrium distribution at third order is defined, it can be used to obtain a continuous MRT operator. Indeed, by reminding the results of the moments obtained from the collision operator given in the equation (eq. (IV.14)); a natural construction of the collision operator would be to include the relaxation time directly inside the distribution, and would read:

Ω � (f ) = ( f � (ρ, v, ωP P (0) , T, 0) -f � (ρ, v, ωP P, T, ωqq) (IV.29)
where f � (and respectively f (0) ) is expressed with the moments of the current distribution f (and respectively f (0) ) multiplied by their relaxation times.

f

Immediately, the moments yield by the operator Ω � are:

r Ω � (f ) dξ =ρ -ρ = 0 (IV.30a) R D r ξΩ � (f ) dξ = = ρv -ρv = 0 (IV.30b) R r c ⊗ cΩ � (f ) dξ =ωP (P (0) -P ) = -ωP P (1) (IV.30c) R r c 2 Ω � (f ) dξ =0 (IV.30d) R r c 2 cΩ � (f ) dξ =0 -ωqq (1) (IV.30e) R
Thus, it appears clear that this complexifi of the collision operator allows to obtain a continuous MRT operator. In other words, with Ω � (f ) the Prandlt number is released from the fixed value of 0.7.

The last point to verify would be the variation of entropy.

IV.4 Entropic approach of quasi-equilibriums

Before to evaluate the suggested method, a discussion about the quasi-equilibrium is carried. Indeed, similarities and differences can be noted, which leads questions.

IV.4.1 Quasi-equilibrium theory

Introduced by Gorban and Karlin [GK94; GK06] in the middle 90 the quasi-equilibrium distribution functions f * are defined as the solution of the problem

min{H(f )} s.t. M (f ) = 0 (IV.31)
where H is the H-functional from the Boltzmann H-theorem, and M is a given operator link to some constrains bring by the studied problem. In the general case the f * can be determined through a convex optimization with the Lagrange multipliers.

From [START_REF] Ansumali | Quasi-Equilibrium Lattice Boltzmann Method[END_REF] the Quasi-Equilibrium (QE) collision operator reads

Ω(f, f ) = -ω * (f -f * ) -ω0 ( f * -f (0) (IV.32)
and the entropy production must remain non-positive semi-definite, which imposed the relaxation time hierarchy

ω * ≥ ω0 (IV.33)
This last relation is the reason why in recent papers [START_REF] Frapolli | Simulations of Heated Bluff-Bodies with the Multi-Speed Entropic Lattice Boltzmann Method[END_REF], distinguished cases have to be introduced for the different Prandtl regime. Another remark about the QE distributions is that under the usual constrains of imposed moments up the order n for the distribution, a classical result about the maximal entropy distribution family with a Lagrange multipliers optimization yields that the distribution should have the form

f * (ξ) = exp (Pn(ξ)) (IV.34)
where Pn is a n-order polynomial function. Unfortunately, this optimization fails for odd n numbers, since the polynomials reaches +∞ no matter its leading coeffi t.

c

IV.4.2 Link with GOED function and questions

Since the QE is requiring an extra condition about minimum of entropy, it seems clear that the QE are specific cases of GOEDF. But questions arise, with the multivariate skewnormal distribution.

A first question could be: is this distribution minimizing the entropy under the constrains over the moments expressed by the equations eqs. (IV.14) to (IV.17)? The answer to this question is relatively easy since it is known that the minimization of the entropy fails to the odd numbers, so multivariate skew-normal distribution is not minimizing the entropy with three moments. However, a finer question emanates: is the multivariate skew-normal distribution minimizing the entropy for a large number (maybe infinite) of moments?

This would prove that the multivariate skew-normal distribution is a good QE distribution. And this would have the advantage of avoiding the optimization process through the Lagrange multipliers, plus the distinguishing for Prandtl numbers could also be avoided by using the analytical MRT introduced in the section IV.3.2.

A second noticeable aspect of the GOEDF appears with the quantification of the perturbations. Trying the determine the first order distribution of the GOEDF, all the terms different from the Gaussian can be assumed as a first order. This assumption leads to

µ = v -Ek 1q θ (IV.35) Λ = pI + Ek 2q θ ⊗ q θ (IV.36) α T λ -1 = E k 3 q θ .
(IV.37) p Using the Sherman-Morison-Woodbury formula of inversion of a matrix sum and the derivation of the determinant operator, the first order Taylor expansion gives

f � = f (0) r 1 + E ( cσ -1 c 2 - k1c • q θ p k2σ -1 q θ q θ σ -1 -2 + 2k2q θ σq θ Tr(σ) k2 Tr(q θ q θ ) 2 q θ • c \ 2 ) -2p - + k3 2p π p + O ( E (IV.38)
But the moments yield by this approximation do not fit with the assumption. This could signify that the perturbation orders are more complicated to capture with a GOEDF approach.

IV.5 Solid Equilibrium

IV.5.1 Quasi-static state with the GOEDF

According to the previous sections, it appears that the GOEDF can be tuned to describe a quasi-static equilibrium for solids. Indeed, an equilibrium described by the GOEDF, according to its definition, has moments given by eqs. (IV.14) to (IV.17) on page 73. Plus, considering the following collision-interaction operation

Ω(f ) = ω ( f -f � (IV.39)
and taking ω = 1, it implies that the distribution f is going the faster to the distribution f � . This also means that if the distribution f is equal to the distribution f � , then the previous c ∇ collision operator is null and the Boltzmann equation yields the following macroscopic equations

∂ρv ∂ρ ∂t + ∇x • (ρv) = 0 (IV.40a) ∂t + ∇x • (ρv ⊗ v + σ) = ρg (IV.40b) ∂E k + E θ + ∂t x • ((E k + E θ )v + σ • v + q θ ) = ρg • v (IV.40c)
In order to simulate such situation, it is necessary to ensure that the two distributions f � and f are equal. To do so, a Gauss-Seidel procedure can be employed. But in a first approach, this optimization procedure will be skipped. This leads to differences in the evolution equations. Therefore, in first approach only the steady state is regarded, when the distributions are equals. Thus, one can rather expect to observe the following equations

∇x • (ρv) = 0 (IV.41a) ∇x • (ρv ⊗ v + σ) = ρg (IV.41b) ∇x • ((E k + E θ )v + σ • v + q θ ) = ρg • v (IV.41c)
which describes a quasi-static state.

To evaluate the adequacy of the model with the reality, a numerical verifi is conducted. To simplify the model, in a first approach only isothermal problems are considered.

IV.5.2 1D spring test case

In order to validate the suggested solution, the example of a linear continuous spri ng in 1D is studied. In both compression and decompression, the numerical method is compare d with analytical formulation. So, the following constitutive law is used to describe the l i ne ar elasticity of the matter σ = k∇x u (IV.42)

where σ is the Cauchy stress tensor, k is the scalar rigidity and u is the displacement. Therefore, the eq. (IV.6) with the relation ∂tu = v becomes

∂ t (ρ) + ∇x • (ρv) = 0 (IV.43a) ∂ 2 2 t (ρu) + ∇x • (ρv ⊗ v) + k∇ x (u) = 0 (IV.43b)
In 1D it is easy to follow the displacement using the discretized integration u (t + δt) = u(t) + δtv. Moreover, the computation of σ needed that the equilibrium distri bution function be obtained with a classical second order centered finite difference schema.

Only the steady state is studied. So, only the macroscopic variables ρ and u(x) are relevant to the steady state, and completely determined by the initial state and boundary conditions.

IV.5.2.1 Compression

The first test case is a compression of a spring. See a density of 1, the speed and displacement are null. At the initial time the left borde r is pushed to impose constant velocity and so a growing displacement occurs until a displacement of 5 length units is reached, while the right border is kept fixed with a bounce-back condition [START_REF] Malaspinas | General Regularized Boundary Condition for Multi-Speed Lattice Boltzmann Models[END_REF]. When the wanted displacement is reached by the left border, the bounce-back condition is also imposed on left border to maintain the compression. 3 show the resulting ρ and u after 25 000 time steps which can be considered as a steady state for the system since the maximum velocity after such time is lower than 0.000 35 in the units of the lattice.

Convinced by the results obtained, the decompression is also studied.

IV.5.2.2 Decompression

The decompression stage takes as initial condition the exact final step of the compressi on. The left border is released and a constant negative velocity is imposed until the zero displacement is obtained from the left border; after what the fixed condition imposed is again a bounce-back. During this time, no change to the right border conditions is performed. So, the results of the 1D test case are in agreement with the theory for compression and decompression of spring in static situations.

IV.5.3 2D test case

To illustrate the ability of the suggested method to simulate quasi-static state for any constitutive laws in 2D, the study of pure shear deformation is carried.

A square Eulerian grid is simulated. The bottom of the square is fixed while the top of the square has an imposed displacement during a given period. Then it is al so fi xe d. Both sides are free to move with periodic boundary conditions. The Figure IV.5 on the next page shows a schematic of the simulation suggested.

So, the following constitutive law, used to describe the isotropic linear elasticity of the matter under small displacements, is

σ = λ Tr (ε) I + 2µε (IV.46)
where ε is the small deformation strain tensor. The quasi-static condition over the equations eq. (IV.43) on page 79 becomes

λ ( ∂ 2 (ux) + ∂ ( ∂ (ux) )) + µ ( 2∂ 2 (ux) + ∂ 2 (ux) + ∂ ( ∂ (uy ) )) = 0 (IV.47a) x x y x y x y λ ( ∂ 2 (ux) + ∂ x ( y (ux) )) + µ ( ∂ 2 (uy ) + ∂ x ( y (ux) )) + 2∂ 2 (uy ) = 0 (IV.47b)
Plus, symmetry in the problem leads to the conditions ∂ y (uy ) = 0 and ∂ y (ux) = 0. After 20 000 time iterations, a quasi-static state is obtained. Figure IV.6a represents in colour map the L 2 -norm of the displacement resulting. This displacement is only carri e d by the transversal axis conformity to a pure shear situation. The dashed lines on Figure IV.6b represent the displacement in the y direction along the x-axis, and shows that whole the displacement is carried by the y-direction. On the same figure, the red lines are the theoretical The results seem in good agreements with the theory also for shear displacement in 2D.

x u * y

x

IV.6 Conclusion of the chapter

This chapter introduced the notion of Generalized Out-Equilibrium Distribution Function, in order to obtain arbitrary stress tensor and heat flux. After studying a multivariate skew normal distribution, it appears that it can be turned into a GOEDF. This given distribution provides a convenient example of GOEDF. These functions are then proved to be able to provide an analytical MRT operator. This theory is also numerically tested on simple examples. Through examples in 1D and 2D it is shown how these functions can simulation any constitutive law in quasi-static state with promising results.

If it brings part of answers on the ability of the LBM is simulated solids and theoretical foundations of the MRT, many question remains yet. The ability to reproduce solid dynamics through this approach is not very clear, and would deserve more investigations. A crucial part of these investigations should focus on the boundary conditions. Also, the links between the GOEDF and the QE still have the answers to give.

Despite the results yield by the GOEDF approach, other methods can be imagined to capture the dynamic behaviour of solids.

after some reminding about the theory, the Vlasov equation is adapted to solids. Then, the model is numerically tested and compared to theoretical results in 1D and 2D.

Résumé du chapitre

L'utilisation de l'opérateur de collision est totalement remise en question dans ce chapitre. En effet, l'opérateur de collision est une partie primordiale de la LBM. Cependant, il peut être contesté lorsqu'il s'agit de solides. Aussi, la suppression de ce terme dans l'équation Boltzmann conduit à l'équation Vlasov. Sur cette base, une approche originale est développée dans ce chapitre. Ainsi, cette dernière équation, sous certaines hypothèses, peut être ajustée pour s'adapter aux exigences de la simulation des solides. Ainsi, après quelques rappels sur la théorie, l'équation Vlasov est adaptée aux solides. Ensuite, le modèle est testé numériquement et comparé aux résultats théoriques en 1D et 2D.

V.1 Introduction

As the previous chapter states it, the simulation of solids with the LBM is not straightforward. Indeed, even if the use of QE-GOEDF enables to simulate static mechanical equilibrium, the dynamics of the problem are somehow not accessible with these methods. This is probably due to the manner of treating the problem as a fluid, i.e. keeping the assumption of molecular chaos.

A noticeable attempt of solid simulation with LBM done by Chopard and his team [START_REF] Marconi | A Lattice Boltzmann Model for a Solid Body[END_REF], gives dynamic results for solids. Nevertheless, their approach lies on a numerical change of the distribution of density by a distribution of forces, but the drawback is the loss of the theoretical relation with the Boltzmann equation. Indeed, the Boltzmann equation is based on particle movement and not on forces. In practice, the method suggested by Chopard is closer to the lattice-spring method than the LBM. Thus, the Chopard approach is not able to obtain all the macroscopic moments yield by the LBM. In particular, the continuity equation is not accessible, and the authors also mentioned some difficulties with the shear stress.

The method suggested by Chopard seems a direction capable of dealing with the dynamics of solids. Inspired by this approach, a similar direction but aiming to recover the theoretical link with the Boltzmann equation is attempted in this chapter.

One of the main concerns when using the Boltzmann equation to simulate solids is the use of the collision-interaction operator. Indeed, the collision operator described in the equation eq. (I.37) on page 17 is the one obtained for the model of billiard balls whose collisions would be elastic. Although questionable for fluids (in which internal forces lead to cohesion and surface tension), this model seems to be very unsuitable for solids. A classic vision of a solid material is the one of a network of masses and springs, in which each mass vibrates around an average position. This reflects a simplified monolithic crystal solid model, in which the spring force is a simplification of a more complex model of interatomic forces such as the Lennard-Jones type (see eq. (I.27) on page 13) or others, and whose vibration frequency and amplitude are related to temperature.

Nevertheless, the Boltzmann equation describing the evolution of particles density remains an interesting approach. Obviously, the best way to work with this equation would be to possess a generalized collision-interaction operator that would transcribe all interactions between the particles composing matter -elastic (non-linear) "springs" in solid phases, elastic collisions for fluids, plus the phase transitions.

Unfortunately, such an operator is not easy to exhibit, and let alone to manipulate to work with. The simple use of a potential of Lennard-Jones to model inter-particle interactions in the Boltzmann equation turns out to be a complex task which is always prone to mathematical-physical research. One of the difficulties lies in the fact that the displacement of a single particle impacts the movement of a whole set of particles in its vicinity, which makes the two-by-two study of the particles impossible.

So, a simpler approach is experienced in this chapter: use a version of the BE without collisions.

V.2 Vlasov equation and theory

In order to find a suitable framework for the use of the Boltzmann equation for solids, it is possible to note that at rest (at thermodynamic equilibrium) a particle distribution having an absolute zero temperature could be described by a Dirac distribution. This can be seen as the masses points of the lattice-spring models.

Starting from this state, a vibration of one edge of the network, will propagate freely without loss. This is reminiscent of Vlasov equation [START_REF] Vlasov | The Vibrational Properties of an Electron Gas[END_REF], which is none other than Boltzmann's equation deprived of its collision term:

∂ t (f )(x, ξ, t) + ξ • ∇x (f ) + g • ∇ ξ (f ) = 0. (V.1)
This equation is also called the optical flow equation, in some other contexts like in image processing. Then the evolution of the particles described by this equation is a propagation of the particles without further modifications. Usually this equation is used to describe collision-less plasma. In such cases, it is coupled with the Poisson equation to add the electrical interactions that drives the motion of this mixture of electrically charged particles. The Vlasov-Poisson equations are an approximation of the Vlasov-Maxwell equations and read:

∂fα + v ∂t α ∂fα • ∂x + qαE m α • ∂f α = 0, (V.2) ∂v
where qα is the particle electric charge, mα is the particle mass and E(x, t) is the selfconsistent electric field resulting from the Poisson equation. The Vlasov-Poisson equations are also useful to describe the famous Landau damping. In these phases, the collision occurs less often and the interaction are strong; which leads to a version of the Boltzmann equation where the collision-interaction operator is null and where all the interactions are put in the external forces term.

The interesting idea behind this equation is that complex systems with few collisions and strong interactions can be modelled by a Vlasov equation. This is the case of solids, even if the ergodic hypothesis is probably not true in such state of matter. The use of an adapted Vlasov equation seems a viable option to capture the dynamics of solids.

In a very dense phase like in solids, these particles are not animated by an eternal movement but rather by a movement that they pass on to their neighbours, a bit like the pendulum of Newton.

But in contrast, with the Newton cradle, at microscopic scale, in solids there is no real shocks between atoms and in order to be more accurate, one could imagine the same cradle with magnets providing momentum transmission with no contacts [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. Just as a remark this collision without any contact is also what is happening under the "elastic collision of hard spheres" in the diluted gases as the case originally considered in the Boltzmann equation.

Thus, in such circumstances, the crossing bridge from the microscopic to macroscopic scale (for any A, since Ω = 0) allows the following equations to be obtained

f (x, ξ, t) = ρδv (ξ) (V.3a) ∂ t (ρ) + ∇x • ρv = 0 (V.3b) ∂ t (ρv) + ∇x • (ρv ⊗ v) = ρg (V.3c) ∂ t (E k ) + ∇x • (E k v + κ∇x θ) = ρg • v (V.3d)
where δv (ξ) is the Dirac distribution centred on v. This system of equations is close to the one wanted at the internal stress excepted.

As in the Vlasov-Poisson equation, it is then necessary to add to this model the interaction forces. To do this, it is possible to use an equivalent of force per mass unit (or to go through the theory of mean fields). In any case, by noting that the divergence of the stress tensor acts as a mass force (this amounts to nothing less than using the Green-Ostrogradsky theorem), it is possible to write

1 g = ρ ∇x • (σ) (V.4)
Therefore, one can obtain the desired system of undefined equation motion for solids submitted to no external forces

f (x, ξ, t) = ρδv (ξ) (V.5a) 1 g = ρ ∇x • (σ) (V.5b) ∂ t (ρ) + ∇x • ρv = 0 (V.5c) ∂ t (ρv) + ∇x • (ρv ⊗ v) = ∇x • (σ) (V.5d) ∂ t (E k ) + ∇x • (E k v + κ∇x θ) = (∇x • (σ)) • v (V.5e)
Thus, the Vlasov equation leading to the previous system of equations seems suitable for the solid simulations. In order to perform numerical simulation, a discretization needs to be done. The implementation of a discretized version is done in exactly the same way as for the Boltzmann equation (cf. section I.4 on page 21). All of these concepts are synthesized in the schematic given in Figure V.1.

A remaining methodological difference with the LBM is about the volume force. Indeed, in previous cases such as those of section II.4.1 and section II.4.2 on page 46 the external forces were constant, where here it is varying with the internal stress. To hope to simulate solids with this method, it is necessary to be able to evaluate numerically the divergence of the stress tensor. To do this, a relatively simple numerical method, moreover, compatible with the LBM, is the use of finite differences.

V.3 Application to solids

To study the eff eness of the method proposed in the previous section, it is useful to resort numerical simulation. This makes possible to observe whether the method generates situations similar to those encountered by the experiment. The first study case real i ze d is that of the compression of a 1D spring modelled by a segment. In the next section, a she ar wave, and the localized compression of a 2-dimensional plate will be studied.

V.3.1 1D case: Spring compression

A first and simple test to confront the proposed model with the theory is a test made of a 1D spring whose right end is fixed and immobile; while the left end is mobile and controlled in speed. By imposing a speed in the material on the left, it indirectly imposes a compre ssi on of the spring described with an Eulerian point of view. This speed is maintained constant for a defined period of time ∆t (see Figure IV.2 on page 80). In practice, this time is chosen relatively short compared to the time required for the compression wave generated in the spring to travel the Eulerian length of the spring. This makes it possible to obse rve the dynamic behaviour of the spring in compression. After, the left border impose s a ze ro speed on its edge with bounce-back. So, when it is reached by the compression wave it allows rebounds without losses on the wave during its "way and way-back" on this edge.

Since the spring is perfect (without friction or viscosity), perpetual back and forth of the wave is expected. Nevertheless, the waveform may possibly change in such a manne r which reduces the local mechanical stress in the material, in accordance with the law of behaviour chosen for the spring.

For this first simulation, the different parameters were used:

• Eulerian length if the spring lx = 200

• time length of the impulse ∆t = 50

• linear constitutive law: ∇x • (σ) = k∇x • (ε) = k∂ 2 (u)
• dimensionless stiffness of the spring k = 0.15

The Figures V.2 to V.3 does well exhibit a round-trip behaviour of a compression wave according to theoretical predictions. Furthermore, a slight inflection of the wavefront is visible, which tends to show good behaviour of the model. A last notable aspect that reinforces the good capabilities of the model is its temporal stability. Temporal stability shows no significant defects in the 20 000 time steps used for this simulation. However, the rear part of the wave has many unrealistic oscillating disturbances, more visible on Figure V.6b (blue line). These oscillations appear to vibrate at the spatial sampling frequency and appear on the left edge of the simulation. A plausible explanation is that the schema used inappropriate finite differences. Possibly, a Weighted Essentially Non-Oscillating (WENO) schema would offer some better results. Another explanation would be a perfectible boundary condition. The last possibility would be linked to the use of volume forces whereas the stress is a surface force.

V.3.2 Surface forces

Indeed, this last explanation seems quite plausible because a finer study of the origin of these oscillations shows that at the start-up, when the incident wave deforms the material, an unbalanced force is created. This force accelerates the particles before they are directly affected by the wave, allowing a decrease in local stress upstream of the wave, which is physically correct. But this force also accelerates in the same direction the particles located just after the wavefront, which is incorrect, because one would expect them to be slowed down in order to conserve energy locally, as represented on the noise is reduced by surface forces. Unfortunately, this leads to an unstable numerical scheme in time as the values booming on Figure V.6a shows it. However, this approach brings an element of answers about the observed numerical oscillations. Although promising, this approach will not be used in the later stages, as it will require further development to achieve satisfactory stability. Data :

f 2 f 1 f 2 f 1 f 2 f 1 f 2 f 1 Chapter V.
f i s , f i (0) Result : f i c V.1.2 f i c ← f i s -ω( f i s -f i (0) + 0 V.1.3 EndFunction V.1.4 Function stream(): Data : f i c , gi Result : f i s V.1.5 f i s ← Stream(f i c ) V.1.6 g c ← (gi) V.1.7 g s ← Stream(g c ) i i V.1.8 g s ← BounceBack(g s < 0) i i V.1.9 f i s ← f i s + ( 1 -1 ω ) g s 2 i
V.1.10 EndFunction Algorithme V.1 : Pseudo-code of the streaming and collision steps of a LBM with surface forces Since only the y-axis is carrying the information, it will be considered in the colour map from now rather than its norm. Thus, going further in time, This succession of state, described by the Figures V.8 to V.10 is repeating after, leading to the same situations over tens thousands of time iterations. Thus, propagation seems to be without lost and exciting the first shear mode of the system. Even if no comparison with analytical or other numerical methods are provided, the suggested method seems to yield the expected results, in these simple cases.

V.3.3.2 Stamping

A more complicated 2D study case, but also physically more interesting, is brought here: a stamping example. Thus, the study of localized compression (stamping) offers the possibility of observing both the dynamics of the compression wave in the direction of the applied force normal to the imposed velocity and the dynamics of the shear wave in the perpendicular direction.

A scheme of this study case is given on Figure V.11a. One can observe that this Eulerian stamping test is realized by imposing a speed normal to the surface of a portion of the top surface during a given time. Then the speed is forced to zero to maintain the constrains. On all the other boundaries of the area of interest, the displacement is imposed to zero through a bounce-back condition. Since, from chapter III, all simulations are performed on images, Figure V.11b represents the input image on which the simulation done. • time length of the impulse ∆t = 50

• behaviour law: ∇x • (σ) = k1 ( ∂ 2 (ux) + ∂ 2 (ux) ) + k2 ( ∂ 2 (uy ) + ∂ 2 (uy ) ) x y x y
• dimensionless constants of the shell k1 = 0.025, k2 = 0.15.

The figures Figures V.12a to V.16 on pages 97-98 show the results obtained at different instants.

The Figure V.12a is a colour map of the L 2 -norm of the displacement. This latter is mainly carried by the x-axis, which is due to the compression. On Figure V.12b the speed norm of the impulse going down is visible. Nevertheless, on the same fi gure , two smal l shear waves rising from the side of the speed packet is very clear, because it makes the si de edges blur. After 1555 time steps, the impulse did some round-trips. Thus, the Figure V.13a shows the displacement to which the plate is subjected. The wave character of the move ment is well observable on the Figure V.13b, where the compression and shear wave are streami ng and superposing in the shell.

After some hundreds of time steps, it becomes harder to distinguish the compression and shear wave clearly, since they are diluted and going in all directions. Nevertheless, in these mixed conditions, some patterns appears periodically. The Thus, the suggested method seems to be able to capture the dynamics of the solids. However, it has to be noticed that for the last simulation after 58 300 time steps, the simulation diverges. A very likely explanation lies in the accumulation of errors, already presented for the 1D cases, even if the 2D induces a reduction of their visibility, these perturbations are still present in the 2D simulations.

V.4 Conclusion of the chapter

This chapter studied a completely different approach to tackle the solid dynamics with the LBM. The approach imagined is to remove the collision operator, so use the Vlasov equation, plus to introduce the divergence of the stress tensor as external forces. This approach has the advantage to be relatively simpler and to have better physical meaning. The suggested approach offers promising results. This Vlasov method for solids is confronted to 1D and 2D situations. The dynamic behaviour of a 1D spring is studied for both volume and surface forces. Then the shear wave in a 2D ring is simulated. Also, the dynamic evolution of a 2D shell stamping shows results predicted by the modal theory. The results given in 1D and 2D seems very close to theory, even if further investigations should be accomplished to definitively validate this model. This model should also win in numerical stability, to be more widely applicable. Of course, this stability remains largely conditioned by boundary conditions.

Presented works in a nutshell

As it is stated in this manuscript, the numerical simulation for image-based diagnostic is a demanding scientific field; demanding in terms of scientific specialities, hardware equipment and computational time. Therefore, in order to reduce these demands, the possibility to use the LBM numerical method is questioned. Indeed, the LBM appears to be able to recover the Navier-Stokes-Fourier equations from solving the Boltzmann equation at the macroscopic scale. This approach is known to be faster but require more computational memory. When combined with the fact that the LBM is able to manage image processing operations, it seems that the LBM could bring major advantages to the image-based diagnostic. Also, through this manuscript different questions are investigated, around the building of a pipeline able to perform numerical simulations for image-based diagnostics only with the LBM.

Thus, before to imagine the development related to these investigations, an implementation and validation of a new LBM are made. The implementations realized provide features such as fast and generic tools, ready for parallelization and certifying the nonregressivity of the code. This non-regressivity is also used to validate the implementations through a test cases benchmark. The test cases presented are a hydrostatic pressure profile, a Poiseuille flow under gravity and Zou-He boundary conditions, and Von Kármán vortex street. All these tests validate the initial code for further modifications.

Pursuing the idea of the pipeline from images to numerical simulation, a new concept is introduced. This concept is based on the use of a single LBM network to perform at the same time the image processing, the mechanical and biological simulation. To enhance this concept, the LB3M is introduced. This adaptation of the LBM equations is proved to be equivalent to the MM. The whole is illustrated on a giant cerebral aneurysm. This example emphasizes the powerful concepts using the LBM, but also highlights the missing element of the solid dynamics with the LBM.

Two different approaches are investigated to tackle the problem of solids with the LBM. The first one defines a GOEDF to provide an analytical MRT operator. This operator allows to control arbitrary stress tensor and heat flux in the equations while keeping a true link with the theory. So, the stress tensor is driven in numerical examples to reproduce solid behaviour with any constitutive law. An example of compression and decompression of a 1D spring is furnished. This is followed by an example of a 2D shell sheared. These two examples reveal the same conclusion, a good convergence to the quasi-static state.

The second approach is radically different and suggests using an adaptation of the Vlasov equation. This equation is a collisionless Boltzmann equation where the divergence of the stress tensor is placed in the external forces. This simpler approach keeps physical meaning. The suggested method is tested on 1D and 2D test cases. The 1D example of spring compression is reused, and with this method the dynamic behaviour is captured. A variation of the expression of forces to go from volume forces to surface forces is also tested on this example. Then, 2D examples are carried. The 2D shell shear is also reused, with the same observation about the dynamic capture made. To finish, a 2D stamping example is studied. This example allows to visualize on the same simulation the compression and shear waves. This simulation also seems to present a modal behaviour of the shell, in accordance with the elastic theory.

Contributions and answers about the image-based diagnostic

Thus, the works carried in the manuscript bring element of answers about the initial questioning. The questions about the feasibility of a pipeline from image to numerical simulation through only LBM were far from being solved. Also, questions about the use of the LBM for the image processing were asked. And many interrogations about the solid simulations with the LBM were presented.

More than a simple answer, the feasibility of a pipeline from image to simulation is demonstrated through a complete example. This demonstration brings the LBM closer than ever to the image-based diagnostic, and offers a real save of time and complexity of pre-processing.

The pipeline is enhanced by the incorporation of the MM through the LB3M. The LB3M answers about the ability of the LBM to do MM operations through a mathematical proof and numerical examples. This proved links enlarge the possibilities of the LBM for the image-based diagnostics.

A bit less directly and a bit more theoretically, the new definition of the GOEDF also do its bit of image-based diagnostic. Indeed, the GOEDF allows among other things to produce an analytical MRT, which is useful to control arbitrary stress tensor and heat flux in the LBM. This control numerically shows to be able to simulate quasi-static state of solids. So, it is obviously a partial answer to the solid simulations with LBM, but it also implies that the LBM can hope to perform numerical simulations for image-based diagnostic of solid matter in quasi-static state.

The last contribution brought by this manuscript is the use of the Vlasov equation for solids. This method shows numerically to be able to capture complex transient dynamic phenomenon and seems to obtain modal regimes. These results bring additional information about the ability of the LBM to capture solid dynamics. So, it also implies hopes for the LBM to do simulations for image-based diagnostic of solid dynamics.

Despite all the elements of answers about the LBM and the image-based diagnostic associate to any continuum media, more questions emanate.

Opened questions and future work

All contributions, brought by this manuscript work, open new questions on its way.

To extend the use of the LB3M outside of the image-based diagnostic, in pure image processing for example, more information about its computation time cost are necessary.

Romain Noël

Of course, the question of the use of other image processing methods with the LBM remains open and full of interests. Especially, the image clustering with the LBM could be used to simplify images before numerical simulations by reducing the number of colours to manage.

The GOEDF brings its pack of questions also. What could be a GOEDF at fourth order, and how to use it? The introduction of a p-norm in the exponential function could be a lead to follow, as it could a variable kurtosis. Also, questions about the ability to the third order GOEDF to capture transient solid dynamics remains entire. And the impact of boundary condition in the previous question could play an important role. Other questions about the link between the GOEDF and the QE request more mathematical investigations.

The adaptation of the Vlasov equation for solids is simpler, but not necessarily are the questions about it. Since, the adapted Vlasov equation seems to present similitude with the lattice-spring method, is it possible to find a mathematical path to link them together? This question could resort to the phonon quasi-particle to investigate the path. The same vibrational particles could probably add a thermal information to the Vlasov equation which is another remaining issue with this approach.

The introduction of higher spatio-temporal dimensions could also lead to useful applications. Indeed, the possible to solve 3D problems enables to have more realistic and sophisticated simulations. But also, in 4D, i.e. 3D spatial dimensions plus 1 temporal dimension, it is possible to challenging to perform calibration and to tackle inverse problems. In the same philosophy, a major question is the use of the adapted Vlasov equation and the optical flow to determine the stress evolution directly from a sequence of images.

For both the adaptation of the Vlasov equation and the GOEDF, is it possible to improve the stability and precision of the results? The fuse of the adapted Vlasov equation and the GOEDF could probably help in this direction, combined with an appropriate expression of the boundary conditions. And for better performances, is it possible to recover a completely local numerical method with these approaches? Such improvement would increase the parallel power of the computations with these methods.

The last but not least group of questions could concern the Fluid-Structure Interactions. In order to continue in the image-based diagnostic direction, investigation of the multiphase-multicomponent approach to manage the Fluid-Structure Interactions would be a great step. This investigations could use a formulation of interface forces at a mesoscopic scale. 106 for multi-class and heterogeneities, has been introduced here. Its ability to reproduce the fundamental diagram is proved here, for both single-class and multi-class flows. This al l ows easily simulating complex and realistic cases of mixture of multi-class traffic flow. The se simulations are able to capture jamming in various traffic situations such as road merging, reduction of number of lanes or change of speed limits.

Résumé de l'article

La modélisation du trafic routier garde souvent l'échelle mésoscopique dans la sphère théorique à cause de la nature intégro-différentielle de ses équations. Dans le travail présent, nous suggérons d'utiliser la méthode de Boltzmann sur réseau pour surmonter ces difficultés. En particulier, cette méthode a l'avantage d'avoir de solides fondations théori ques. Une version améliorée de la méthode de Boltzmann sur réseau, permettant de gérer l e s hétérogénéités et les mélanges multi-classes, est introduite ici. Sa capacité à reprodui re le diagramme fondamental est prouvée ici, pour les trafics à simple classe comme ceux à multiples classes. Ainsi, elle permet de simuler aisément des cas complexes et réaliste s de flux routier multi-classes. Les simulations sont capables de capturer des phénomènes de bouchons et ralentissements routiers dans de nombreuses situations telles que les insertions, la réduction du nombre de voies ou les changements de limitation de vitesse.

A.1 Introduction

Since the end of the twentieth century, with the increase of personal car and road jam, traffic modelling has become a subject in the centre of economic, ecological and social considerations. In this context, constructors and road operators look forward to havi ng a better understanding and anticipation of the traffic flows, in order to promulgate the optimal practical solutions.

In this perspective, the models should handle the heterogeneous characteristics of the traffic flow. One of its major heterogeneity is its composition by multiple classes (or categories) of vehicles. Road operators and car drivers agree on the role of lorries. Due to their longer dimensions and heavier weight, they are incriminated for faster deterioration of the structures, densification of traffic situations, slowing of flows and faster creation of traffic jams for a longer time.

Models have to be able to capture various nonlinear phenomena, such as those responsible for the growth of traffic jams. Commonly traffic models are classifi in three groups: microscopic, mesoscopic and macroscopic. First, microscopic models descri be the behaviour of each driver individually. So, the interaction between two vehicles can be studied finely. This level of detail allows including the psychological aspect of the dri ve rs. At the microscopic scale, the lane-changing and the overtaking questions are cruci al . The se models are very adapted to simulate the evolution of the traffic inside cities' streets but also what could be the new traffic behaviours with intelligent or unmanned vehicles. The widely used method in these approaches are the cellular automaton [CLQ96; NS92; SS93; Nag+98] and the car-following theory [Pip67; New02; New93]. The main counterpart of the microscopic scale is the numerical resources needed to simulate large areas or numerous vehicles. This is due to the number of differential equations to solve with the car-fol l owing and the repeat simulations in order to improve signal-to-noise ratio with the cellular automaton.

At the opposite of the spectra, the macroscopic models have high computational effi ciency. This effi is mainly obtained thanks to the reduced number of macroscopic variables describing the traffic flow and the nature of the partial differential equations ruling these variables [LW55; Ric56; Pay71]. Such models are the mixture results of the empirical behaviour and hydrodynamics equations. It is an appropriate model for long roads where the position and velocity of each vehicle are not at stake nor significant.

Mesoscopic models use a statistical description to recover the macroscopic equations with a finer level of detail. These approaches adapt the kinetic theory of ideal gases to traffic situations [START_REF] Prigogine | Kinetic Theory of Vehicular Traffic[END_REF]. They are based on Boltzmann-like equations, therefore they lie on the microscopic interactions, which give "physical" explanation and foundations to the resulting behaviours. On the one hand, the main mesoscopic variable is the distribution function, reducing the numbers of variables and the computation time. On the other hand, as a price to pay to have fine details and numerical effi , is a work in functional spaces plus the complexity of the Boltzmann-like equations and their solutions.

The lattice Boltzmann method (LBM) is born after the statistical averaging for lattice gas cellular automaton, in order to have meaningful results when simple simulations are subjected to numerical noise. But the LBM has a strong theoretical base lying on the Boltzmann equation [START_REF] He | Theory of the Lattice Boltzmann Method: From the Boltzmann Equation to the Lattice Boltzmann Equation[END_REF]. Thus, it seems to be a natural framework to build numerical schemas for mesoscopic models.

The multi-class eff are challenging and have been studied at the macroscopic scale [LI08; JW04]. Nevertheless, there are fewer studies at a microscopic scale with cellular automaton [START_REF] Lan | Inhomogeneous Cellular Automata Modeling for Mixed Traffic with Cars and Motorcycles[END_REF] or with car-following [Hid02; PZZ05]. The available work for multiclass traffic flow at mesoscopic scale is mainly theoretical [HB00; HB01; Cha05; SH99].

In the present work, a statistical description of continuous kinetic model is introduce d, before discussing the eff of multi-class heterogeneity on continuum models. Then the construction of the lattice Boltzmann methods and the formulations associated, for both the single class and heterogeneous multi-class frameworks. This is followed by nume rical simulations in order to validate the suggested model. Finally, discussions and concl usi ons are performed in the last section.

A.2 Statistical Description

A.2.1 Continuous Kinetic Models

To build a continuous mesoscopic and kinetic model, one would study the density of vehicles ρ at a spatial position x along a highway, at a time t. Obviously, on multi-lane roads, several vehicles can be at the same time and position with a different velocity. This remark is also true, with single lane road, since the representative elementary volume (REV) used for the mesoscopic scale is such that it contains many vehicles. Naturally, this leads to introduce f (x, ξ, t) the distribution of vehicles with the velocity ξ. By working in the velocity space, the density, average speed v and flow q are trivially given by: r

ρ(x, t) = f (x, ξ, t) dξ (A.1) r q(x, t) = ρ(x, t)v(x, t) = ξf (x, ξ, t) dξ. (A.2)
In order to describe the evolution of any distribution of particles, Boltzmann suggested using the sum of a transport term with interactions between these particles. Since the ∂t ∂t -vehicle can be considered as particles interacting, the Boltzmann equation seems perfectly adapted to vehicle transport problems linked to traffic flow. Thus, Prigogine, Herman and Andrews [PA60; PH71], suggest describing the evolution with no "external forces" and by decomposing the collision-interaction operator into two terms. The Prigogine-Boltzmann equation for traffic flow is

( ∂ ∂ \ + ξ ∂t ∂x f = Ω(f, f ) = ( ∂f \ ∂t rel ( ∂f \ + ∂t int (A.3)
where the left-hand side member is the transport term, Ω is the collision-interaction operator, The relaxation term is the reckoning that drivers will toward a certain speed called desired velocity. This desired velocity for all drivers is described by a distribution function f d , and it is reached after a certain time τ called the relaxation time. Therefore, the relaxation term can be expressed as:

( ∂f \ ∂t rel f (x, ξ, t) f d (x, ξ, t) = τ . (A.4)
The interaction term renders the overtaking process. When a fast incoming vehicle attains a slower one, the slow one is not affected but if the fast one cannot overtake it has to slow down in a short period. The interaction term can be given by:

( ∂f \

∂t int r = (1 -p)f (x, ξ, t) (ξ 1 -ξ)f (x, ξ 1 , t) dξ 1 (A.5)
where p denotes the average probability to overtake slow vehicles. This model is criticised by some authors [START_REF] Munjal | An Analysis of the Boltzmann-Type Statistical Models for Multi-Lane Traffic Flow[END_REF]. Paveri-Fontana highlighted some problems with this model when τ and p are constant. He also suggests a close enhanced model [START_REF] Paveri-Fontana | On Boltzmann-like Treatments for Traffic Flow: A Critical Review of the Basic Model and an Alternative Proposal for Dilute Traffic Analysis[END_REF]. Moreover, he shows that Prigogine-Boltzmann equation and the Paveri-Fontana improved Boltzmann equation yield same macroscopic equation up to the second order.

Thus, as a first approach the Prigogine-Boltzmann equation is a simpler and good description of the traffic flow if τ and p are linked to the vehicle density. A common relationship between the parameter p and the density ρ is γ = τρ(1p), where γ is a constant. In order to have a model closer to the conventional lattice Boltzmann method, the Prigogine-Boltzmann equation can be rewritten

( ∂ ∂ \ + ξ ∂t ∂x f (x, ξ, t) = f (0) -f τ 1 (A.6) with f (0) = τ 1 f d + ( 1 τ τ 1 \ ρ(ξ τ τ v) τ 1 = . 1 + γ
In the following, the distribution f (0) is called equilibrium distribution function.

Some authors continue this path to create generic kinetic traffic model adapted to many situations [START_REF] Helbing | Traffic and Related Self-Driven Many-Particle Systems[END_REF]. Despite the discussions about the validity of their assumptions (specially about the vehicular chaos), the continuous kinetic models seem to be good tool s to deal with traffic flows. The LBM in its classical form can solve various problems related to transport of particles, under certain conditions like compressible or viscous flows. The equilibrium distribution function has to be in adequacy with the physical phenomena underlying that one desires to capture. This distribution function is, in the Euler conservation case, uniquely found through the application of mathematical theorem [START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF] and physical conservation laws.

When applying the LBM for road traffic [START_REF] Meng | Lattice Boltzmann Model for Traffic Flow[END_REF][START_REF] Shi | Revised Lattice Boltzmann Model for Traffic Flow with Equilibrium Traffic Pressure[END_REF], the lack of conservation laws makes the appreciation of the equilibrium distribution function harder. As some authors suggest [START_REF] Meng | Lattice Boltzmann Model for Traffic Flow[END_REF] a workaround by constructing it from the observable data. The other noticeable difference with the LBM applied to traffic problems is the phase space. It is R n (symmetric) for fluids or gases but it becomes R + (asymmetric) with roads since the dimension of roads is one and vehicles almost never use reverse gear for something else than to park. Therefore it is suggested to use a D1q6 (see Previous studies [START_REF] Meng | Lattice Boltzmann Model for Traffic Flow[END_REF] suggest choosing an equilibrium distribution function that offers interesting capacity to model empirical phenomena. Here a completeness of this empirical model to deal with heterogeneities is introduced which allows defining the equilibrium distribution function with

\ ∀ i � 1 ρ(x) � 0 ∀ i ∈ [vm + 1, m -1]] vm ρ(x) = � , ρ(x + ei) i=0 vm + 1 (A.20)
where vm is the maximum speed of vehicles (it can be the speed limit, like on Figure A.1, or higher values, if ones wants to capture over-speeding) and � parameter. This latter can be seen as the reachable forward occupation rate. In other words, it is the density, that the drivers feel in front of them, in which they will have to navigate. Obviously, the maximum speed of vehicles vm is an integer lower than the maximum speed used to model the system m -1. This also means, that the variations of vm along the road can only be multiple of the lattice step δx/δt = e. This equilibrium distribution function expressed the fact that drivers under a constant relaxation time, change their desired speed with the variation of the forward reachable occupation rate. It is a manner to express that they are adapting their speed with the traffic. This model has to be completed with "virtual boundary condition" to prevent having density higher than one when too many vehicles from different cells want to reach the same cell. Naturally, the vehicles from the further cells have to slow down. A virtual boundary condition inspired by [START_REF] Meng | Lattice Boltzmann Model for Traffic Flow[END_REF], can be expressed through the function Bv (see eq. (A.21)), in which the quantity ρtest (see eq. (A.22)) represents the density that would happen if all the drivers could stream as they wish, i.e. if the density could be greater than one and so that car crash could happen.

To avoid any confusion, since in the definition of Bv the comparison is made wi th 1, the definition of ρ or ρtest used here is the occupation rate varying from 0 to 1. Thus, is suggested to rewrite the streaming step that takes into account the virtual boundaries: fi(x + eiδt, t + δt) = Bv (fi(x, t + δt)).

(A.23)

A.3.2 Lattice Boltzmann Model for Multi-Class Traffic

Strengthened by its impressive results of the LBM to model the Navier-Stokes equations, researchers quickly tackled more complex problems like mixture of fluids. A noticeable work has been done in case of immiscible mixture [START_REF] Shan | Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components[END_REF]. At the same time, some pioneers initiate the work on miscible fluids [START_REF] Holme | Lattice-G as and Lattice-Boltzmann Models of Miscible Fluids[END_REF], improved some years after by adding the thermodynamics [START_REF] Inamuro | A Lattice Boltzmann Method for a Binary Miscible Fluid Mixture and Its Application to a Heat-Transfer Problem[END_REF].

When it comes to traffic, despite heavy machinery vehicles tend to form a continuous lane on highways, the second lane is most of the time full of personal car. So, the hypothesis of immiscibility seems not to be the most relevant in a first approach. Even i f, fe w publications deal with the use of cellular automatons for multi-class traffic flows [START_REF] Ez-Zahraouy | The Effect of Mixture Lengths of Vehicles on the Traffic Flow Behaviour in One-Dimensional Cellular Automaton[END_REF]; the use of the lattice Boltzmann method to solve heterogeneous multi-class traffics has never been studied.

In order to take into account the necessity for the equilibrium density function to depend on the global density and the class density, and since the variable � the density felt forward drivers; it seems natural for a first approach to link the dependency ρ. Thus, the macroscopic variable related to the eq. (A. The virtual boundary conditions also have to follow the same logic. One should notes that since the Bv function changes the values of the distribution at a point backward from a given point where variables are evaluated, its algorithmic application should be backward recursive. The variable ρtest has to be the reflection of the global density, while the modifications over densities has to be accomplished per class. Therefore, the virtual boundaries can be expressed by:

A.4 Numerical Validating Simulations

In all the following simulations, we take a space step of 5.5 metres and a time step of 1 second. This leads to a speed step of almost 20 kilometres per hour. Moreover, we use d a D1Q6 schema, consequently the maximum speed is close to 100 kilometres per hour. But for the sake of generality, the following results will be presented in their non-dimensional form, i.e. expressed in space or time step units.

All the following simulations are obtained with a relaxation frequency of 0.9, at the exception of the simulations linked to Figure A.5 where the value of the relaxation time is varying. The presented model is here confronted to different simulations: fundamental diagrams, road merging, speed limit and number lane changes and variation of truck concentration.

A.4.1 Fundamental Diagrams

In order to evaluate the behaviour of the suggested model, the study of density-flow fundamental diagram is performed, which allows estimating the relationship between the flow and the density [START_REF] Kerner | Cellular Automata Approach to Three-Phase Traffic Theory[END_REF][START_REF] Ngoduy | Multiclass First-Order Traffic Model Using Stochastic Fundamental Diagrams[END_REF]. The study of a ring road, of 1000 cell length, simulated for different average density. Each simulation starts with a density spatially varying around the average density by following a random noise of 10%. The road has a speed restriction of 5 cells per unit of time (see Figure A.3). The relaxation time is set to 0.9. After 2000 time steps, the averaging over time and space is made to obtain the Moreover, the eff of the relaxation time on the fundamental diagram is major. Thi s eff is represented on the Figure A.5. The relaxation time interval in which the numerical schema remains stable is directly linked to the equilibrium distribution function suggested. And the fact that its eff is only noticeable for congested-flow can be i nte rpre ted as the too slow ability of the drivers to change (or reach) their desire speed. In other words, the y have a speed incompatible with the current density of traffic.

For values of the relaxation time higher than 1.20 and lower than 0.65, the numerical model becomes unstable for high densities.

A.4.2 Road Merging

Changes of the traffic conditions can be caused by road merging. To study the eff of road merging on traffic conditions, the simulation of 2-lane highway is considered. This road is constituted of 5000 cells and the speed limit is 5 cells per unit of time. The complete highway is empty at the time t = 0 of the simulation. At the beginning of the road, the density is set to a constant while 2000 cells after a road merging add another constant density except if the sum is higher than one (see

A.4.3 Number of Lane

To evaluate the eff of a change of lane number on traffic conditions; the study of a road of 5000 sites is proposed. It starts with 3 lanes and a constant density, while a reduction to 2 lanes is located at site number 2500.

A.4.4 Speed Limit Change

The last external source of modification for traffic conditions studied here, is the change of speed limit. The same length of the road is used than the previous numerical investigations for a road of 2 lanes and the same relaxation frequency. As Figure A.12 represents it a reduction of speed limits from 5 to 4 cells per unit of time is imposed at site number 2500.

A.4.5 Truck Concentration

To evaluate the eff of heterogeneous multi-class traffic, a ring road is studied. The simulation of the appendix A.4.1 is reused but a variable ratio of 2 vehicle classes is put in the model. These two classes have different speed limits: one has a speed limit of 5 ce l l per unit of time, while the other has a speed limit of only 4. The ring road is made of 1000 cells, 3 lanes (see Moreover, the latter figure highlights that the maximum possible flow is lower, with slower vehicles. This is even more drastic when re-dimensioning the fundamental diagram. Indeed, if the slower vehicles are lorries since they are more represented through the eq. (A.12).

A.5 Discussions & Conclusions

The suggested method allows to reproduce macroscopic situations from a mesoscopic formulation through a Boltzmann-like equation. The present lattice-Boltzmann method is capable of modelling the ideal traffic behaviour described in the literature (such as Drake model), and the basic psychological behaviour of the drivers through the relaxation time parameter. The improved formulation allows simulating various road situations, unprecedented with such method. The method is able to deal with road merging, change of lane numbers or speed limits in both free-fl w and congested-fl w conditions in accordance with the macroscopic previsions. The influence of a multi-class traffic found here is close to those presented by some authors [START_REF] Ez-Zahraouy | The Effect of Mixture Lengths of Vehicles on the Traffic Flow Behaviour in One-Dimensional Cellular Automaton[END_REF].

To conclude, the lattice Boltzmann method is an effi ient numerical method to ove rcome the integro-differential difficulties introduced by statistical models and might be a practical manner to solve numerically the Prigogine-Boltzmann like equation. This remains a good compromise between, on the one hand, the level of detail but time consuming provided by the microscopic description and, on the other hand, the loss of informati on but the much faster computation yield by the macroscopic description. Moreover, the macroscopic results are bounded to the choice of the equilibrium density function that can be tuned to reproduce various eff and models. The present results prove for the first time the capacity of this method to solve he terogeneous multi-class traffic flows with precision. And suggested formulations gi ve e asy treatment to deal with numerous road situations. Plus, statistical aspects in the me thod could be useful for stochastic simulations.

This work could be extended by adding of new parameters involved to model cer-

ix I. 1

 1 Graph of the potential energy resulting from the attractive and repellent forces after the Lennard-Jones model. ............................................................ I.2 Sequence of the main concepts composing the Lattice Boltzmann Method (LBM)................................................................................................................... I.3 Scheme of the discretization of the velocity space. ............................................. I.4 Lattice schemes in 1D. ........................................................................................ I.5 Lattice schemes in 2D. ........................................................................................ I.6 Lattice schemes in 3D. ........................................................................................ I.7 Collision and streaming sub-step illustration. ..................................................... I.8 Graphical representation of the entropic LBM. .................................................. II.1 Organization of the computational tool............................................................... II.2 Organization of the warping environment. ......................................................... II.3 Schematic connections between warping environment and the GUI. .................. II.4 Preview of the GUI. ............................................................................................ II.5 Pressure field along a 2D column of water under gravity ................................... II.6 Norm L 2 of the speed field along a 2D flow of air under gravity........................ II.7 Graph of the speed in the direction of gravity along the orthogonal direction. 47 II.8 Norm L 2 of the speed field along a 2D flow of air under gravity........................ II.9 Graph of the speed in the direction of gravity along the orthogonal direction. 49 II.10.................................................................................. P araView visualization of the streamlines and norm of the speed for Von Kármán vortex street. ....................................................................................... III.1 Illustrative computation of blood flow through an aortic aneurysm from medical image. .................................................................................................... III.2 Schematic of the probability of streaming. .......................................................... III.3 Result of segmentation with a LBM anisotropic diffusion reaction [Wan14]. . III.4 Morphological dilation in greyscale using the LBM applied to Lena's image. III.5 Morphological erosion in greyscale using the LBM applied to Lena's image. III.6 Use of mathematical morphology in LBM to simulate bone growth from images. ................................................................................................................ III.7 Schematic representing the order to use the different LBM layers ..................... III.8 Result of the thrombus segmentation with a LBM.............................................. III.9 Computed velocity of the blood flow inside the aneurysm and its parent vessel 65 III.10Schematic of the thrombosis reaction system ..................................................... xvii List of Figures A.14 Modelling of a road containing a reduction of speed restriction under congestedflow traffic conditions. ........................................................................................ 120 A.15 Fundamental diagram: flow-density relationship for different lorry concentration. ............................................................................................................. Values of the weight and roots of the usual Gauss-Hermite quadratures. . II.1 Parameters taken for the hydrostatic pressure simulation................................. 45 II.2 Parameters taken for the Poiseuille flow under gravity simulation. ................. 46 II.3 Parameters taken for the Von Kármán vortex street different regimes. . . III.1 Parameters taken for the anisotropic diffusion reaction for segmentation. . . xxiii

I. 4

 4 .1.3 Truncation order . . . . . . . . . . . . . . . . . . . 23 I.4.2 Gauss-Hermite quadrature . . . . . . . . . . . . . . . . . . 23 I.4.2.1 Optimization . . . . . . . . . . . . . . . . . . . . 23 I.4.2.2 Usual formulations . . . . . . . . . . . . . . . . . . 24 I.4.3 Discretization in space and time . . . . . . . . . . . . . . . . 27 I.4.4 Algorithm: Streaming and colliding . . . . . . . . . . . . . . . 29 I.5 The numerical developments of the LBM . . . . . . . . . . . . . . . . 29

Figure I. 1 -

 1 Figure I.1 -Graph of the potential energy resulting from the attractive and repellent forces after the Lennard-Jones model.

Figure I. 2 -

 2 Figure I.2 -Sequence of the main concepts composing the LBM.

Chapter

3 -s

 3 Figure I.3 -Scheme of the discretization of the velocity space.

FigureFigure I. 7 -

 7 Figure I.7 -Collision and streaming sub-step illustration.

  Chapter I. Context and State of the Art 3 0

Figure I. 8 -

 8 Figure I.8 -Graphical representation of the entropic LBM.

II. 1 . 1

 11 Function LBM_explicit_time(): Data : Lattice, boundary & initial conditions Result : Density and speed fields EndFunction Algorithme II.1 : Pseudo-code LBM

Figure II. 1 -

 1 Figure II.1 -Organization of the computational tool.

Figure II. 2 -

 2 Figure II.2 -Organization of the warping environment.

Figure II. 3 -

 3 Figure II.3 -Schematic connections between warping environment and the GUI.

  Figure II.4 gives a preview of the GUI.

Figure II. 4 -

 4 Figure II.4 -Preview of the GUI.

  ∇x • (p) = ρg. (II.1)This leads to the famous hydrostatic relation∆p = ρgH (II.2)where H is the height of the column of water considered. The parameters considered for the simulation are given in TableII.1TableII.1 -Parameters taken for the hydrostatic pressure simulation.

Figure

  Figure II.5 -Pressure field along a 2D column of water under gravity.

  Figure II.6 -Norm L 2 of the speed field along a 2D flow of air under gravity.

Figure II. 7 -

 7 Figure II.7 -Graph of the speed in the direction of gravity along the orthogonal direction.

Figure II. 6

 6 Figure II.6 and Figure II.7 give the results obtained by the simulation run from the implemented LBM algorithm after 20 000 time steps. Figure II.6 represents the norm of the speed in the pipe, and the non-variation of the speed along the x-axis is following the theory. Figure II.7 shows the parabolic profile of the flow in the pipe for both theory and simulated values. A very thin difference is visible on Figure II.7 which is confi the good implementation for this test case. The maximum speed resulting is 3.113 38.10 -2 , this means a relative error of 0.54%.

  Figure II.8 -Norm L 2 of the speed field along a 2D flow of air under gravity.

Figure

  Figure II.8 and Figure II.9 show the results given by the simulation run for an inlet profile flow of maximum dimensionless speed of 0.15 and after 10 000 iterations to observe the steady flow.The maximum speed resulting in the outlet is 1.4757.10 -1 , this means a relative error of 1.62%. This relative error is larger than the previous test case. A likely explanation is the fact that the periodic boundary condition was not used in the simulation.

Figure II. 9 -

 9 Figure II.9 -Graph of the speed in the direction of gravity along the orthogonal direction.

  Figure II.10 -ParaView visualization of the streamlines and norm of the speed for Von Kármán vortex street.

Figure III. 1 -

 1 Figure III.1 -Illustrative computation of blood flow through an aortic aneurysm from medical image.

Figure III. 2 -

 2 Figure III.2 -Schematic of the probability of streaming. The upper row is the result of the collision step where the two columns are adjacent cells and the colours are identifying the origins of the f i ; the second row is during the streaming step where a portion of the f i are not propagating; and the lower row is the result after the streaming step.

( a )

 a Figure III.3 -Result of segmentation with a LBM anisotropic diffusion reaction [Wan14].

  ) ) (eq. (I.92))III. 1. 10f i ← Stream(f i ) (eq. (I.93)) Algorithme III.1 : Pseudo-code LB3M for elementary mathematical morphology operations.

Figure III. 4 -

 4 Figure III.4 -Morphological dilation in greyscale using the LBM applied to Lena's image.

  Figure III.5 -Morphological erosion in greyscale using the LBM applied to Lena's image.

Figure III. 6 -

 6 Figure III.6 represents the growth of cells on the bone substrate at three given times. Figure III.6a shows a slight densification in light grey on the border of the black substrate. Conversely, the Figure III.6b shows a more pronounced densification with a dark grey at the border. While Figure III.6c illustrates well the growth in an area with a high density of extracellular matrix (dark grey) wider than for Figure III.6b, and a slight cellular growth on the edge of the latter (light grey). This application remains illustrative because

Figure III. 7 .Figure III. 7 -

 77 Figure III.7 -Schematic representing the order to use the different LBM layers

Figure III. 8 -

 8 Figure III.8 -Result of the thrombus segmentation with a LBM

Figure III. 10 -

 10 Figure III.10 -Schematic of the thrombosis reaction system

Figure III. 12 -

 12 Figure III.12 -Representation of the onionskin structure of the giant aneurysm

  Chapter IV. Generalized Out-Equilibrium Distribution 70 IV.6 Conclusion of the chapter ........................................................................................ 83

FigureFigure IV. 1 -

 1 Figure IV.1 -Organization of the main concepts composing the LBM

  17

  ) where c = ξv is the centred microscopic speed of the particles and qα = Q αββ . Therefore, it is possible to define through a change of variable the general out of equilibrium distribution function which yields the undefined moments in order to solve the π Chapter IV. Generalized Out-Equilibrium Distribution 74 l s s -desired equations of motion. Thus, with the variables:

-

  Figure IV.2 shows a schematic illustration of the test case. At the initial state, the whole grid is filled with the spring with computed analytical x Density ρ(x)

  Figure IV.2 -Schematic of a 1D spring compressed test case and the image associated.

Figure IV. 3 -

 3 Figure IV.3 -Displacement and density in function of space at the compression steady stage.

Figure IV. 4 -

 4 Figure IV.4 -Displacement and density in function of space at the decompression steady stage.

Figure

  Figure IV.4 shows the resulting ρ and u after another 25000 time step which can be considered as a steady state for the same reasons as the compression stage.So, the results of the 1D test case are in agreement with the theory for compression and decompression of spring in static situations.

Figure IV. 5 -Figure IV. 6 -

 56 Figure IV.5 -Schematic of 2D pure shear test case and the image associated.

Figure V. 1 -

 1 Figure V.1 -Sequence of main concepts composing the LBM

Figure V. 2 -Figure V. 3 -

 23 Figure V.2 -Displacement in lattice units along the 1D spring at different time step.

FigureFigure V. 4 -

 4 Figure V.4 -Schematic differences between surface and volume forces in LBM, where the left boxes are the density post collision and right boxes are those post streaming.

  Function collide():

Figure V. 5 -Figure V. 6 -Figure V. 7 -

 567 Figure V.5 -Displacement and divergence of the stress in lattice units along the 1D spring at 15 time step.

Figure V. 9 -y

 9 Figure V.8 -A shear displacement at 2300 time step.

  Schematic of the local compression of a 2D (b) Associate image of the local compression simshell.ulation.

Figure V. 11 -

 11 Figure V.11 -Schematic of the local compression simulation and its associate image.

Figure V. 12 -Figure V. 13 -

 1213 Figure V.12 -Norm of the displacement and the speed in a 2D-shell locally compressed after 75 time step.

  Figure V.14, Fi gure V.15 and Figure V.16 shows the displacement of the main patterns appearing. The Figure V.14, looks like a mode 3 in the transverse direction, while the Figure V.15 seems to be the fundamental mode of the plate. The Figure V.16 could be a mode 3 in the transverse direction combined with a mode 2 in the normal direction.

  Figure V.14 -A displacement at time step.

  relaxation term and the interaction term.

  Figure A.2).

Figure A. 2 -

 2 Figure A.2 -Schematic of an asymmetric D1Q6 network.

  (xeiδt) =ρ can be understood as

  m

  eiδt) + fi-1(xeiδt) if ρtest -m \ fj (xejδt) n (xe δ )

  10) can be � to the global density to * -respect of conservation equations leads to a new formulation of the equilibrium density function:(x) \ -1ρc(x) \ for i = 0 c,i = i exp e 2 ρc(x) -i � 1ρc(x) ρc(x) (A.27) 1 + vm,c , e 2 exp � ( e i (x) \ -∀ i ∈ [[1, vc,m ] i=1 1ρc(x) ρ (x) = � 0 ∀ i ∈ [vc,m + 1, m ] vm,cis the maximum speed of the class c vehicles, and ρc is the reachable forward occupation rate for the class c.

Bv

  (fi(x)) = 0 if ρtestfj,c(xejδt) nl(xej δt) > 1 nl(x) * i+1≤j≤m fi,c(xeiδt) = f (xe δ ) else 0≤c≤Nc nl(xej δt) i-1,c (xeiδt) = fi,c(xeiδt) + fi-1,c(xeiδt) if ρtestn (x e δ ) ρtest(x) = \ \ fi,c(xeiδt) l i t . (A.30) i=0 c=0 nl(x)

Figure A. 3 -

 3 Figure A.3 -Schematic of a ring road.

Figure A. 5 -

 5 Figure A.4 -Fundamental diagram of some macroscopic models.

FigureFigure A. 6 -Figure A. 7 -Figure A. 8 -

 678 Figure A.6 -Schematic of merging roads.

Figure A. 9 -

 9 Figure A.9 -Schematic of a road with number of lane change.

Figure A. 10

 10 Figure A.10 represents density evolution in time and space for a simulation with the entrance density of 0.10. This situation leads to a density increase after the reduction of lanes, which remains in the free-fl w domain. The Figure A.11, has for incoming density 0.20, which means an increase of its density will bring it into the congested-fl w domain and might create important interactions. On Figure A.11, the entrance in the congested-fl w domain leads to a slowing-down situation waving backward.

Figure A. 11 -

 11 Figure A.10 -Density for reduction of lanes in free-fl w conditions.

Figure A. 12 -Figure A. 13 -

 1213 Figure A.12 -Schematic of a road with a change of speed limit.

  Figure A.3). The injection point imposes a global density ρp distri buted in two classes, through a coeffi t α. Thus, ρ1(xp) = αρp and ρ2(xp) = (1α)ρp. The Figure A.15 shows the influence of class of slower vehicles (with index 2), with speed limit of 4 cells per unit of time (to simulate the heavy weighted machines), on a faster class (with index 1) having a speed limit of 5 cells per unit of time (to model the personal cars). The density ρp is set to 0.725, in order to run through almost all the traffic domains. The Figure A.15 shows the fundamental diagrams obtained.

Figure A. 15 -

 15 Figure A.15 -Fundamental diagram: flow-density relationship for different lorry concentration.

  Chapman-Enskog expansion . . . . . . . . . . . . . . . . . 18 I.3.10 Study of the first order: the NSF equations . . . . . . . . . . . . 19 I.4 Lattice Boltzmann Method . . . . . . . . . . . . . . . . . . . . . 21 I.4.1 Projection on Hermite basis. . . . . . . . . . . . . . . . . . 21 I.4.1.1 Study of the basis ........................................................................... 21 I.4.1.2 Projection of the passage formula................................................... 22 I.5.1 The origins from the Lattice-Gas Cellular Automata . . . . . . . . 30 I.5.2 The pioneers . . . . . . . . . . . . . . . . . . . . . . . . 30 Segmentation . . . . . . . . . . . . . . . . . . . . 34 I.6.1.4 Morphological image processing . . . . . . . . . . . . . 34 I.6.1.5 Image processing with the lattice Boltzmann method . . . . 34 I.6.2 Image-based diagnostic with macroscopic methods . . . . . . . . . 35 I.7 Emanating questions . . . . . . . . . . . . . . . . . . . . . . . . 35 I.8 Conclusion of the chapter . . . . . . . . . . . . . . . . . . . . . . 36

	Chapter I. Context and State of the Art
	I.3.9
	Chapter
	6

I.5.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . 30 I.5.4 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 I.5.5 Multiphase-Multicomponent . . . . . . . . . . . . . . . . . . 31 I.5.6 Multiple Relaxation Time . . . . . . . . . . . . . . . . . . . 32 I.5.7 Entropic . . . . . . . . . . . . . . . . . . . . . . . . . . 32 I.5.8 LBM for solids . . . . . . . . . . . . . . . . . . . . . . . 33 I.6 Image-based diagnostic . . . . . . . . . . . . . . . . . . . . . . . 33 I.6.1 Some notions about image processing . . . . . . . . . . . . . . 33 I.6.1.1 Denoising . . . . . . . . . . . . . . . . . . . . . . 34 I.6.1.2 Classifi . . . . . . . . . . . . . . . . . . . . 34 I.6.1.3

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 II.2 Specifi and paradigms . . . . . . . . . . . . . . . . . . . . . 38 II.3 Structure and capacities. . . . . . . . . . . . . . . . . . . . . . . 40 II.3.1 Computational code . . . . . . . . . . . . . . . . . . . . . 40 II.3.1.1 Choices: programming language and compilers.............................40 II.3.1.2 Libraries . . . . . . . . . . . . . . . . . . . . . . 40 II.3.1.3 Optimization . . . . . . . . . . . . . . . . . . . . 41 II.3.2 Environment . . . . . . . . . . . . . . . . . . . . . . . . 42 II.3.3 Graphical User Interface. . . . . . . . . . . . . . . . . . . . 43 II.4 Benchmark to validate the initial code .....................................................................44 II.4.1 Hydrostatic pressure . . . . . . . . . . . . . . . . . . . . . 45 II.4.2 Poiseuille flow under gravity . . . . . . . . . . . . . . . . . 46 Conclusion of the chapter .........................................................................................51

	Chapter II. Implementation & validation of a LBM
		Computer programming is an art, because it applies
		accumulated knowledge to the world, because it requires
		skill and ingenuity, and especially because it produces
		objects of beauty.
		Donald Knuth, "Computer Programming as an
		Art". 1974
	Chapter	II
	Implementation & validation of a
		LBM
	Contents of the chapter
	38

Abstract of the chapter . . . . . . . . . . . . . . . . . . . . . . . 38 Résumé du chapitre . . . . . . . . . . . . . . . . . . . . . . . . 38 II.1 II.4.3 Poiseuille flow with Zou-He conditions. . . . . . . . . . . . . 47 II.4.4 Von Kármán vortex street . . . . . . . . . . . . . . . . . . 49 II.5

1. 5 end III. 1. 6 for tStep = 1 to tMax do III. 1. 7

	III. 1. 2	InitiateLattice()
	III. 1. 3	if erosion then
	III. 1. 4	S

f ← ApplyBounceBack(S f ) (eq. (I.94)) III. ρ ← ComputeMacros(f i ) (eq. (III.6) or eq. (III.10)) III. 1. 8 f i (0) ← ComputeFeq(ρ, S f ) (eq. (III.7) or eq. (III.11)) III. 1. 9

Je tiens à remercier mon directeur de thèse Laurent Navarro de m'avoir donné l'opportunité d'effectuer cette thèse, son encadrement, ses conseils précieux et sa patience qui m'ont permis de mener à bien ce travail. Je tiens également à exprimer tout aussi chaleureusement ma gratitude à mon co-directeur Guy Courbebaisse pour ses conseils stratégiques et avisés. Merci à vous deux pour votre amitié et votre temps.En guise de transition entre l'encadrement et le laboratoire, je vais commencer par Stéphane Avril le directeur qui fut également dans l'encadrement de cette thèse aux premières heures, que je remercie pour ses aides et conseils. Je voudrais également exprimer ma gratitude à toutes les personnes du laboratoire CIS qui m'ont rendu la vie plus facile pendant cette période, et elles ont été de plus en plus nombreuses avec le temps. En pensant à tous ces moments de convivialité, je suis obligé de penser aux pauses café avec toutes leurs discussions ; mais comment penser à ces pauses sans penser au roi des mots fléchés et gestionnaire logistique du café Thierry. Alors merci à toi Thierry pour ta convivialité et sympathie à peine arrivé et que tu as su entretenir. Ensuite, dans le désordre je voudrais remercier Vincent le geek qui ne s'en cache pas, Raksmey toujours disponible et prêt à jeter un oeil sur le code, puis Benjamin la relève pour le quota lyonnais. Merci également aux permanents de l'équipe StBio, à Jérôme l'ancien grimpeur toujours prêt à vous indiquer un spot à vos risques et périls, à Pierre B le champion du vélo et des bons conseils, Baptiste P. que l'on dit plus prudent en montagne depuis qu'il est papa, Claire au grand coeur qui n'a pas peur des grosses équations, Woo-Suck dont les questions sont parfois déroutantes. Il me faut aussi exprimer ma gratitude à nos assistantes Françoise et Amélie, sans qui

III.6 Conclusion ............................................................................................................... 67

Maintenant que les fondations sont validées, ce chapitre apporte un nouveau concept: l'idée de coupler la LBM pour la CFD avec le traitement des images dans un seul outil. Aussi, avant de plonger directement dans ce nouveau concept, il sera introduit progressivement en étudiant chacune de ses parties constitutives. Ainsi, en utilisant le constat du lien entre le support des images et le réseau LBM comme point de départ, cela permet d'effectuer des simulations à partir des images. Ensuite, un premier regard sur LBM pour le traitement d'images est donné avec utilisation de diffusion anisotrope. Cette utilisation de la LBM pour le traitement d'images est étendue avec le MM à travers une preuve théorique. Cette extension est pratique puisque la MM peut également être utilisée pour la simulation numérique telle que la croissance biologique. Pour finir, le nouveau concept apparaît naturellement comme le regroupement de tous ces développements sur un seul réseau LBM. Ce concept est ensuite illustré sur un exemple d'anévrisme cérébral.

www.thrombus-vph.eu

Remerciements

Je tiens ici à remercier toutes les personnes qui ont apporté, à leurs manières, leurs soutiens et contributions à cette thèse. J'ai conscience de la chance que j'ai, car ils sont nombreux. Pour autant, ils ne me rendent pas la tâche si aisée, car il ne faut oublier personne. Je vais essayer.

Je suis en premier lieu très reconnaissant aux membres du jury : leur validation est un élément rassurant et confortant mes choix. Merci à Thierry Fournel d'avoir accepté rapidement de participer à ce jury et pour ses questions pertinentes. Je souhaite remercier les rapporteurs Yu Chen et Marianne Clausel d'avoir examiné mes travaux, lu attentivement ce manuscrit et formulé des critiques constructives. Je voudrais bien sûr exprimer ma gratitude au président du jury Claude-Henri Lamarque pour la confi et le soutien qu'il m'apporte depuis plusieurs années, ainsi que ses remarques, questionnements et conseils.

Abstract of the chapter

Now that the foundations are validated, this chapter brings a new concept: the idea of coupling the LBM for CFD with image processing in a single tool. Also, before to dive directly in this new concept, it will be introduced gradually by studying each of its composing parts. So, the notice of the link between image support and LBM network is used as a starting point, which allows to perform simulation from images. Then, a first look at LBM with anisotropic diffusion use for image processing is added. This use of LBM for image processing is extended through theoretical proof with the MM. This extension is convenient since the MM can also be used for numerical simulation such as biological growth. To finish, the new concept appears naturally as the gathering of all these developments on a single LBM network. This concept is then illustrated on a cerebral aneurysm example.

Résumé du chapitre

Conclusion and future works

A

Lattice Boltzmann Method for Heterogeneous Multi-class Traffi Flow

Contents of the chapter

Abstract of the appendix

The traffic modelling often keeps the mesoscopic scale in the theoretical sphere because the integro-differential nature of its equations. In the present work we suggest to use the lattice Boltzmann method to overcome these difficulties. In particular, the method has a strong theoretical foundation. An improved version of the lattice Boltzmann method

Appendix

A.2.2 Continuum Kinetic Multi-Class Traffic Approach

Some authors [HB00; CD15] construct continuum kinetic model for multi-class flows, through numerous inter-class interaction. However, some authors [Cha05; BR03] remind that the use of multi-class distribution is easier.

In this work, an alternative formulation is used. Indeed, guided by the i ntui ti on that every driver interacts in the same manner with all the other vehicles without independentl y of the class of the other vehicles. So, the construction of the equilibrium distribution for each class should be the function of the global density over the class and the local densi ty of its own class to respect the continuity equation.

Let us denote with the subscript c one class of vehicle out of Nc the number of classes. Thus, the mesoscopic multi-class Boltzmann-like equations are

To be exact, densities ρ and ρc considered here should be interpreted as non-dimensional densities or equivalently as occupation rates. Indeed, let us consider the case where there is two classes, the first one is made of vehicles of length L1 and respectively the second of length L2. If the length of the Representative Elementary Volume (REV) is Dx, therefore the maximum number of vehicles from the first class at the position x is denoted ρ1(x) and has the maximum value of Dx/L1. This value is reached when there is no free space between vehicles, which leads to the definition of ρc as the occupation rate. Thus for any class, the relationship between ρ and ρ is

Moreover, the global density can only be expressed in terms of vehicles number equivalent to vehicles of one class. In the previous example, the number of vehicles from the second class equivalent to the ones of the first class is ρ 2/1 = ρ2L2/L1. And so, if ρ tot/1 is the total number of vehicles equivalent to those of the first class, this number is given by

which is more convenient to compute with rate occupations (see eq. (A.10))).

The major drawback with kinetic models (single or multi-class) is the integro-differential nature of the equations. This fact, combined with the lack of physical properties, is responsible for absence of analytical solutions. This is why numerical methods are for now necessary.

L

A.3 Lattice Boltzmann Method

A.3.1 Lattice Boltzmann Model for Heterogeneous Traffic Flows

The lattice Boltzmann Method (LBM) is a discretisation of the continuous Boltzmann equation. This discretisation is performed on all three space, time and velocity space. The time discretisation gives an explicit schema. The phase space discretisation is often characterised by the appellation DnQm where n specifi the physical dimension of the problem and m is the number of points used to condense the phase space. The denomination lattice is linked to the regular spatial discretisation and the DnQm schema used to connect the spatial points; to coincide with the velocities [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF].

Commonly in the LBM approach, the assumption of a collision-interaction operator composed only of a relaxation term is made. This assumption is named BGK in reference to Bhatnagar, Gross, Krook [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. Thus, under the BGK approximation the lattice Boltzmann equation (LBE) reads [Wol00; Cho+02]

where fi is the distribution evaluated in ei the discrete velocities (associate to the phase space), δt is the time step, therefore the space or lattice step is given by δx/δt = e.

For the numerical resolution the LBE is divided in two steps. The first step is the collision-interaction step defined by the equation

This is usually followed by the streaming step defined by: fi(x + eiδt, t + δt) = fi(x, t + δt).

(A.15)

Nevertheless, the last equation must be adapted to take into account possible lane number changes (see Figure A.1). Therefore, to face these possibilities, it is suggested to turn the eq. (A.15) into:

with n l is the number of lanes for a given spatial point.

The macroscopic variables are recovered by using the classic following summations: 122 tain traffic situation. An example could be to incorporate the psychological behaviour of drivers. Other aspects such as the sinuosity of the road or the traffic pressure in multi-class flows could also be interesting like stochastic aspects. Investigations on the equilibrium distribution functions and their justifi would be a major work.