Bernard Rousseau

M Guillaume

Mme Irina

Martin Graur

Ce travail de thèse a été réalisée au sein du laboratoire Modélisation et Simulation Multi Echelle (UMR 8208 CNRS) de l'université Paris-Est et il n'aurait probablement jamais été mené à terme sans le soutien d'un grand nombre de personnes que je tiens vivement et très sincèrement à remercier.

Je commencerai d'abord par remercier Monsieur Guy Lauriat qui a dirigé cette thèse dans la continuité de mon stage de master. La pleine conance qu'il m'a accordée dès le début m'a permis de progresser régulièrement. Je voudrais aussi le remercier pour le temps qu'il m'a consacré et la patience avec laquelle il a accompagné mon travail.

Je remercie particulièrement Monsieur Quy Dong To, mon co-encadrant de thèse, pour sa grande contribution à ce travail, son soutien, sa bonne humeur ainsi que l'amitié qu'il m'a témoignée tout au long de ces trois années.

Je voudrais également remercier Madame Céline Léonard pour sa collaboration sur la partie de la

Résumé

Dans les systèmes micro ou nano-uidiques, le glissement de vitesse à l'interface uideparois est un phénomène important. Pour modéliser cet eet à l'échelle macroscopique, les conditions aux limites de Navier ont été introduites, avec la longueur de glissement comme paramètre. Lorsque le uide est un gaz, cette longueur est liée au coecient d'accommodation tangentiel (TMAC) et au libre parcours moyen, selon le modèle de Maxwell. Le but de ce travail est de traiter systématiquement ce modèle par une approche multi-échelle et de l'étendre en incorporant la morphologie et l'anisotropie de la surface.

La thèse est composée de cinq chapitres.

Après l'Introduction, les notions de base de la théorie cinétique des gaz, l'équation de Boltzmann et les solutions associées (Navier-Stokes-Fourier, Burnett, Grad, Direct Simulation Monte Carlo...) sont rappelées dans le Chapitre 1. Les modèles d'interaction gaz-paroi ainsi que les modèles de glissement introduits dans le cadre de la mécanique des uides sont aussi rappelés. Le chapitre se termine par la description de la méthode de calcul par dynamique moléculaire (MD) utilisée dans ce travail.

Le Chapitre 2 est dédié au développement d'une technique simple an de simuler les écoulements induits par la pression. Le principe est de se baser sur les formules atomistiques du tenseur des contraintes (Irving Kirkwood, Méthode de Plan, Contraintes Virielles) et de modier les conditions périodiques, tout en maintenant la diérence entre l'énergie cinétique des atomes à l'entrée et à sortie du domaine de calcul. Plusieurs types de conduite sont étudiés avec cette technique. Les résultats (température, vitesses...) sont discutés et comparés.

Le Chapitre 3 concerne l'étude du potentiel d'interaction gaz-paroi par la méthode abinitio. Le code CRYSTAL 09 est utilisé pour obtenir le potentiel entre un atome d'argon (Ar) et une surface de platine (Pt) <111> en fonction de la distance. Ensuite, le potentiel atome/surface est décomposé en potentiel binaire et approché par une fonction analytique. Cette fonction est ensuite implémentée dans un code MD pour simuler les collisions gazparoi et déterminer le coecient TMAC.

Dans le Chapitre 4, l'eet de morphologie est étudié. Le potentiel multi-corps Quantum Sutton Chen (QSC) est utilisé pour le solide Pt <100> et un des potentiels binaires étudié dans le chapitre précédent pour le couple Ar-Pt. Le potentiel QSC est nécessaire pour reproduire l'eet de surface qui aecte le résultat nal. Diérentes surfaces sont traitées : surface lisse, surface nanostructurée, surface aléatoire obtenue par déposition de vapeur (CVD). Le coecient TMAC est déterminé de façon généralisée, c.à.d en fonction de l'angle du ux d'atomes incidents sur la surface. Les anisotropies de surface et le noyau de collision sont également examinés.

Dans le Chapitre 5, on propose un modèle de glissement anisotrope pour les uides en fonction du tenseur d'accommodation. Le modèle est obtenu par les calculs analytiques approchés développés dans le cadre de la théorie cinétique. On a ainsi généralisé l'équation de Maxwell en montrant que le tenseur de longueur de glissement est directement lié au tenseur d'accommodation. Le modèle est en bon accord avec les résultats de la méthode MD. Concernant la simulation MD, on a développé une technique pour reproduire l'anisotropie du tenseur d'accommodation.

Le mémoire de thèse se termine par une synthèse des résultats obtenus. Des perspectives pour de futures études sont proposées.

)>IJH=?J

In most applications concerning a uid owing over a solid surface, the no-slip velocity condition was widely used because it is simple and produces the results in agreement with experiments. However, this dynamical boundary condition is not appropriate when the ow under consideration is at a micro or nano length scale. In order to model this eect at the macroscopic scale, the Navier boundary conditions have been introduced, with the slip length as a parameter. When the uid is a gas, this length is related to the tangential momentum accommodation coecient (TMAC) and the mean free path, according to the Maxwell model. The aim of this work is to systematically address this model using a multi-scale approach and to extend it by incorporating both the morphology and the anisotropy of a surface.

The thesis consists of ve chapters.

In Chapter 1, the basics of the kinetic theory of gases, the Boltzmann equation and related solutions (Navier-Stokes-Fourier, Burnett, Grad, Direct Simulation Monte Carlo ...) are briey presented. The models of gas-wall interaction and slip models introduced in the uid mechanics are also recalled. The chapter ends with a description of the computational method used for the molecular dynamics simulations performed in this work.

Chapter 2 is dedicated to the development of a simple technique to simulate the pressure driven ows. The principle is to rely on the atomistic formulas of the stress tensor (Irving Kirkwood, Method of Plane, Virial Stress) and to modify the periodic conditions by maintaining the dierence between the kinetic energy of the ingoing and outgoing particles of the simulation domain. Several types of channels are studied with this technique. The results (temperature, velocity ...) are discussed and compared.

Chapter 3 deals with the study of the gas-wall interaction potential by the ab-initio method. The code CRYSTAL 09 is used to obtain the potential between an atom of argon (Ar) and a surface of platinum Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] as a function of distance. Then the gas-wall potential is decomposed into binary potential and approached by an analytic function. This function is then implemented in a MD code to simulate the gas-wall collisions and determine the TMAC coecient. In Chapter 4, the eect of morphology is studied. The multi-body Quantum Sutton Chen (QSC) potential is used for Pt [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF] solid and the binary potential proposed in the previous chapter for the Ar-Pt couple is employed. The QSC potential is needed to reproduce the surface eects that aect the nal results. Dierent surfaces are treated : smooth, nanostructured surface and, random surface obtained by Chemical vapor deposition (CVD). The TMAC is determined using a generalized approach, i.e. depending on the angle of incident ux of gas atoms on the surface. The surface anisotropy and the scattering kernel are also examined.

In Chapter 5, we propose a model of anisotropic slip for uids based on accommodation tensor. The model is obtained by the analytical approximate calculations developed in the framework of the kinetic theory. We thus generalize Maxwell's equation by showing that the slip length tensor is directly related to the accommodation tensor. The model is in good agreement with the MD results. Thanks to our MD simulations, we develop a suitable technique for reproducing the anisotropy of the accommodation tensor .

The thesis ends with a conclusion section in which we suggest some perspectives for a continuation of this work.

Introduction

Over the past decades, micro/nanouidic systems have developed rapidly as the result of continuous progresses in micro/nanofabrication technologies. Such systems contribute to a technological revolution in many areas such as the cooling of electronic components, micro electro mechanical systems (MEMS), biotechnology and medicines. One of the most successful commercial application of microuidics is the inkjet printhead which consist of an ink reservoir fabricated from silicon, a heating element to drive the uid, and a nozzle (Fig. 1). Droplets created by this system is generally about 50 µm. According to Tabeling, tens of millions of inkjet printers use MEMS and billions of documents are written and read thanks to microuidics [START_REF] Tabeling | Introduction to microuidics[END_REF].

Figure 1: An inkjet printhead (left) and the visualization of droplets of ink projected onto a target (right) [START_REF] Tabeling | Introduction to microuidics[END_REF], [START_REF] Ho | ASME Int. Mechanical Engineering Congress and Exposition[END_REF].

Another examples of microuidic systems which can be found in computer science, mechanical and chemical engineering are chemostats, pumps, conducts, mixers, lters, and valves (Fig. 2). Besides, the actual technologies permit also the construction of microdevices containing multiple components with dierent functionalities. A single integrated chip could perform many works concerning the biological or chemical processing from the sampling, pre-processing, and measurement in an assay. This leads to save time and reduce the reagent consumption. Nowadays, lab-on-a-chip devices are used for abundant elds including analysers of air and water quality, diagnostics of illnesses, and devices that replace many functions of components of human body. The economic possibilities of this eld have been estimated at tens of billions of dollars per year [START_REF] Tabeling | Introduction to microuidics[END_REF]. As many investigators pointed out, physics that apply to nanotechnology scale is dierent from that on a human scale. Indeed, in the macroscopic scale, the surface/volume ratio is on the order of unity while in nanoscale, it rises to nearly a million. Such large ratios allow the surface eects dominate the physics and quantum eects begin to show and thus become important. The understanding of these physical eects at such small scales is therefore a vital interest in the manufacture and use of micro/nano electrodynamic systems.

At the same time, improvements in the numerical modeling are also needed to analyze the reliability of the systems. Depending on the nature of the problem, there are dierent suitable computational methods, such as those based on the Navier-Stokes and energy equations, the Burnett or Super Burnett models [START_REF] Wu | A slip model for rareed gas ows at arbitrary knudsen number[END_REF], or the molecular simulations. In the great family of molecular simulation methods, this is the Monte-Carlo method, which rst appeared in 1953 with the work of Metropolis et al. [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. This method is purely statistical and is well suited to the study of thermodynamic properties. The method is to successively change the movement of molecules randomly following a conditional probability by the thermodynamic state of the considered system. However, the classic Monte-Carlo method does not allow to study the dynamic aspects of systems. First proposed by Alder et Wainwrigth [START_REF] Alder | Phase transition for a hard sphere system[END_REF] in 1957 , the molecular dynamics (MD) simulation can overcome this obstacle. In general, there are two main types of MD simulation, which can be distinguished according to the way of interaction forces calculations. In classical MD simulations, the interaction forces are derived from the xed empirical potentials, while in ab initio MD simulations, the quantum nature of the chemical bond is explicitly taken into account and the interaction forces are determined from the quantum chemistry. By integrating the equations of motion at time t, the molecular dynamics allows to evolve over time a system of particles dened by their positions, velocities and accelerations. As the user desires, the molecular properties can be sampled, averaged, and used to compute ow quantities, such as velocity, density, and viscosity. For micro/nanoows, the MD method is one of the most accurate method since realistic interactions between particles or between particles and boundaries may be accounted for.

In most applications concerning a uid owing over a solid surface, the no-slip velocity condition was widely used because it is simple and produces the results in agreement with experiments. However, this dynamical boundary condition is not appropriate when the ow under consideration is at a micro or nano length scale. In order to model this eect at the macroscopic scale, the Navier boundary conditions have been introduced, with the slip length as a parameter. When the uid is a gas, this length is related to the tangential momentum accommodation coecient (TMAC) and the mean free path, according to the Maxwell model. The aim of this work is to systematically address this model using a multi-scale approach and to extend it by incorporating both the morphology and the anisotropy of a surface.

The thesis consists of ve chapters.

In Chapter 1, the basics of the kinetic theory of gases, the Boltzmann equation and related solutions (Navier-Stokes-Fourier, Burnett, Grad, Direct Simulation Monte Carlo ...) are briey presented. The models of gas-wall interaction and slip models introduced in the uid mechanics are also recalled. The chapter ends with a description of the computational method used for the molecular dynamics simulations performed in this work.

Chapter 2 is dedicated to the development of a simple technique to simulate the pressure driven ows. The principle is to rely on the atomistic formulas of the stress tensor (Irving Kirkwood, Method of Plane, Virial Stress) and to modify the periodic conditions by maintaining the dierence between the kinetic energy of the ingoing and outgoing particles of the simulation domain. Several types of channels are studied with this technique. The results (temperature, velocity ...) are discussed and compared.

Chapter 3 deals with the study of the gas-wall interaction potential by the ab-initio method.

The code CRYSTAL 09 is used to obtain the potential between an atom of argon (Ar) and a surface of platinum (Pt)⟨111⟩ as a function of distance. Then the gas-wall potential is decomposed into binary potential and approached by an analytic function. This function is then implemented in a MD code to simulate the gas-wall collisions and determine the TMAC coecient.

In Chapter 4, the eect of morphology is studied. The multi-body Quantum Sutton Chen (QSC) potential is used for Pt ⟨100⟩ solid and the binary potential proposed in the previous chapter for the Ar-Pt couple is employed. The QSC potential is needed to reproduce the surface eects that aect the nal results. Dierent surfaces are treated : smooth, nanostructured surface and, random surface obtained by Chemical vapor deposition (CVD). The TMAC is determined using a generalized approach, i.e. depending on the angle of incident ux of gas atoms on the surface. The surface anisotropy and the scattering kernel are also examined.

In Chapter 5, we propose a model of anisotropic slip for uids based on accommodation tensor. The model is obtained by the analytical approximate calculations developed in the framework of the kinetic theory. We thus generalize Maxwell's equation by showing that the slip length tensor is directly related to the accommodation tensor. The model is in good agreement with the MD results. Thanks to our MD simulations, we develop a suitable technique for reproducing the anisotropy of the accommodation tensor.

The thesis ends with a conclusion section in which we present some perspectives for continuations of this work.

+D=FJAH

Theories and computation methods in microuidics A uid is an ensemble of a large number of molecules that interact and move in space. The macroscopic behavior of the uid can then be completely determined from the knowledge of the particle positions and velocities. However, such a description is not feasible when the number of molecules achieve the order of the continuum limit. A more cost eective way to model uids, especially suitable for gases, is to use a single statistical distribution function [START_REF] Liou | Microuid Mechanics[END_REF].

Consider at time instant t and at location x(x 1 , x 2 , x 3 ) in the physical space, there are N uid molecules occupying a dierential volume dx = dx 1 dx 2 dx 3 . In addition to time and location, these molecules can also be classied by their velocity v(v 1 , v 2 , v 3 ). As a result, the probability density function f (x, v, t) for the velocity can be dened by

dN = f dv, (1.1)
where dN is the number of particle in the sample with velocity in the range of v and v + dv. The term dv = dv 1 dv 2 dv 3 denotes a dierential volume in the velocity space (see Fig. 1.1).

Knowing the distribution function f , the macroscopic properties Q (stresses, temperature, Here, m and ρ = nm are respectively the molecular mass and the mass density

Q Q Stream velocity u v Stress tensor σ -ρ(v -u) ⊗ (v -u) Heat ux q ρ 2 |v -u| 2 (v -u)
Figure 1.2: Phase space element mean velocity, etc...) can be determined by means of average of the molecular quantities Q (velocity, momentum, energy, etc...)

Q = 1 n ∫ Qf dv, n = ∫ f dv, (1.2)
where n is the number density. Generally, these macroscopic quantities are moments of the distribution function f , as listed in Table 1.1.

The Boltzmann equation

The evolution of f with time in physical and velocity spaces can be studied on the basis of the conservative equations. For an element in the combined phase space dvdx shown in Fig. 1.2, the local rate of change of the number particles at time t is

∂f ∂t dxdv, (1.3)
The change in the number of particles within dxdv can be caused by inows of particles through each side of the phase space element. First, the net inow of particles through the physical space element dx is due to the convection of particles across the surface S r of dx.

Considering class v particles, the convection inow reads

- ∫ Sr (f v) • (e r dS r )dv, (1.4)
with e r being the unit normal vector of surface element dS r . Using the Gauss theorem, the surface integral can be rewritten as the volume integral :

∫ dx ∇ • (f v)dxdv .
(1.5)

Since f and v are constant within dx and only particles of class v are considered, the equation becomes

-v • ∂f ∂x dxdv.
(1.6)

Secondly, the ux of particle across sides of dv can be caused by external force per unit mass F. Similar to particle uxes across side of dx, the net inow through the velocity space dv is written as

-F • ∂f ∂v dxdv.
(1.7)

In addition to the above processes, the local number of particles can be changed by intermolecular collisions in the velocity space. This eect is expressed via the integral operator

{ ∂f ∂t } coll dxdv.
(

The collision term can be evaluated by considering the collision between two particles of dierent class, say v and v 1 and their post collision velocities v ′ and v ′ 1 . Skipping the proof which is presented in many previous works [START_REF] Struchtrup | Macroscopic Transport Equations for Rareed Gas Flows: Approximation Methods in Kinetic Theory[END_REF][START_REF] Liou | Microuid Mechanics[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF], we admit the following nal expression

{ ∂f ∂t } coll dxdv = (∫ ∞ -∞ ∫ 4π 0 v r I(v r , Ω)(f ′ f ′ 1 -f f 1 )dΩdv 1 ) dxdv, (1.9)
where f 1 , f ′ and f ′ 1 is the velocity distribution function at v 1 , v ′ and v ′ 1 class and v r is the relative velocities between the particles before and after collision,

v r = |v -v 1 | = |v ′ -v ′ 1 |.
(1.10)

The term I(v r , Ω) is the dierential cross section of the collision, in which the relative velocity of the colliding particles falls within the solid angle dΩ.

Finally, the equation for the change of particle number of class v becomes

∂f ∂t + v • ∂f ∂x + F • ∂f ∂v = ∫ ∞ -∞ ∫ 4π 0 v r I(v r , Ω)(f ′ f ′ 1 -f f 1 )dΩdv 1 .
(1.11) which is the Boltzmann equation for a simple dilute gas.

Approximate solutions

For the gas in equilibrium, the velocity class must be constant with time and position. As a result, the collision term in the Boltzmann equation must vanish

{ ∂f ∂t } coll = 0.
(1.12) Navier-Stokes-Fourier Kn 2 Burnett That condition implies that

f ′ f ′ 1 -f f 1 = 0 or ln f ′ + ln f ′ 1 = ln f + ln f 1 (1.13)
or ln f is a collision invariant. Thus, it is natural to express ln f as linear combination of linear momentum and kinetic energy since they are both collision invariant. Finally, we can derive the famous Maxwellian distribution f 0 .14) with k B being the Boltzmann constant. For nonequilibrium gases, the Boltzmann equation is mathematically dicult to solve due to the integral form of the collision term. In the literature, approximate solutions have been proposed based on the simplied treatment of the latter [911].

f 0 (v) = n ( m 2πk B T ) 3/2 e -m|v-u| 2 2k B T . ( 1
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] replaced the collision operator by a dierence between the sought distribution function f and a local equilibrium distribution f 0 { ∂f ∂t

} coll = 1 τ (f 0 -f ) , (1.15)
where τ is the characteristic time of the evolution from the particle distribution to the equilibrium distribution. The basic idea of the BGK model is to see the collision process as a relaxation phenomenon and the main role of the collision term is to relax the distribution function into an equilibrium distribution. It is also assumed that the actual distribution function is not very far from the equilibrium distribution.

Using the Knudsen number, Chapman and Enskog (CE) [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] employed an asymptotic expansion of

f f = f 0 (1 + Knψ 1 + Kn 2 ψ 2 + ...) , (1.16)
where ψ 1 , ψ 2 , ... are functions of gas density, temperature and macroscopic velocity. By truncating the series at order k of Kn, we obtain the k-order constitutive equations from the Boltzmann equation as shown in Table 1.2.

An alternative approach is proposed by Grad [START_REF] Grad | Principles of the kinetic theory of gases[END_REF], the Grad13 and Grad26 equations, taking higher order moments as variables in addition to the usual quantities temperature, density, etc.... Although these theories are very promising, the development is still underwork, dealing with issues such as boundary conditions, solution existence, etc... At nite Knudsen number, the Boltzmann equation can be also solved numerically with Direct Simulation Monte Carlo (DSMC) method. Representative particles are used in a DSMC simulation instead of real molecules and the collisions between these particles are treated randomly. To date, the DSMC method is considered as the most accurate and ecient numerical method for rareed gas ow simulation and has been widely used in practice.

Flow regimes and validity of models

A key nondimensional parameter for microows is the Knudsen number, which is dened as the ratio

Kn = λ H , (1.17)
where H is the characteristic length scale of ow and λ is the mean free path estimated by

λ = 1 √ 2nπσ 2 .
(1.18)

The quantity σ in Eq. 1.18 is the diameter of the uid molecules. In the case of Lennard Jones (LJ) uid, this diameter is approximately taken equal to the reference length σ appearing in the LJ interaction potential.

In rareed gas dynamics, the Knudsen number is normally used as a measure of the degree of rarefaction. Rarefaction eects become more important as the Knudsen number increases and thus pressure drop, shear stress, heat ux and corresponding mass ow rate cannot be predicted by using models based on the continuum hypothesis [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF]. A classication of the dierent ow regimes and the corresponding governing equations are shown in Fig. 1.3.

• Kn < 10 -3 : the ow is considered as a continuum and the transport phenomena can be modeled by Navier-Stokes and energy conservation equations.

• 10 -3 < Kn < 10 -1 : it is considered a slip ow. The Navier-Stoke and energy equations remain applicable but it is necessary to take into account the boundary conditions with slip velocity and temperature jump in the walls.

• 10 -1 < Kn < 10 : transition ow.

• Kn > 10 : free-molecular ow.

For gas ow with high Knudsen number (> 10 -1 ), the number of particle in a volume unit of gas decreases and the number of intermolecular collisions is probably not sucient to form an equilibrium state. The behavior of individual particle become more important for the macroscopic behavior of gas ow as the Knudsen number increases. The continuum approximation used in the ows with small Knudsen number is then invalid and it is necessary to apply microscopic approaches like the Boltzmann equation, that recognize the molecular property of gases.

Gas-surface interface mechanisms

At the micro/nanoscale, the gas-surface interactions can have signicant impacts on the overall ow behavior. In the near surface region, the gas cannot be considered as in equilibrium and the energy and momentum transfer between the gas particle and the surface should be examined. This area is called the Knudsen layer : a thermodynamically nonequilibrium region extending a distance of a few mean free path from the surface. For ows with small Knudsen number, the eect of Knudsen layer is not signicant, the continuum approximation can be used by extrapolating the bulk gas ow towards the surface. At high Knudsen number, the collision frequency may be insucient to reach an equilibrium state and the distribution function of gas will deviate from the Maxwellian distribution in the Knudsen layer.

The gas-surface interaction can be accessed via some averaged parameters such as the momentum and energy accommodation coecients. While the tangential and normal momentum accommodation coecients are used to characterize the momentum exchanges parallel and perpendicular to the surface, the thermal accommodation coecient is used for the energy exchange.

The Tangential Momentum Accommodation Coecient σ t (abbreviated as TMAC) is dened as [START_REF] Cao | Molecular momentum transport at uid-solid interfaces in MEMS/NEMS: A Review[END_REF] σ

t = v it -v rt v it -u w , (1.19)
where v it and v rt are respectively the tangential velocity of incident and reected particles and u w is the wall velocity.

When σ t = 0, the reection is purely specular. As a result, the tangential velocity component is unchanged, but the normal velocity component is reversed (see Fig. 1.4). This corresponds to an extreme case where the wall is ideally smooth. In the other extreme, the wall is ideally rough and the reection is purely diusive, i.e σ t = 1. The reection occurs at a random angle and the average tangential velocity of reected particles is zero. In Maxwell's work, a mixed reection mechanism is proposed: a portion of the particles reect specularly and the remaining reect diusively.

Since it is impossible to predict theoretically the trajectories or distribution of the incident and reected particles, σ t is usually determined from experiments or numerical simulations As shown in the later sections, the coecient depends on many factors such as gas-surface interaction potentials, the morphology of surface and temperature.

Unlike the TMAC, the Normal Momentum Accommodation Coecient (NMAC) σ n has various forms of denitions as follows [START_REF] Cao | Molecular momentum transport at uid-solid interfaces in MEMS/NEMS: A Review[END_REF],

σ n = |v in | + |v rn | |v in | + |v wn | , (1.20)
σ n = |v in | -|v rn | |v in | -|v wn | , (1.21)
σ n = |v in | -|v rn | |v in | , (1.22)
σ n = |v in | + |v rn | |v in | , (1.23)
in which v in and v rn are the normal velocity of incident and reected particles, respectively, and v wn is given by :

v wn = √ πk B T w 2m , (1.24)
where T w is the wall temperature.

In Equation (1.21), σ n is singular when |v in | = |v wn |, while in the other equations (1.20), (1.22) and (1.23), a singular σ n can be avoided. However, the expressions of (1.22) and (1.23) only show the relative normal momentum of incident and reected particles but cannot represent the accommodation with the wall [START_REF] Cao | Molecular momentum transport at uid-solid interfaces in MEMS/NEMS: A Review[END_REF].

The energy transfer between the monoatomic gas particles and the surface can be characterized via the Energy Accommodation Coecient (EAC) σ T dened by

σ T = v 2 i -v 2 r v 2 i -v 2 w , (1.25)
At rst glance, the coecients σ t , σ n , σ T appear to be identical for all collisions. However, with all the complexity of the wall structure at the atomic scale, σ t , σ n , σ T should be understood as eective quantities computed from a given uxes of momentum and energy. One of the earliest uid-solid boundary conditions was proposed by Navier [START_REF] Navier | Mémoire sur les lois du mouvement des uides[END_REF], postulating that the tangential uid velocity relative to the wall u slip is proportional to the ow shear rate

u slip = u f luid -u wall = L s ( ∂u ∂n ) s , (1.26)
where L s is the slip length, u wall the wall velocity and ∂ ∂n the normal derivative at the wall surface. The slip length L s is shown in Fig. 1.5 with three cases : no-slip, partial slip and perfect slip. For a pure shear ow, L s can be interpreted as the ctitious distance below the real surface where the no-slip boundary condition is satised.

According to Maxwell [START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF], the slip length can be estimated by a simple relation

L s = λ 2 -σ t σ t .
(1.27)

The slip equation is usually written in dimensionless form. When combined with thermal transpiration eect [START_REF] Smoluchowski | Über wärmeleitung in verdünnten gasen[END_REF], it becomes

U slip = 2 -σ t σ t Kn ( ∂U ∂n ) s + 3 2π C p -C v C p Kn 2 Re Ec ∂T ∂s , (1.28)
In (1.28), ∂ ∂s implies tangential derivative at the wall surface, C p and C v the heat specic at constant pressure and volume, Re and Ec the Reynolds number and Eckert number.

The slip velocities in second order have also been proposed to extend the validity range of continuum equations for Knudsen numbers corresponding to the start of the transition regime. For the isothermal ows, they could be written in the general form as :

U slip = C 1 Kn ( ∂U ∂n ) s -C 2 Kn 2 ( ∂ 2 U ∂n 2 ) s , (1.29)
where C 1 and C 2 are the rst and second order slip coecients, respectively. For Maxwell's model in Equation 1.28, C 1 = 2-σt σt and C 2 = 0.

The slip coecients C 1 and C 2 available from the literature are listed in Tables 1.3 and 1.4 [START_REF] Zhang | A review on slip models for gas microows[END_REF][START_REF] Cao | Molecular momentum transport at uid-solid interfaces in MEMS/NEMS: A Review[END_REF]. Maxwell(1879) [START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF] First-order model 1 0 Schamberg (1947) [START_REF] Schamberg | The fundamental dierential equations and the boundary conditions for high speed slip ow, and their application to several specic problems[END_REF] Second-order model 1

5π/12
Albertoni et al. (1963) [19] BGK model 1.1455 0 [START_REF] Deissler | An analysis of second-order slip ow and temperature-jump boundary conditions for rareed gases[END_REF] [START_REF] Deissler | An analysis of second-order slip ow and temperature-jump boundary conditions for rareed gases[END_REF] Second-order model 1

9/8
Chapman and Cowling (1970) [START_REF] Chapman | The mathematical theory of nonuniform gases[END_REF] Linearized Boltzmann model 1 0.5 [START_REF] Loyalka | Approximate method in the kinetic theory[END_REF] [START_REF] Loyalka | Approximate method in the kinetic theory[END_REF] BGK model 0.7252 0 Hsia and Domoto (1981) [START_REF] Hsia | An experimental investigation of molecular rarefaction eects in gas lubricated bearings at ultra-low clearances[END_REF] Second-order model 1 0.5 Loyalka and Hickey (1989) [START_REF] Loyalka | Plane poiseuille ow: Near continuum results for a rigid sphere gas[END_REF] BGK based kinetic model 1.1019 0.0449 [START_REF] Mitsuya | Modied Reynolds Equation for Ultra-Thin Film Gas Lubrication Using 1.5-Order Slip-Flow Model and Considering Surface Accommodation Coecient[END_REF] [START_REF] Mitsuya | Modied Reynolds Equation for Ultra-Thin Film Gas Lubrication Using 1.5-Order Slip-Flow Model and Considering Surface Accommodation Coecient[END_REF] 1.5-order model 1

2/9
Pan et al. (2000) [START_REF] Pan | A modied direct simulation monte carlo method for low-speed microows[END_REF] First-order model 1.125 0 Karniadakis and Beskok (2002) [START_REF] Karniadakis | Micro ows-fundamentals and simulation[END_REF] Second-order model 1 -0.5 [START_REF] Wu | New rst and second order slip models for the compressible reynolds equation[END_REF] [START_REF] Wu | New rst and second order slip models for the compressible reynolds equation[END_REF] New second-order model Using a similar approach made for the slip velocity, von Smoluchowski also proposed a boundary condition for the temperature that can be written in the dimensionless form as,

T jump = T f luid -T wall = 2 -σ T σ T [ 2C p C p + C v ] Kn P r ( ∂T ∂n ) s , (1.30)
where σ T is the thermal accommodation coecient and P r is the Prandtl number given by :

P r = C p η κ , (1.31)
in which η is the dynamical viscosity and κ is the thermal conductivity. 2-σt σt

9 4π P r(γ ′ -1) γ ′ Wu (2008) [35], f = min[1/Kn, 1] 2 3
[

3-σtf 3 σt -3 2 (1-f 2 ) Kn ] 1 4 [ f 4 + 2 Kn 2 (1 -f 2 ) ] 1.

Molecular Dynamics method

Molecular Dynamics was rst introduced by Alder and Wainwright in the late 1950's to study the interactions of hard spheres [START_REF] Alder | Phase transition for a hard sphere system[END_REF]. This method consists in studying the trajectory of particles, simulating the spatial evolution of system over the time. In the simulation, each particle is considered as a point mass whose movement is determined by all the forces exerted on it by other particles. The laws of Newtonian classical mechanics applied to a particle i can be written as :

v i (t) = ∂r i ∂t , (1.32)
m i ∂v i ∂t = F i = - ∂V pot ∂r i
, (1.33) in which v i , r i and m i are the velocity, position and mass of particle i, respectively. F i is the force acting on the particle i and V pot is the total potential energy of the system consisting of N particles.

The equations of motion could be integrated by breaking the computations into a series of very short time steps (normally about some femtoseconds). At each step t i , the forces on the atom i are calculated and combined with the current position and velocity to determine its new position and velocity at time t i+1 . An updated set of forces is then computed for the new position, and so on the trajectory of particle versus time is totally known.

Interatomic potential

The interatomic potentials play an important part in the MD simulations since they govern the dynamics of the system and thus the accuracy of the results. The best known of pairwise potentials is the Lennard-Jones (LJ) 12-6. It represent the van der Waals forces with one repulsive term in r -12 and one attractive term in r -6 . For a pair of atoms i and j the potential energy is :

V (r ij ) = 4ϵ [ ( σ r ij ) 12 - ( σ r ij ) 6 ] , r ij =| r i -r j | (1.34)
where the parameters ϵ and σ are the scales of energy and length, respectively. The interaction repels at close range, then attracts, and is eventually cut o at some limiting separation r c .

Traditionally the Lennard-Jones potentials are the most used since the mathematical expression allows to derive forces between particles eciently. However, although the attractive part depending on r -6 describes correctly the dispersion interaction between two particles, the repulsive part in r -12 has no theoretical justication. To overcome this drawback, the repulsive part can be replaced by V 0 e -α(r ij -R 0 ) which reproduces correctly the Pauli repulsion. The total pairwise potential should then be described by the Buckingham form [START_REF] Buckingham | The Classical Equation of State of Gaseous Helium, Neon and Argon[END_REF]:

V (r ij ) = V 0 e -α(r ij -R 0 ) - C 6 r 6 ij , (1.35)
where R 0 is the position of the minimum of V (r ij ). Increases in R 0 decrease the corrugation of the gas-surface potential. V 0 and α must allow to reproduce the Pauli repulsion at short ranges, and C 6 is the van der Waals dispersion parameter.

Besides, for a metallic cluster or a nanoparticle, the Sutton-Chen potential is often employed to calculate the interaction between solid atoms. It has the form :

V pot = ϵ     1 2 N ∑ i=1 N ∑ j=1 j̸ =i ( a r ij ) n -c N ∑ i=1 ρ 1/2 i     , ρ i = N ∑ j=1 j̸ =i ( a r ij ) m , (1.36)
where c is a dimensionless parameter and ρ i is the local density of atom i. The parameters ϵ and a determine the scales of energy and length, respectively and, n and m the range and shape of the potential.

To save computational time, the interaction potential is often truncated at a cut-o distance of r c ≈ 2.5σ to 5σ. This cut-o radius corresponds to the distance beyond which the interactions between the particles are negligible and not taken into account for the calculation of the forces. For example, at a typical cuto distance of r c = 2.5σ the interaction energy of LJ potential is just about 0.0163 of the well depth. However, the use of truncated potential can pose a problem at small cut-o radius. It may cause a discontinuity in both the potential energy and the force near the cuto value that eect to the energy conservation of system.

Boundary conditions

In order to conserve the macroscopic behavior of simulated system and save the computational time, the boundary condition should be treated in a reasonable manner. In Figure 1.6: Periodic boundary conditions for a simulation box in two dimensions molecular dynamic simulations, the most used method is periodic boundary conditions where the simulation region is divided into boxes and the central box is replicated in all directions to give a periodic array. An example of two dimensional box is plotted in Fig. (1.6). As seen in the gure, the central box is bounded by 8 neighbors; in three dimensions, the number of neighboring boxes would increase to 26. During the simulation, if a particle leaves the box, another particle having the same physical properties go through the opposite side. Hence, there are in fact no physical walls in boundaries of simulated box and the number of particles thus remains constant.

Another problem concerning the boundary condition is the interaction eld between the particles. It is important to ensure that the particle cannot interact with its own image and, hence, inuences its own behavior. This condition is known as the minimum image convention where the smallest box dimension must be greater than two times of the interaction range. In the minimum image convention, each particle detect at most just one image of every other particle of system which is repeated innitely by the periodic boundary conditions [START_REF] Leach | Molecular modelling: principles and applications[END_REF]. As the minimum image convention is respected, the periodic boundary conditions allow to maintain the energy and momentum of the system.

Initial state

Before a simulation can be performed, it is necessary to set initial positions and velocities to all particles in the system. As the aim of the simulation is to study the equilibrium uid state, the choice of the initial state should not have eects on the results of the simulation. Usually, the initial conguration of system can be started with simple lattices for the the gas atoms, such as the square or cubic lattice. On the other hand, the solid atoms are placed at the sites of a crystalline lattice, for example cubic, body-centered cubic, face-centered cubic, as shown in the Fig. 1.7. The lattice size should be chosen appropriately according to the type of system under study. The initial velocities are drawn from Maxwell-Boltzmann distribution corresponding to the average temperature of the system considered. The time integration of the Newton equations of motion allows the determination of the evolution of a system as function of time. We have the positions, velocities, and accelerations at time t and we want to calculate these values at time t + δt, where δt is the time step. A popular integration scheme which was widely used is the Verlet method and the Leapfrog algorithm [START_REF] Verlet | Computer "experiments" on classical uids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF]. The position of particle i at time t + δt is derived from the Taylor expansion of the position at time t :

r i (t + δt) = r i (t) + δtv i (t) + δt 2 2 a i (t) , (1.37)
where v i (t) and a i (t) are the velocity and acceleration of particle i at time t, respectively. The next step is to evaluate the velocity at half time (i.e. t + δt/2) and to calculate the acceleration at time t + δt from the interaction forces :

v i (t + δt 2 ) = v i (t) + δt 2 a i (t) , (1.38)
a i (t + δt) = - F i (t + δt) m i .
(1.39) Finally, the velocity at time t + δt is given by

v i (t + δt) = v i (t + δt 2 ) + δt 2 a i (t + δt) .
(1.40)

The process is continued with successive time step δt and the evolution of a system is therefore completely calculated. In molecular dynamic simulations, typical values of δt are close to 1 f s.

Thermostats

During the simulation, at equilibrium or non-equilibrium state, it is helpful to work at a constant temperature, at which a thermostat should be used. The thermostat can be applied either on wall particles, uid particles or all particles of the system. The thermostat methods are separated into global and local natures. The global thermostats, such as velocity scaling or the ones of Nosé-Hoover [START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] and Berendsen [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF], control the system temperature and dissipate the energy uniformly within the system. The local thermostats, including Andersen [START_REF] Andersen | Molecular dynamics simulations at constant pressure and/or temperature[END_REF], Langevin and Dissipative particle dynamics (DPD) thermostats, dissipate energy on the spatially localized scale from a stochastic point of view [START_REF] Yong | Thermostats and thermostat strategies for molecular dynamics simulations of nanouidics[END_REF].

In the present work, the Langevin thermostat is used for the solid walls. The motion of a particle i is governed by the equation

m i dv i (t) dt = -ξv i (t) + f i (t) + R i (t) .
(1.41

)
where ξ is a damping coecient, v i is the velocity and m i is the mass of particle i. f i and R i are the interaction force and the random force applied on the particle i, respectively. In statistical physics, one of the basic concept is the ergodic hypothesis, in which the average values over time of the physical quantities are equal to the statistical average values of the quantities. This hypothesis is expressed by

Q ≡ 1 τ ∫ τ 0 Q(t)dt = ∑ i p i Q i ≡ ⟨Q⟩ , (1.42)
where Q and ⟨Q⟩ are the time average value and the ensemble average value of the quantity Q, respectively. p i is the probability of nding the system in microscopic state i and Q i is the value of Q at microscopic state i. Physical systems for which the ergodic hypothesis is valid are called ergodic systems.

Macroscopic systems containing many particles are complex systems. In principle, to completely characterize a system with N particles, we have to know three components of position and three component of velocity for each particle. It is practically impossible to obtain these 6N real numbers. However, not knowing all the information about all particles does not prevent to calculate the macroscopic properties of the system. In fact, quantities like energy, temperature, pressure ... are computed as statistical averaged quantities of a large enough number of particles.

Dierent macroscopic environmental constraints lead to dierent types of ensembles, with particular statistical characteristics [START_REF] Kuzemsky | Theory of transport processes and the method of the nonequilibrium statistical operator[END_REF]. As the Newton equations lead naturally to energy conservation, molecular dynamics simulation is traditionally performed under the microcanonical or N V E ensemble where the number of particles N , volume V and energy E of the systems are kept constant. However, other ensembles are sometimes used in molecular dynamics as N V T ensemble, N pT ensemble and µV T ensemble. The ensemble widely used for Monte Carlo simulation is canonical or N V T ensemble in which the number of particles, volume, and temperature are xed. Generally, in large N V E systems, the uctuations of temperature are small, and temperature may be considered approximately constant. Therefore, a rather small part of a system may be considered as a canonical system.

Closely related to the N V T ensemble is the isothermal-isobaric, or N P T ensemble, with the pressure P , rather than the volume V , kept xed. While the microcanonical, canonical and isothermal-isobaric ensembles describe closed systems for which there is no change in the number of particles, the grand canonical or µV T ensemble is for an open system in which the number of particles can change and the chemical potential (µ), volume and temperature remain constant. Dierent from the kinetic theory, Molecular Dynamics methods describe explicitly all atoms in the system and the physical properties of the system can be computed accurately from their positions and velocities.

Energy

The internal energy is the sum of potential and kinetic energies

E = E p + E c .
(1.43)

For a system with N particles, the potential energy E p and kinetic energy E c are given by

E p = N ∑ i=1 E e i + V pot , (1.44)
E c = 1 2 N ∑ i=1 m i (v i -v) 2 , (1.45)
where E e i is the external interaction potential applied to particle i and V pot is the internal total potential of the system as mentioned at the beginning of the present Chapter. In (1.45), v i is velocity of particle i and v is barycentric velocity of ow.

Temperature

Unlike the canonical ensemble, in microcanonical ensemble, the temperature is not constant. It can be deduced from the kinetic energy of a system thanks to the theorem of the equipartition of energy as follows :

E c = k B T 2 (3N -N c ) , (1.46)
where N c is the number of constraints on the system and k B is the Boltzmann constant.

Pressure

For homogeneous uids, the external pressure is expressed via the virial theorem as the sum of two terms, one that corresponds to the part of ideal gas and the other involves the interactions between the particles (virial) :

P = 1 V   N k B T - 1 3 N ∑ i=1 N ∑ j=i+1 r ij f ij   , (1.47)
where r ij is the vector joining particles i and j, and f ij is the interaction forces between i and j.

For inhomogeneous uids, the pressure may be calculated from the formulae of Irving and Kirkwood [START_REF] Irving | The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics[END_REF], Todd et al. [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF] and the associated method of plane (MOP) by Evans et al. [START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF]. These issues will be dealt with in Chapter 2.

Optimization method for interaction computations

There are dierent methods to compute the interactions between particles, like : all-pairs, cell subdivision and neighbor list (or Verlet list). All-pairs method is easy to perform but is not ecient, at least when the interaction eld r c is small compared with the size of simulation region. Although the use of cuto radius allows to reduce the number of calculated forces, the test number of interatomic distances remains unchanged in comparison with the case without cuto radius (see Fig. 1.8). In the present work, we have combined cell subdivision and neighbor list methods in order to optimize the computational eort.

Figure 1.8: All-pairs method for interaction computations : in a simulation box, all pairs of particles should be examined.

Cell subdivision

This method consists in dividing the simulation region into many identical cubic cells whose side size is greater than the cuto radius. Thus a given particle interacts only with the particles located either in the same cell or in the adjacent cells. An example of cell subdivision method is shown in Fig. 1.9 in the 2D case. The computational zone of interactions is bounded by the central cell and 8 neighboring cells. In three dimensions, the number of adjacent cells of a given cell is increased up to 26. However, because of symmetry, only half of cells needs be examined for the interatomic forces computations (i.e. 13 in 3D and 4 in 2D).

Neighbor list

During the simulation, the microscopic environment change slowly : the particles almost remain in the same region of space for some integration steps. Therefore, in parallel with the cell subdivision method, we have used the neighbor list method to create the lists of neighboring particles with new radius r v . This radius is chosen such that a particle located outside the sphere bounded by r v radius does not have the time, during the rebuild time of neighboring list, to enter into the interaction zone of the studied particle. The radius of the neighbor list is written as :

r v = r c + ∆r .
(1.48) Figure 1.9: Representation of a simulation box in cell subdivision method in two dimensions.

Figure 1.10: Truncated sphere interactions related to r c and neighbor list related to r v , of a given particle.

The value ∆r is inversely related to the speed at which the list must be rebuilt, and it also determines the number of noninteracting pairs that are included in the list. It therefore has a certain inuence on both processing time and storage. Typically, for the fastest computation at liquid densities, ∆r ≈ 0.3σ to 0.4σ [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF].

Parallelization of the calculations

A molecular dynamic code can be written using two programming models: sequential or parallel. In the sequential model, a program is carried out by a single process. This process running on a processor of a machine and has access to processor memory. However, for complexes problems with important number of particles, the memory of a single processor may be not large enough and/or the computation time is too high. To overcome these problems, we used the parallel calculation of MPI type (Message Passing Interface). This allows to distribute computation charges on multiple processors.

In a parallel programming model by exchanging messages (MPI), each processor executes a copy of the calculation program and has access to its own memory. Therefore, the program variables become local variables at each processor. In addition, a processor cannot access the memory of neighboring processors. However, it may send/receive information to/from other processors by passing messages between processors via library functions MPI. The communication cost consists in two parts: the time to get started a message transfer, and a transfer time that depends on the message length. In order to make the message transfer working eciently, it is important that the communication cost should be low compared with the computation time of a processor. There is also the problem of charge balancing;

obviously the computation is optimal if all processors can be kept busy doing useful computing. If some processors have less work to do than the others, the overall eectiveness is reduced [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF].

The MD code used in this thesis is the parallel version described in Ref. [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. The original code has been developed to match the aim of the present works. For the problems studied, the ow eld is divided into several subregions which each processor handles. The number of subregions depends on the number of processors that will be used. In our simulations, the number of processors was only varied from 4 to 16 but it can easily be extended to higher numbers processors according to the machine architecture. Amongst the processors, there is a processor numbered 0 having the responsibility of collecting the values from all the other processors, after each of theme has calculated its local values. Processor 0 calculates then the total sums and send these values to each of the other processors.

The communication between two neighboring processors is shown in Fig. 1.13. The subroutines Detection detects if a particle leaves the eld of processor 1, it registers the particle number in Table A and then records the corresponding data as : position, velocities, accelerations of particle in Table B by subroutine P ackData. Since a particular particle is no longer associated with a xed memory location, so unlike the sequential version of the Figure 1.12: Spatial decomposition scheme in two dimensions : each cell is handled by a dierent processor. program, each particle is numerated with a unique identier. The subroutine M sgP ack sends the information stored in Table B to a buer called buf f send. This buer is often referred to as envelope because the data is being packed into a single message before transmission. After being sent and received between two processors via the MPI sending and receiving functions, the information is stored in the buf f receive. Then unpack the data, we nally have the information particle which just enters into the neighboring processor 2.

+D=FJAH

Molecular dynamics simulation of Poiseuille ows

Introduction

One of the fundamental and practical problems in uid mechanics is the Poiseuille ows between two parallel plates. In the case of Newtonian uids, the problem has a simple analytical solution which has been used extensively as benchmarks for numerical methods, constitutive models and studies of boundary conditions (see [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF] and the references cited therein). Generally, ows and related simulation methods can be classied as acceleration driven or pressure driven. Conventional MD simulations with periodic boundary conditions [4851] and gravity-like force eld are typical examples of the former class. Regarding the latter case, dierent techniques have been developed to mimic ows induced by pressure dierence [5254]. Lupkowski and Swol [START_REF] Lupkowski | Computer simulation of uids interacting with uctuating walls[END_REF] placed two rigid walls at the inlet and outlet and applied a piston external forces on them. Li et al. [START_REF] Li | Coupling continuum to molecular-dynamics simulation: Reecting particle method and the eld estimator[END_REF] used a ctitious membrane which allows atoms to pass from one direction and forces atoms from the other direction to be elastically reected with a given probability. Sun and Ebner [START_REF] Sun | Molecular-dynamics simulation of compressible uid ow in two-dimensional channels[END_REF] simulated 2D uid ows in a long box while controlling the temperature at the two ends.

The present chapter presents a simple alternative procedure to simulate pressure driven ows with Molecular Dynamics method. Section 2.2 gives the problem statement, presents some approximate analytical solutions reported in the archival literature for Poiseuille ows. The pressure expression and the pressure dierence from the atomistic point of view are also discussed in Section 2.3. Section 2.4 is focused on the implementation of periodic boundary conditions linked to the pressure dierence concept, and its implementation in a MD code. The results of simulations are then compared with approximate analytical solutions of the Navier-Stokes and energy equations in Section 2.5. The application of the present method to uids other than ideal gases and to rib-roughened channels are also discussed.

Poiseuille ows and analytical solutions

In the Cartesian coordinate system Oxyz, we assume that the ow direction is along the x-axis. The length, width and height of the channel are denoted by L, B and H in direction x, y and z, respectively. The channel width B being very large with respect to the other dimensions, it will be treated as innite. The uid can exchange energy with the channel walls which are maintained at constant temperature T w .

The range of Kn chosen in the present study, 0.02 ≤ Kn ≤ 0.1, belongs to the slip ow regime according to [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF], i.e the uid continuum mechanics is still useful provided that slip boundary conditions are applied at the walls. The validity of the Navier-Stokes equations, as shown by many authors (see [START_REF] Bhattacharya | Molecular dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases[END_REF] for example), can therefore be extended to Poiseuille and Couette ows for the present range of Kn.

In what follows, the terms "pressure driven ow" and "acceleration driven ow" are employed to dene two particular cases: i) the body force ργ applied on a unit uid volume is zero and the pressure p decreases along the ow direction, ργ = 0, dp/dx ̸ = 0 ii) the inverse situation where ργ ̸ = 0, dp/dx = 0. We shall also consider some available analytical solutions of the Navier-Stokes equations reported in the literature. For both ows, the uid velocity v in a cross section which accounts for slip eects (Navier boundary conditions) at the wall admits the parabolic form

v(ẑ) = v(0) [ 1 -a 2 ẑ2 ] , ẑ = z H . (2.1)
where a 2 is a positive dimensionless constant related to pressure gradient (or body force), viscosity, slip length, etc... It is noted that Eq. 2.1 is generally obtained for incompressible uids but it is also valid for compressible uid ows in a long channel [START_REF] Arkilic | Gaseous slip ow in long microchannels[END_REF]. In the latter case, all the quantities v, a 2 are functions of the streamwise coordinate x.

An approximation solution of the energy equation can be obtained by neglecting the convective term [START_REF] Landau | Fluid mechanics[END_REF]. The temperature T becomes then a quartic function of the coordinate

z T (z) = T w [ b 0 -b 4 ẑ4 ] .
(2.2) Equation 2.2 is based on the incompressibility assumption. For compressible uid ows in a long channel, Cai et al. [START_REF] Cai | Gas ows in microchannels and microtubes[END_REF] used a perturbation technique to derive the following temperature distribution

T (x, z) = T w [ c 0 -c 2 ẑ2 -c 4 ẑ4
] .

(2.

3)

The dimensionless coecients c i , b i in Eq. 2.2 and Eq. 2.3 are related to pressure gradient (or body force), viscosity, conductivity, specic heat etc... In Eq. 2.3, the coecients c i depend on the position x. We also note that for acceleration driven ow, Todd and Evans [START_REF] Todd | Temperature prole for poiseuille ow[END_REF] suggested a correction to the temperature prole given by Eq. 2.2 as a sextic function of z:

T (z) = T w [ d 0 -d 2 ẑ2 -d 4 ẑ4 -d 6 ẑ6
] .

(2.4)

In section 2.5, MD simulations are used to reexamine the validity of the velocity and temperature proles given by Eq. 2.1-2.4 and to determine the coecients a i , b i , c i , d i by curve tting.

Pressure tensor and pressure dierence

Before proceeding the molecular dynamics simulation, let us look at the denition of pressure tensor from the atomistic viewpoint. Using statistical mechanic theory, Irving and Kirkwood [START_REF] Irving | The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics[END_REF] derived the following pressure tensor decomposition (IK):

P = P K + P V .
(2.5)

The kinetic term P K depends on the molecular square velocity while the potential term P V depends on pairwise interactions between molecules. In a system of N molecules, P K and P V read

P K (r) = N ∑ i=1 ⟨m(v i -v) ⊗ (v i -v)δ(r i -r)⟩ P V (r) = - 1 2 N ∑ i N ∑ j̸ =i ⟨ r (ij) ⊗ f (ij) O ij δ(r i -r) ⟩ (2.6)
The term inside the angular bracket ⟨.⟩ denotes the ensemble average and ⊗ denotes the dyadic product. The term r i , v i , v are respectively the position vector, the velocity of the particle i and the mean velocity. The distance vector and interaction force vector between two molecules (i) and (j) are denoted respectively as r (ij) and f (ij) . For LJ uids, the two vectors r (ij) and f (ij) are colinear and the force magnitude f (ij) is derived from the interaction potential in function of the distance between the two molecules

f (ij) = ∂V(r (ij) ) ∂r (ij) , r (ij) =∥ r i -r j ∥ .
(2.7)

In Eq. 2.6, the notation δ(r i -r) is used for the delta Dirac function at r i and the expression of the operator O ij reads

O ij = 1 - 1 2! r (ij) ∂ ∂r + ... + 1 n! [ -r (ij) ∂ ∂r ] n-1 + ... (2.8)
Equations 2.8 and 2.6 show that the derivation of IK pressure tensor involves an innite sum of high order derivatives of the delta function and ensemble average, not suitable for MD computations. A more convenient form of the pressure tensor and an associated calculation method, the method of plane (MOP), was proposed in Todd et al. [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF] and Evans and Morriss [START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF]. When the uid density is uniform, the operator O ij is reduced to unity and equation 2.6 becomes (see [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF][START_REF] Hoover | Microscopic and macroscopic stress with gravitational and rotational forces[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF]):

P K (r) = 1 V ∑ i m(v i -v) ⊗ (v i -v), P V (r) = - 1 V ∑ i<j r (ij) ⊗ f (ij) (2.9)
where V is the volume of the uid element located at r. In the literature, Eq. 2.9 is sometimes referred to as IK1 pressure [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF][START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF]. In Eq. 2.6 and 2.9, the pressure can be decomposed into two scalar quantities p K and p V as follows

p = p K + p V , p K = 1 3 tr P K , p V = 1 3 tr P V .
(2.10)

It is clear that for ideal gases, the potential part p V vanishes and the pressure is reduced to the kinetic part p K . For the other uids, the potential part p V cannot be neglected. The pressure dierence ∆p between two points at distance ∆l reads

∆p = ∆p K + ∆p V = k B ∆(nT ) + ∆p V .
(2.11)

The pressure components P xx , P yy , P zz , number density n and stream velocity v are invariant in direction y. From a microscopic point of view, taking ∆l = ∆y, we must have

∆⟨v 2 x ⟩ = ∆⟨v 2 y ⟩ = ∆⟨v 2 z ⟩ = 0, ∆⟨v x ⟩ = ∆⟨v y ⟩ = ∆⟨v z ⟩ = 0 .
(2.12)

As it will be discussed in the following, all these conditions are satised in our MD simulation algorithm if we apply periodicity for the boundary conditions in the y direction. However, n, p, T can all vary in the ow direction x. We note that in the steady state regime, all quantities like T , p, p V , n are stable with time, and thus their dierences ∆T , ∆p, ∆p V , ∆n are also stable. In particular, equation (2.11) shows that the nite variation of temperature ∆T along x is related to those of other quantities via the expression

∆T = ∆p -∆p V -k B T ∆n k B n + k B ∆n .
(2.13)

If the stream velocity is much smaller than the thermal velocity, we can obtain a simplied expression for the temperature dierence as

∆T = m 3k B ∆⟨v 2 -v 2 ⟩ ≃ m 3k B ∆⟨v 2 ⟩ .
(2.14)

The objective of the MD method discussed in what follows is to maintain the dierences in squared velocity and, thus in temperature, in direction x. By that way, we can indirectly generate the pressure dierence ∆p. In the case where density ρ is invariant in x-direction (incompressible ow assumption), the temperature dierence is proportional, for an ideal gas, to the pressure dierence since

∆p = ρ 3 ∆⟨v 2 ⟩ .
(2.15)

Modied boundary conditions

In Molecular Dynamics, periodic boundary conditions (PBC) applied to velocities are techniques to reduce the study of a large system to a smaller one far from the edge. Considering a simulation domain as a cube, PBC requires that if a molecule passes through one face, it reappears on the opposite face with the same velocity. Obviously, there is no dierence in pressure, density or temperature between any two opposite faces of the simulation domain.

Consequently ∆p = 0, as it will be shown in the next section.

In our problem, all the pressure components P xx , P yy , P zz , do not vary along y, which can be satised by the traditional PBC applied on the faces y = ±B/2. However to create a pressure dierence along the x-axis, we must develop a strategy that produces a constant dierence between the squared velocity of the molecules crossing the face x = 0 and those crossing the face x = L.

Usually, to keep unchanged the total number of molecules in the domain whenever a molecule goes out of the domain, we must insert another molecule inside. In this chapter, we generalize the PBC to account for the pressure dierence by maintaining the dierence in squared velocity at the two opposite faces x = 0 and x = L.

From equations 2.13 and 2.14, we know that a constant pressure dierence is related to a constant dierence in squared velocity. A constant parameter δ is thus used to create a constant dierence in squared velocity between the inlet and outlet faces (see Fig. 2.1).

In order to apply a velocity square dierence equal to a constant 3δL between the faces x = 0 and x = L, we modify the periodic boundary conditions so that if a molecule goes trough the face x = L with a velocity v, we insert a molecule at x = 0 with a velocity v ′ satisfying the following conditions Analogously, if a molecule crosses the face x = 0 with a velocity v, it reappears at x = L with a velocity v ′ so that

v ′ = αv, α = √ 1 -3δL/v 2 .
(2.17)

Since 3δL is positive, it may happen that v 2 -3δL < 0. It is then impossible to nd v ′ Figure 2
.1: Modied periodic boundary conditions. The velocity v between the outgoing molecules are dierent from the ingoing molecules (Eq. 2.16 and 2.17) Thermal walls satisfying Eq. 2.17. However, this eect can be assumed negligible for ows with relatively small speed (of the order 1% of the thermal speed), as seen in most MEMS/NEMS devices.

In the present simulations, the number of molecules that do not satisfy Eq. 2.17 were found less than 1% of the outgoing molecules. In order to simulate low speed ows, small values of δ were used.

2.5 Numerical results

Simulation of ideal gas ows

To model the interaction between the molecules, the 6-12 Lennard Jones potential is used:

V (r) = 4ϵ [ ( σ r ) 12 - ( σ r ) 6
] ,

(2.18

)
where ϵ is the depth of the potential well, σ is the nite distance at which the inter-particle potential is zero, and r is the distance between the particles.

The potential well depth parameter is ϵ = 0.1ϵ 0 with ϵ 0 being the reference potential well depth. Two cuto distances were used in this study : r c = 2 1/6 σ [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF] corresponds to repulsive interactions between molecules while r c = 2.5σ [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] corresponds to attractive-repulsive interactions. In general, the present method is valid for any pair potential. However, the simulation conditions of Table 2.2, i.e. (ϵ = 0.1ϵ 0 , 2 1/6 σ), lead to a very small potential part of the pressure p V (less than 1 % of the kinetic part p K ) making the behavior of the uid close to an ideal gas.

Table 2.2: Input data for the rst series of MD simulations. The two walls are modeled as thermally diusive walls at the same and constant temperature. Whenever a molecule arrives at a wall, it is reected with the velocity corresponding to the wall temperature T w = 2ϵ 0 /k B and with a random direction. In this work, we do not use thermostats, the uid can exchange energy with the wall due to the wall model described previously.

N L B H Kn δ γ [σ] [σ] [σ] [ϵ 0 /mσ] [ϵ 0 /mσ]
In our simulations, the global number density n is kept xed, n = 0.1σ -3 in every congurations. Other geometric parameters like L, B, H and the number of the molecules N are changed as shown in Table 2.2. The values of δ are chosen such that the dierence of the square of the molecular velocity between the inlet and outlet is small enough so that the stream velocity is much smaller than the thermal speed. From the microscopic viewpoint, the local stream velocity is the local average velocity of the molecules and from the macroscopic viewpoint, it is the ow speed v in the Navier Stokes equations and thus connected to the pressure gradient ∇p (or volume force ργ). On the other hand, the thermal speed is the root mean square of the molecular velocity and is directly connected to the local temperature T . The slope of the simulated total pressure distribution along of the channel (see Fig. 2.2 and 2.3) is used to determine the pressure gradient -∆p/∆x and the equivalent external acceleration ργ used in the accelerationdriven ow simulations. The quantity ρ is the average mass density of the system, ρ = mn. The computation with Leap-Frog Verlet integration scheme is carried out for 10 8 time steps from equilibrium, each of which is equal to 0.005 unit time τ = √ mσ 2 /ϵ 0 . The height and length of the channel are divided into 50 layers to determine accurately the distribution of local velocity v, temperature T and pressure tensor P. The local velocity v measured in one bin located at r is given by

v(r) = 1 n bin ∑ r i ∈bin v i , (2.

19)

n bin being the particle number in the bin. The local pressure P (IK1-model) is:

P(r) = 1 V bin   ∑ r i ∈bin m(v i -v) ⊗ (v i -v) - 1 2 ∑ r i ∈bin r (ij) ⊗ f (ij)   .
(2.20)

The pressure components were also computed by using the method of plane [START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF][START_REF] Han | Method for calculating the heat and momentum uxes of inhomogeneous uids[END_REF]. With α, β being x, y or z, the pressure component P βα along direction β and acting on the area element A α normal to the axis α is dened by:

P βα (r) = 1 τ A α ∑ r i ∈Aα m(v iβ -viβ )(v iα -viα )/v iα + 1 τ A α ∫ dt ∑ r ij ∩Aα,i<j f (ij) β . (2.21)
In Eq. 2.21, the rst sum is for all molecules crossing the area element A α over the period τ , and the second sum is for all pairs whose distance vectors cut the area element A α . The temperature T of one bin is then calculated by

T (r) = m 3k B ⟨∥ v -v ∥ 2 ⟩ bin .
( along the channel for pressure driven ow case (Kn = 0.05,

r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ). The pressure unit is [ϵ 0 /σ 3 ].
In what follows, the x-distributions of pressure, temperature and density correspond to an average on y and z of the local quantities p(r) and n(r). They are denoted as p(x) = ⟨p(r)⟩ y,z and n(x) = ⟨n(r)⟩ y,z , respectively. Analogously, the prole along z of velocity, temperature and density correspond to an average on x and y of local quantities with the notations: v(z) = ⟨v(r)⟩ x,y , T (z) = ⟨T (r)⟩ x,y and n(z) = ⟨n(r)⟩ x,y , respectively.

From the simulations, the axial pressure distribution is plotted in Fig. 2.2 for the pressure driven and the acceleration driven ows with Kn = 0.05, r c = 2 1/6 σ. As expected, for the acceleration driven ow case, the pressure is constant. In contrast, the pressure distribution for pressure driven ow decreases linearly. We note here that the slope of the pressure curve, equal to ∆p/∆x, is used to compute the value γ in the acceleration ow simulation according to the formula γ = -∆p/(ρ∆x). Under the simulation conditions n = 0.1σ -3 , T w = 2ϵ 0 /k B and r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , the kinetic part p K of the pressure dominates. In order to know how the pressure dierence is distributed in the x, y, z-directions , the pressure components Pxx , Pyy , Pzz along x are plotted in Fig. 2.3. All the three pressure components decrease with x and the slopes are almost the same. The values of components Pxx , Pyy are very close whereas the value of component Pzz is smaller because the uid is conned in z-direction. Although the dierence between the pressure components is small (less than 4%), it implies that the velocity distribution deviates from the equilibrium distribution. Because the pressure p is the average value of Pxx , Pyy , Pzz , the p-curve lies between the others and have the same slope. It should be noted that equations 2.20 and 2.21 leads to relatively close results when considering the average axial pressure component Pxx (see Fig. 2.4).

v(ẑ) = v(0)(1 - a 2 ẑ2
) that best t the numerical results (a 2 = 0.18 for Kn=0.05 and a 2 = 0.29 for Kn=0.1).

channel width varies slightly around the average density n, except near the walls where it can be as high as 1.4n. However, the phase diagram for LJ uid [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] shows that the uid near the wall is still in gas state. In our simulations, n and ϵ are chosen to be very small (n = 0.1σ -3 , ϵ = 0.1ϵ 0 ) and the uid temperature near the wall is relatively large (the temperature of the wall is T w = 2ϵ 0 /k B ). The uctuation of the density prole shows that the incompressibility assumption is no longer valid, especially near the channel walls. However, the uctuations being localized, it is interesting to check if the analytical solutions presented in Section 2 still agree with the MD solutions.

Figure 2.6 shows that acceleration and pressure-driven ows exhibit parabola-like velocity proles in the center region of the channel, in agreement with Eq. 2.1. Near the wall where the Knudsen layer dominates, the velocity distribution tends to deviate from the solution given by Eq. 2.1 and a velocity slip at the wall is predicted. Based on the kinetic theories, dierent models have been derived to explain the slippage [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF] and generally, the slip eect becomes important when the Knudsen number increases. The little dierence in velocity prole between the two types of ow show that, despite their dierent microscopic natures, volume forces can be seen as equivalent to pressure gradients at the macroscopic scale. Regarding the temperature prole across the channel width, there is no visible dierence between acceleration and pressure-driven formulations whatever Kn (Fig. 2.7). The temperature prole cannot be tted by any of the three approximate analytical solutions (Eq. 2.2-2.4). The temperature is minimal at the center of the channel and increases rapidly near the wall, which does not agree with the approximate solutions obtained with the incompressible assumption (Eq. 2.2-2.3). With the simulation conditions reported in Table 2.2, this reverse trend cannot thus be explained by using the incompressible ow equations. This anomaly was also observed by using DSMC method [START_REF] Zheng | Comparison of kinetic theory and hydrodynamics for poiseuille ow[END_REF] and Super-Burnett (SB) method [START_REF] Xu | Super-burnett solutions for poiseuille ow[END_REF]. Note that the ow speeds in our simulations are very small, as in [START_REF] Zheng | Comparison of kinetic theory and hydrodynamics for poiseuille ow[END_REF][START_REF] Xu | Super-burnett solutions for poiseuille ow[END_REF].

Next, we study the inuence of γ and δ on the temperature prole. By increasing the acceleration parameter γ or the δ-parameter (squared velocity,) the stream velocity is increased. A change in the form of the temperature prole is observed in Fig. 2.8 and 2.9 as γ (or δ) increases. At high values of γ (or δ), the temperature proles seem to be closer to the approximate solutions which predict a maximum at the center of the channel and a minimum at the walls. At small values of γ (or δ), the temperature proles for the two ow types dier just a little. However, at high values of γ (or δ), the discrepancy becomes more important: the curvature at the center of channel for the pressure driven ow case is higher than that for the acceleration driven ow case. Using curve tting procedure, we nd that the temperature prole of the pressure driven ow case agrees quite well with the quartic expression (Eq. 2.3). The temperature for the acceleration driven ow case does not t well with Eq. 2.2. To explain this discrepancy, Todd and Evans [START_REF] Todd | Temperature prole for poiseuille ow[END_REF] argued that the transport coecients are not constant and there exist an additional cross-coupling between strain rate and heat ux. They proposed a correction to (Eq. 2.3) with a z 6 -extra term for the temperature prole (see Eq. 2.4). The dashed line drawn in Fig. 2.9 shows that the sextic polynomial t very well the present numerical results.

Applications to general cases

The application of the method to cases for which the uid particles interact strongly, i.e the potential pressure p V is of the order of p K , is considered in that section. The case of rib-roughened channels is also briey discussed.

With the algorithm developed in the previous section, two cases were considered with the same parameters, except for the cuto distance r c and density number n. In the rst case, n = 0.1σ -3 and the cuto distance is set to r c = 2.5σ (attractive-repulsive interaction). In the second case, the uid density is increased up to n = 0.8σ -3 while seting r c = 2 1/6 σ as in the previous section. The well depths ϵ are the same in both cases and equal to the reference value, i.e ϵ = ϵ 0 . All the parameters for these two cases are summarized in Table 2

N L B H δ T w [σ -3 ] [σ] [σ] [σ] [ϵ 0 /mσ] [ϵ 0 /k b ]
0.1 2400 86. In the rst case (n = 0.1σ -3 ), gure 2.10 displays the kinetic part pK and the virial pressure p according to Eq. 2.9 along the channel. Except a slight curvature in the pK curve, both pressure curves decrease continuously and linearly. The relative dierence between the slopes of both curves is less than 1%, which means that the gradient of the potential part ∆p V /∆x is also small in comparison with the kinetic part ∆p K /∆x even if pV is not negligible (|p V /p K | ≃ 0.2). In the second case (n = 0.8σ -3 ), Figure 2.11 shows a totally reverse trend: both the virial pressure prole and its potential part decrease with x.

The variation of the potential part contributes to 80% of the virial pressure variation with ∆p/∆x = -0.0177 and ∆p V /∆x = -0.0159. We conclude that the method is convenient for generating virial pressure proles which decrease in the axial direction.

Next, we consider the velocity proles across the channel section for these two cases. The velocity proles (Fig. 2.12) agree then with Eq. 2.1. The slip length for the case n = 0.8σ -3 is smaller than for the case n = 0.1σ -3 because the mean free path is decreased when the density number is increased. When computing the average pressure component Pxx along the channel, gure 2.13 shows that there is no visible dierence between equations 2.20 and 2.21 Figure 2.12: Dimensionless velocity prole in half of the channel cross-section for n = 0.1σ -3 and n = 0.8σ -3 . The dashed lines represent the equations a 0a 2 z 2 that best t the numerical results.

It should thus be emphasized that the present method is relevant to generate various uid ows even if we control only the dierence in squared molecular velocity. It does not imply that the potential part p V , which depends on the interatomic interaction, does not vary along the ow direction. As shown numerically, the variation of p V contributes considerably to the pressure gradient. The third case aimed at modeling the ow within a rib-roughened channel. The parameters are the same as in the second case (n = 0.8σ -3 , r c = 2 1/6 σ ) but the channel height is suddenly reduced in its middle part by inserting on both walls ribs of height h and length l.

The rib-to-channel width ratio is equal to h/H = 0.8 and l/L = 0.3 (Fig. 2.14). Figure 2.14 exhibits the onset of two vortices close to the rib corners at the upstream section. It is also shown that the uid ow is highly non uniform and characterized by wavy streamline patterns within the downstream region. From Fig. 2.15, the axial pressure variations are rather dierent from those predicted for smooth channels. Although the variations of the axial pressure predicted by the two expressions are quite similar, there are considerable dierences between them. The MOP pressure is 8% smaller than the IK1 pressure. The changes in pressure at the downstream and upstream sections display also smoother pressure variations. This considerable dierence is due to the strong inhomogeneity of the uid within the channel. Generally, the MOP, based on the hydrodynamics equation, must be used in such situation. +D=FJAH !

Gas-wall potential and eects on the accommodation coecient

Introduction

As mentionned in Chapter 1, Maxwell introduced a gas-wall momentum accommodation parameter σ t (TMAC) to quantify the slip eects [START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF]. He postulated that after collision with the wall, a gas atom rebounds either diusively or specularly, with the associated portions of σ t and 1σ t , respectively. Although Molecular Dynamics simulations showed that the reection mechanism is more complicated than Maxwell's postulate, the coecient σ t is still widely used due to its simplicity. For example, in the case of isothermal ows, Eq. (1.28) is reduced to

U slip = 2 -σ t σ t Kn ( ∂U ∂n ) s . (3.1)
With MD method, the accommodation parameter can be computed either indirectly or directly. The former concerns the simulation of conned ows [START_REF] To | A slip model for micro/nano gas ows induced by body forces[END_REF][START_REF] Cao | Temperature dependence of the tangential momentum accommodation coecient for gases[END_REF], nding the slip velocity and using slip equation (3.1) to calculate σ t . The latter method consists in projecting a gas atom with velocity v i into a surface, measuring the reected velocity v r , and computing σ t [START_REF] Finger | Molecular dynamics simulation of adsorbent layer eect on tangential momentum accommodation coecient[END_REF]. However, the success of the MD method relies on the accuracy of the interaction potential of the gas-wall couple into consideration. The choice of such potential for a MD simulation is crucial. This interaction potential can be either determined by quantum electronic structure calculations or elaborated to reproduce experimental data with some empirical parameters.

Based on periodic density functional theory computations, a Pt-Ar potential is developed for MD simulations. Other pairwise potentials from the literature will be also considered to discuss what kind of interaction potentials should be used for correct MD simulations and what are the consequences of the choice of a given potential.

We have chosen the Pt-Ar couple because a number of references exist about the interaction potential between the Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface and an Ar atom. In particular, the Pt(111)-Ar couple has been widely involved in MD simulations of atomic beam experiments [6774] or of slip ows [START_REF] Cao | Temperature dependence of the tangential momentum accommodation coecient for gases[END_REF]7578], as examples. A section of the present chapiter is devoted to a brief review of the existing Pt-Ar pairwise potentials. Experimentally, only few data exist about the interaction of an Ar atom with a Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface. Zeppenfeld et al. [START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF] measured the vertical vibrational frequency of the adsorbed atom at about 5 meV from inelastic He scattering on Ar monolayer physisorbed on Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF].

To the best of our knowledge, the rst reference to the total adsorption energy appeared in the work of Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] who reported a Ar on Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] binding energy of about 80 meV (7.5 kJ.mol -1 or 900 K). Zeppenfeld et al. [START_REF] Zeppenfeld | Structure of monolayer Ar on Pt (111): Possible realization of a devil staircase in two dimensions[END_REF] mentioned a value of 96 meV from unpublished results of Zeppenfeld, Kern and Becher. Most of the empirical potentials have been optimized in order to reproduce the 80 meV binding energy and experimental results on trapping, desorption and scattering of Ar on Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF].

In the following, we give a brief description of the previous interaction potentials between Ar and Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF]. r represents the position of the Ar atom, r i the position of the ith Pt atom and z is the normal distance of Ar above the surface. Most of the potentials given below are empirical and Table 3.1 contains the associated parameters.

Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] have employed Morse type functions to reproduce the experimental measurements about the trapping and desorption of the Ar atom on the Pt surface, as well as the binding energy of 80 meV. They assumed that the incorrect form of the attractive part of this potential does not play a crucial role for trapping or scattering:

V Head-Gordon = D e ∑ i (exp[-2σ(|r -r i | -R 0 )] -1) 2 , (3.2)
where D e is the dissociation energy parameter which controls the strength of binding and σ is the strength parameter which is linked to the eciency of energy transfer between gas and surface. It must be noticed that Yamamoto [START_REF] Yamamoto | Slip ow over a smooth platinum surface[END_REF] used this Morse potential for the study of the slip ow of Ar on Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF], but with a twice too large value for σ.

The most common empirical potentials used for MD are Lennard-Jones 6-12 potentials, in order to obtained a better representation of the long-range attractive part that can be important in the simulation of the fragmentation of clusters or droplets:

V LJ = 4ϵ ∑ i [ ( σ |r -r i | ) 12 - ( σ |r -r i | ) 6 ] , (3.3)
ϵ is the depth of the potential well and σ is the nite distance at which the pair potential is zero. As examples, Svanberg and Pettersson [START_REF] Svanberg | Survival of noble gas clusters scattering from hot metal surfaces[END_REF] have optimized parameters such that the interaction potential reproduces the one of Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] in order to model the scattering of Ar clusters from the Pt(111) surface. Maruyama and Kimura [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF], and Spijker et al. [START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF] have suggested values for σ and ϵ for the simulation of nucleation of Ar liquid droplets on Pt surface. Ramseyer et al. [START_REF] Han | Thermally-driven nanoscale pump by molecular dynamics simulation[END_REF] have also derived a Lennard-Jones type potential from the experimental data derived by Zeppenfeld et al. [START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF][START_REF] Zeppenfeld | Structure of monolayer Ar on Pt (111): Possible realization of a devil staircase in two dimensions[END_REF].

More elaborated interaction potentials have also been derived. For a MD study of the trapping and scattering of Ar on Pt(111), Smith et al. [START_REF] Smith | A molecular dynamics study for the trapping and scattering of Ar/Pt(111)[END_REF] used an interaction potential between the Ar atom and the Pt surface dened as a sum of pairwise repulsive terms and a non-corrugated van der Waals attractive part:

V Smith = ∑ i V 0 e -α(|r-r i |) - -C vdW √ (z -z 0 ) 6 + (C vdW /V min ) 2 .
(3.4)

The parameters were chosen such that C vdW = 2.3697 eV.Å 2 , evaluated from the Ar polarizability and the Pt dielectric function. Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF] determined an empirical potential energy function for the Pt(111)-Ar system based on atomic beam scattering data, the measured well depth, and the vibrational frequency of the absorbed atom by Comsa et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF][START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF]. This function is composed of a pairwise sum of non-central potentials. The repulsive pair potential has the Buckingham form [START_REF] Buckingham | The Classical Equation of State of Gaseous Helium, Neon and Argon[END_REF]:

V Kulginov = ∑ i ( V 0 e -α(|r-r i |-r 0 ) -D(r) C 6 |r -r i | 6 ) , (3.5)
and the attractive part corresponds to the van der Waals attraction multiplied by a damping function D(r). C 6 is chosen such that the total potential reproduces the van der Waals attraction calculated by Smith et al [START_REF] Smith | A molecular dynamics study for the trapping and scattering of Ar/Pt(111)[END_REF]. They suggested also that their and the Head-Gordon functions [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] should by shifted by ≈ -0.7 Å to better agree with the sum of the van der Waals radii of the Ar and Pt atoms. They conrmed that the position of the potential well can not be determined by the scattering data alone and that electronic structure calculations should be made. From the ab initio point of view, only Lahaye [START_REF] Lahaye | Site dependent energy loss in Ar scattering from Pt(111)[END_REF] have used Hartree-Fock calculations on Ar approaching a Pt 4 cluster to derive a Pt(111)-Ar interaction potential with a similar form to the one of Smith et al. [START_REF] Smith | A molecular dynamics study for the trapping and scattering of Ar/Pt(111)[END_REF] (see equation (3.4)).

Computational approach

As explained in the previous section, several potentials exist for the Pt(111)-Ar interaction, most of them being obtained empirically since the computation of such interaction potential still remains a challenging task. In order to ll this lack, large periodic computations using the density functional theory approach have been performed. All these electronic structure computations were done with the CRYSTAL09 code [START_REF] Dovesi | CRYSTAL: a computational tool for the ab initio study of the electronic properties of crystals[END_REF] using the PBE functional revised for the study of solids (PBEsol) [START_REF] Perdew | Restoring the density-gradient expansion for exchange in solids and surfaces[END_REF]. For the description of the Pt atoms, the pseudo-potential and basis set introduced by Doll [START_REF] Doll | CO adsorption on the Pt (111) surface: a comparison of a gradient corrected functional and a hybrid functional[END_REF] to investigate the absorption of CO on a Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface was used. For the Ar atom, the all electron basis set developed by Nada et al. [START_REF] Nada | Basis sets for ab initio periodic Hartree-Fock studies of zeolite/adsorbate interactions: He, Ne, and Ar in silica sodalite[END_REF] was chosen. As a starting point of these computations, the lattice parameter of a platinum bulk was investigated to check the accuracy of the methodology used. At this level of calculation, the equilibrium lattice parameter is obtained at a = 3.96 Å. That is in relatively good agreement with the experimental value of 3.92 Å [START_REF] Gschneidner | Physical properties and interrelationships of metallic and semimetallic elements[END_REF], and consistent with those computed by Doll [START_REF] Doll | CO adsorption on the Pt (111) surface: a comparison of a gradient corrected functional and a hybrid functional[END_REF] at 4.01 Å and 4.05 Å using the PW91 and B3LYP functionals, respectively.

The pairwise potential cannot be obtained directly from quantum chemistry calculations since the interaction of an Ar atom with an isolated Pt atom is strongly dierent of an interaction of an Ar atom with an Pt atom embedded in a surface slab. A way to obtain the pairwise interaction potential between a gas atom and a surface atom is to decompose of the Pt(111)-Ar potential in two bodies term. To this end, the Pt(111)-Ar interaction was modeled by using a super-cell representing a three-layer Pt slab, the Pt atoms being arranged as in the face-centered cubic lattice, with the PBEsol optimized platinum bulk lattice parameter of 3.96 Å. The Ar atom was adsorbed on one side of the slab, vertical to the [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface in front of a fcc-hollow site. On this site, the Ar atom is placed vertically above a Pt atom of the third Pt layer. We have checked that the adsorption energy is little sensitive to the adsorption site, the variation in the adsorption energy being of 7 meV between the top and fcc-hollow sites. Several slab sizes have been tested. The (4x4x3) slab appears to be a good compromise between accuracy of the interaction energy and computational eort. Using this (4x4x3) slab, the distance between two Ar atoms is about 11.2 Å, about three times larger than the equilibrium distance of the diatomic Ar 2 molecule (R e (Ar-Ar) = 3.76 Å [START_REF] Colbourn | The spectrum and ground state potential curve of Ar2[END_REF]). The relaxation of the rst layer was not considered according to the work of Doll [START_REF] Doll | CO adsorption on the Pt (111) surface: a comparison of a gradient corrected functional and a hybrid functional[END_REF] showing the very small change in the binding energy due to such a relaxation. Finally, a k-point sampling net of size 16 x 32 was used, and the Fermi function was smeared with a temperature of 0.01 E h to improve the convergence of the total energy with respect to k-point sampling and the convergence of the Pt(111) surface energy with respect to the number of layers [START_REF] Dovesi | CRYSTAL: a computational tool for the ab initio study of the electronic properties of crystals[END_REF].

The resulting computed energies have been tted by a sum of pairwise functions in order to reproduce the Pt(111)-Ar interaction potential.

V P t(111)-Ar (z) = ∑ i V P t-Ar (|r -r i |) , (3.6)
z is the distance between the Ar atom at the vertical of the fcc-hollow site and the rst layer of the Pt(111) surface. z = 0 is dened as the position of the Pt nuclei of this rst layer for T = 0K. The sum is made over all of the 121 Pt atoms of the cluster depicted in Figure 3.2.2 which is supposed to reproduce the Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface. According to the previous empirical pairwise potentials, a similar function (see equation (3.7)) as the one used by Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF] has been employed to t the computed energies. The tted parameters are given in Table 3.1.

V P t-Ar (r) = V 0 e -α(r-r 0 ) - C 6 r 6 , (3.7)
In Figure 3.2.2, the circle represents the Ar atom adsorbed on the fcc-hollow site, numbers 1 represent the Pt atoms of the rst layer, numbers 2 the Pt atoms of the second layer, and numbers 3 the Pt atoms of the third layer. By considering the Ar atom at the vertical position above a fcc-hollow site of the Pt(111) surface, the system presents a third order rotation axis. This rotational symmetry is conserved for the construction of the Pt cluster used for the decomposition. The height of this cluster is 4.57 Å and its width is larger than 8.0 Å. When the Ar atom is close to its equilibrium position, the distance between the Ar atom and the more external Pt atoms of the cluster is more than 8.0 Å. For most of the pairwise potentials studied in this work, including the present potential, the pairwise interaction can be neglected for r P t-Ar > 8.0 Å. The convergence of the recomposed potential V P t(111)-Ar (z) relative to the number of Pt atoms has been checked, in particular the addition of a fourth layer has a negligible contribution in the sum of equation (3.6). Moreover, V P t-Ar (r) is expected to be adequate to reproduce the attractive part of the potential: the C 6 value of 64.92 eV.Å 6 obtained is consistent with the value of 68.15 eV.Å 6 chosen by Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF] to reproduce the C vdW value of Smith et al. [START_REF] Smith | A molecular dynamics study for the trapping and scattering of Ar/Pt(111)[END_REF] calculated from the polarizability of the Ar atom and the dielectric function of Pt.

Computation of TMAC by Molecular dynamics

In this section, we describe the MD method to simulate the gas-wall interaction and directly calculate the accommodation coecient. Each gas-wall collision is treated separately as in most direct scattering experiments (see [START_REF] Hurst | Direct inelastic scattering Ar from Pt(111)[END_REF][START_REF] Eichenauer | Interaction potential for one-phonon inelastic hecu (111) and heag (111) scattering[END_REF]). This method is suitable for dilute gas for which it can be assumed that only one atom interacts with the wall in one collision event. In the three-dimensional MD simulation, the Ar atoms beams are projected with dierent incident angles on a clean Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface lying on the Oxy plane. The x direction corresponds to the < 110 > direction associated with the maximum Pt atom density direction, such that the distance between two atoms equals 2.772 Å. The directional σ dir t coecient associated with each θ and φ is dened by the following formula [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF]:

σ dir t (θ, φ) = ⟨v in ⟩ -⟨v rn ⟩ ⟨v in ⟩ , (3.8)
where θ and ϕ are the zenithal and azimuthal angles, respectively. θ is the incident angle of the gas atom with the z axis, n is direction of the projection of the incident and reected velocities on the Oxy plane, making an angle ϕ with the x direction. v in and v rn are the projections of the incident and the reected velocities on n, respectively.

v in = v i sin(θ), v rn = v r sin(θ).
(3.9)

In order to model the Pt-wall and to maintain its temperature, a phantom wall technique [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF][START_REF] Yi | Molecular dynamics simulation of vaporization of an ultra-thin liquid argon layer on a surface[END_REF] is used in this work. Each Pt atom is connected to the nearby ones by harmonic springs with rigidity k w , allowing it to vibrate around its site. The wall structure is made of 924 atoms, arranged in (111) fcc four layers slab of size 40 × 40 × 6.77 Å 3 (see Figure 3.3). The last two layers are the phantom layers, the lowest is xed. The temperature of the bulk part of the wall, T w , is controlled by a Langevin thermostat applied to the second phantom layer [START_REF] Schlick | Molecular Modeling and Simulation: an Interdisciplinary Guide[END_REF]. The motion of an upper phantom atom i is governed by the equation:

m i dv i (t) dt = -ξv i (t) + f i (t) + R i (t) .
(3.10)

In equation (4.6), v i is the velocity of the atom i, f i is the resulting force acting on it by the surrounding ones, m i is the atomic mass and ξ is the damping coecient. The third term in the right hand side of equation (4.6) is then the random force applied on the atom.

In the simulation, it is sampled after every time step δt from a Gaussian distribution with zero average and mean deviation of √ 6ξk B T w /δt. The simulations were carried out by setting rigidity, time step and damping factor at the following values [START_REF] Yi | Molecular dynamics simulation of vaporization of an ultra-thin liquid argon layer on a surface[END_REF]: (3.11)

k w = 46.
The spring constant k w has been optimized by Maruyama and Kimura to reproduce the physical properties of solid platinum [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF]. The phantom atoms technique allows to model the semi-innitely large bulk solid kept at a constant temperature T w with the proper heat conduction characteristics. We also checked that to add on two layers does not change the σ t values. The gas velocity at each time step is calculated by the usual Leap-Frog Verlet integration scheme. For a given pair potential V (R ij ), we can compute the force vector f ij that atom j at position R j exerts on atom i at position R i by:

f ij = - dV P t-Ar dr r ij ∥r ij ∥ , r ij = r i -r j .
(3.12)

The computational costs may be reduced by assuming that an Ar atom only interacts with the Pt wall within a cuto distance r c =12 Å, corresponding to a negligible attraction of about 0.02 meV. On the other hand, r c is very small with respect to the mean free path. Initially, Ar atoms are given with a xed incident velocity v i dened in the Oxz plane and with a norm associated to the gas beam temperature T g , such that

v iτ = √ ⟨v 2 iτ ⟩ = √ k B T g m Ar .
(3.13)

This value corresponds to the average value of the positive part of the incident velocity along the τ direction. The wall temperature T w is kept at 200 K, 300 K and 400 K and the gas beam temperature T g is kept at a slightly higher value than T w , here T g = 1.1 × T w .

The Ar atoms are inserted randomly at the distance r c with respect to the wall surface and their reected velocity v r is recorded when they bounce back beyond this cuto distance. Whenever an Ar atom crosses this boundary, another one will be inserted with velocity v i .

The process continues until we obtain a converged value of ⟨v rn ⟩, i.e. after approximatively 12000 collisions. A typical run requires 4 × 10 7 time steps of 2 fs. The periodic boundary conditions are applied along the x, y directions throughout the simulation.

3.4

Results and discussions

Interaction potentials

The Pt(111)-Ar interaction potential and the corresponding Pt-Ar pairwise potential calculated in the present work are plotted in Figures 3.3 and 3.4, respectively. The potentials described in the above referenced works (see Section 2.1) are also shown. Table 3.1 contains the tted parameters and Table 3.2 presents the equilibrium interaction energies, V e , the equilibrium value, z e , of the Pt(111)-Ar distance and the vertical Pt(111)-Ar harmonic vibrational frequency, ω e . These data are compared with those deduced from other potentials and with available experimental values (V e = -80 meV and ω e = 5 meV [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF][START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF]). ω e has been estimated from the second derivative of the interaction potentials at their respective minimum geometries, using the Ar mass as the reduced mass considering that the Pt surface has an innite mass. [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF] and by Spijker et al. [START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF]. They depict a relatively weak value of V e , of ≈ -40 meV, and a short equilibrium z e < 3 Å. The second group of potentials includes those of Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] and Svanberg and Pettersson [START_REF] Svanberg | Survival of noble gas clusters scattering from hot metal surfaces[END_REF]. They reproduce correctly the experimental value of V e , but with relatively large values of z e compared with the other potentials. Additionally, we can mention the potential of Yamamoto [START_REF] Yamamoto | Slip ow over a smooth platinum surface[END_REF], based on the parameters derived by Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF], but with 2σ = 3.2 Å -2 instead of 2σ = 1.6 Å -2 . The resulting potential depicts a weak value for V e and a relatively large equilibrium distance. The third group concerns the potentials deduced by Smith et al. [START_REF] Smith | A molecular dynamics study for the trapping and scattering of Ar/Pt(111)[END_REF] and Lahaye et al. [START_REF] Lahaye | Site dependent energy loss in Ar scattering from Pt(111)[END_REF]. Both use global van der Waals attractive part of the potential between Ar and the Pt surface. They cannot be decomposed into a pairwise sum. The corresponding V e are lower than -100 meV and the z e are ≈ 3 Å.

The last group includes the present potential, and those of Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF] and Ramseyer et al. [START_REF] Han | Thermally-driven nanoscale pump by molecular dynamics simulation[END_REF]. Our V e value well reproduces the experimental value of approximately -80 meV and is ranged between the Kulginov and Ramseyer ones. The z e = 3.35 Å deduced from our potential is also consistent with those of about 3.25 Å from the Kulginov or Ramseyer potentials. Experimentally, no information exists about the equilibrium distance between the Ar atom and the Pt surface but, such values around 3.2 and 3.4 Å should be correct since the last group of potentials reproduce properly as well the experimental equilibrium interaction energy as the experimental vertical Pt(111)-Ar harmonic vibrational frequency. The values of ω e of this group of potentials are indeed very close to the experimental value of 5 meV.

Finally, the forces corresponding to all considered pairwise potentials are plotted in Figure 3.5. One can notice that the |f (R P t-Ar )| functions (equation (3.12)) decay rapidly as R P t-Ar increases and can be negligible when R P t-Ar ≥ 3R e .

It must be noticed that, although the present and the Kulginov interaction potentials use the same mathematical expression, the repulsive part parameters are strongly dierent. In the Kulginov potential, the V 0 and α parameters have been determined from measured scattering data. This is not possible for R 0 which accounts for the position of the repul- In our model, the potential can be seen as a Morse potential for which the attractive part is replaced by a van der Waals function. R 0 is seen as the minimum position of the potential where the repulsive part equilibrates the repulsive one. The dierences between the present and the Kulginov parameters could come from the fact that the Kulginov parameters, V 0 and α, have been determined before R 0 , whereas in our case, R 0 is almost xed before optimizing V 0 and α. Both potentials are very similar for energies within the minimum and the dissociative region, as shown in Figures 3.3 and 3.4. However, strong divergences appear for very high energies, which are out of the scope of this study.

TMAC results

In order to consider the eect of the wall model, two models have been considered in a preliminary step. The rst model simulates the mutual interaction between a gas atom and the wall (Model 1), whereas in the second model (Model 2), the gas atom exerts no force on the wall atoms. The former is more realistic but very expensive and for a given set of wall atoms, only one gas-wall collision is considered in one simulation. In Table 3.3, the σ dir t (45 • , 0 • ) values have been computed from Models 1 and 2 using our new pairwise potential and the one used by Yamamoto [START_REF] Yamamoto | Slip ow over a smooth platinum surface[END_REF]. For both potentials, σ dir t (45 • , 0 • ) = 0.14 if the simplied Model 2 is considered. , 0 • ) values. Model 2 corresponds indeed to a much harder surface and, the reection tends to be more specular than with Model 1. Numerical results based on Model 1 show clearly that the wall can be locally deformed due to the collision. The gas atoms stay longer near the surface and interact deeply with it, resulting in a more diusive reection. This preliminary study conducted us to use Model 1 for all of the results discussed in what follows.

The values of σ dir t (45 • , 0 • ) computed at 300 K with an incident angle of 45 • are tabulated in Table 3.4. The main characteristics of the pairwise potentials (i.e. the dissociation energies, D e and the equilibrium positions, R e ) are also given. The σ dir t (45 • , 0 • ) values reported in Table 3.4 vary between 0.23 and 0.49. Such dierences are closely related to the type of potential used. A rst group of pairwise potentials, that includes the Maruyama and Kimura [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF] and Spijker et al. [START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF] potentials, present low dissociation energies D e (≈ 6 -7 meV) and short equilibrium positions (< 3.5 Å). For both of these potentials, the

σ dir t (45 • , 0 • ) is about 0.3.
It should be noticed that these Lennard-Jones type potentials are commonly used in MD simulations [START_REF] Cao | Temperature dependence of the tangential momentum accommodation coecient for gases[END_REF] even if they do not reproduce properly the Pt(111)-Ar interaction (see Section 2.2). The second group results from interaction potentials associated with V e close to the experimental value and large z e , i.e. Head-Gordon and Svanberg potentials [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF][START_REF] Svanberg | Survival of noble gas clusters scattering from hot metal surfaces[END_REF]. These pairwise terms are characterized by relatively large values of R e (> 4.5 Å) resulting in small σ dir t (45 • , 0 • ) values, close to 0.25. The third group deals with the pairwise potential of the present work together with the Kulginov and Ramseyer ones [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF][START_REF] Han | Thermally-driven nanoscale pump by molecular dynamics simulation[END_REF]. All these two bodies terms lead to σ dir t (45 • , 0 • ) values in the 0.4 -0.5 range. They present equilibrium geometries around 3.8 -3.9 Å and D e ranging from 9.1 meV to 10.2 meV.

If σ t ≈ 1, the wall is almost diusive whereas if σ t ≈ 0, the wall is almost specular. The wall tends to become diusive if the Ar atoms can approach closer to the wall and interact more with it. This situation is obtain for the third class of pairwise potentials whereas the pairwise potentials of the rst and second classes lead to more specular interaction of Ar with the wall. In the rst class, the pairwise terms present shorter equilibrium , ϕ) has been evaluated for dierent values of ϕ and the resulting values are depicted in Figure 3.6. The orientation of the Oxz plane with respect to the surface has a negligible impact on the σ t value and the surface (111) can be seen as isotropic. Table 3.5 shows that the σ dir t (θ, 0 • ) variations with the incident angle are within 0.40 to 0.45 for the new potential suggested in present work. When averaging over θ, i.e. summing the dierent σ dir t (θ, 0 • ) values determined for the values of θ displayed in Table 3.5 and dividing by the number of the θ values, the resulting σ t is 0.42. This averaged value could be associated to the σ t parameter in the Maxwell's model. Figure 3.4.2 gives the full distribution of the x component of the reected velocity for 12000 collisions using the present pair potential for T w = 300 K, θ = 45 • and ϕ = 0 • . The standard deviation of the tted curve equals 155.4 m.s -1 and is related to the TMAC standard deviation: We notice that the deviations reported in Table 3.5 are very high, especially for small incident angles. These high deviations can be explained through analyzes of the post-collision path of the Ar atoms. After the collision with the wall, an Ar atom can bounce backward or forward with a large reection angle. Therefore, the accommodation coecient for a single Ar atom can be very high and, even negative. Theoretically, the σ t corresponding to a single collision can take any negative or positive value according to the Maxwell model. This phenomenon is in agreement with the work of Finger et al. [START_REF] Finger | Molecular dynamics simulation of adsorbent layer eect on tangential momentum accommodation coecient[END_REF]. The deviation is higher with small θ since the initial tangential velocity in the denominator of equation (3.14) decreases.

σ T M AC = σ ⟨vrx⟩ ⟨v ix ⟩ . ( 3
Although most of the results of this work are obtained using a constant incident velocity corresponding to the gas temperature, we have done separate simulations using the Maxwell-Boltzmann distribution for the normal incident velocity:

f ′ (v iτ ) = ( m g 2πk B T g ) 1 2 × e - mg v 2 iτ 2k B Tg .
(3.15)

The resulting value of σ dir t (45 • , 0 • ) value is 0.39 using the present potential for T w = 300 K. This result is almost identical to the σ dir t value of 0.40 in Table 3.4 obtained by the use of the constant v iτ dened by equation (3.13).

We now compare some previous MD studies on the σ t determination for the Pt-Ar couple with the present model. In their work, Sun and Li [START_REF] Sun | Three-dimensional molecular dynamic study on accommodation coecients in rough nanochannels[END_REF] have simulated Couette ows in smooth and rough nanochannels and have varied the temperature and the types of the surface crystal.

In their model, they have used Pt-Ar pairwise potential derived by Maruyama and Kimura [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF] and the Einstein model for the Pt slab vibration. For an isothermal ow at 200 K, they have obtained a σ t value of 0.348 which can be compared with σ dir t (45 , 0 • ) = 0.47, larger than the Sun and Li one. The low value derived by Sun and Li is probably due to their choice for the wall model since they have chosen a simple spring description based on the Einstein model leading to k = 179.5 N.m -1 , higher than the present k wall value. This implies a more specular behavior of the collision. Spijker et al. [START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF] have explored the computation of accommodation coecients by the use of velocity correlation proles determined from MD simulations. Contrary to the above studies, the Argon gas conned between two Pt walls was not subject to a ow but was at equilibrium. The Pt-Pt interactions in the walls were described by a Lennard-Jones potential. They have reoptimized the parameters of the Lennard-Jones Pt-Ar derived by Maruyama et al. [START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF] to correctly describe the wettability character of Ar on a Pt surface, and at 300 K they derived a σ t value of 0.53 whereas using their potential we obtained σ dir t (45 • ) = 0.28. The dierence of the values can be due to either the wall model or the way of computing σ t . Their model is based on velocities correlations between incoming and outcoming particles colliding with the walls to compute σ t from in equilibrium gas atoms.

Very recently, using the Pt-Ar potential reoptimized by Spijker et al. [START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF], Prabha and Sathian [START_REF] Prabha | Molecular-dynamics study of poiseuille ow in a nanochannel and calculation of energy and momentum accommodation coecients[END_REF] have simulated a Poiseuille ow in a nanochannel by MD. The Ar gas nonequilibrium ow was maintained due to a gravity eld. For a wall temperature of 300 K, a Knudsen number of 0.74 and a gravity eld of 5.10 -7 nm/ps 2 , they obtained a σ t value of 0.290.This result agrees within 0.01 with the present σ dir t (45 • , 0 • ) value of 0.28 computed at the same temperature and using the Spijker potential. This agreement allows to validate the present MD model.

Agrawal and Prabhu compiled theoretical and experimental σ t values in a survey. The averaged value is of 0.983 [START_REF] Prabhu | Survey on measurement of tangential momentum accommodation coecient[END_REF] for Ar atoms independently of the surface material. Experimentally, Cooper et al. [START_REF] Cooper | Gas transport characteristics through a carbon nanotubule[END_REF] deduced a value of 0.52 ± 0.1 from the rst experimental determination of the slip coecient for Ar through tubular carbon structures, and Arkilic et al. [START_REF] Arkilic | Mass ow and tangential momentum accommodation in silicon micromachined channels[END_REF] obtained σ t between 0.75 and 0.85 for measurement of Ar mass ows through silicon micromachined channels due to an imposed pressure gradient. The σ t dierence between these experimental works gives an idea of the importance of the material surface.

To the best of our knowledge, no experiment exists on the σ t determination for the Pt-Ar couple making a good comparison dicult.

In order to conclude this study about the eect of the pair potential on the σ t determination, the results of a scattering simulation are presented. Head-Gordon et al. also reported molecular beam experiments for the Pt(111)-Ar system. Figure 2 of Ref. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] shows the velocity distribution of the scattered ux of Ar atoms with an incident energy of 2.5 kJ.mol -1 and an incident angle of 45 • on a Pt [START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface at 500 K. We have simulated this experiment using the same parameters for the atomic beam with T g = 600 K, corresponding to an energy of 2.5 kJ.mol -1 , and the present potential. The reected velocity has been plotted as a function of the reected angle and the corresponding curve is compared with the experimental and simulated data of Head-Gordon et al. in Figure 3.8.

The simulated curve of Head-Gordon et al. was based on the use of the corresponding potential of Table 3.1, which has been optimized on the measured data of Figure 3.8. The present curve presents the same main trends as the experimental points. The reected velocity is maximum for θ r around 20 • and for θ r larger than 40-50 • , the norm of v r is almost constant. However, the present calculated velocities are higher by approximately 100 m.s -1 in comparison with the measured data. ) recorded for dierent values of the reected angle θ r (in degrees) with respect to the normal direction of the surface. The present pair potential has been used as well as T g = 600 K, T w = 500 K, θ = 45 • and ϕ = 0 • . The dotted curve corresponds to the experimental values recorded by Head-Gordon et al. [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] and the markers to the simulated data from the same work.

+D=FJAH "

Eects of wall surface

Introduction

As mentioned earlier, based on Eq. (5.45) of Maxwell's model, the σ t parameter for a gaswall couple can be determined by either experiments [START_REF] Arkilic | Mass ow and tangential momentum accommodation in silicon micromachined channels[END_REF] or Molecular Dynamics [START_REF] To | A slip model for micro/nano gas ows induced by body forces[END_REF][START_REF] Cao | Temperature dependence of the tangential momentum accommodation coecient for gases[END_REF] in the Navier Stokes slip regime. However, most MD simulations of ows were done at nanoscale [START_REF] To | Molecular dynamics simulations of pressure-driven ows and comparison with acceleration-driven ows[END_REF] and did not have the same conditions as in experiments. In order to compare σ t calculations with measurements for dilute gases, a more relevant MD approach [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF] consists in studying every single gas-wall collision event. Consequently, σ t can be computed directly by projecting gas atoms into the surfaces and nding momentum changes [START_REF] Finger | Molecular dynamics simulation of adsorbent layer eect on tangential momentum accommodation coecient[END_REF]. This approach, which is quite similar to beam experiments [START_REF] Rettner | Thermal and tangential-momentum accommodation coecients for N 2 colliding with surfaces of relevance to disk-drive air bearings derived from molecular beam scattering[END_REF], provides insights into the reection mechanism and can be used to improve Maxwell's model. As far as multi-scale simulations are concerned, the obtained uid-wall interaction results can be coupled with other numerical methods [START_REF] Struchtrup | Macroscopic Transport Equations for Rareed Gas Flows: Approximation Methods in Kinetic Theory[END_REF][START_REF] Liou | Microuid Mechanics[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]9799].

In Eq. (3.1), the use of one parameter σ t as in Maxwell's model means that the slip behavior is isotropic. For anisotropic textured surfaces, more sophisticate models are needed to reproduce the direction dependent slip or gas-wall interaction behavior. Bazant and Vinogradova [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF] suggested using a slip length tensor to quantify this behavior. The tensorial nature of the slip eect was shown to be related to the interfacial diusion [100? 102]. Eective slip tensors with bounds for ows over superhydrophobic surfaces were also obtained [START_REF] Feuillebois | Eective slip over superhydrophobic surfaces in thin channels[END_REF][START_REF] Feuillebois | Transverse ow in thin superhydrophobic channels[END_REF]. As the slip models describe macroscopic behaviors, it is thus relevant to investigate the problem at the scale of uid wall interaction. For gases, Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] generalized Maxwell's scattering kernel by using anisotropic accommodation coecients. The consequences of the model on the slippage have not been studied. Since the anisotropic scattering kernel model does not provide full information about the gas-wall collisions, we shall use MD method to study these interactions in detail with the focus on the surface morphology. The MD code used in this chapter is the parallel version described in Ref. [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. The original code has been enriched (e.g multi-body potentials, statistical tools, etc..) to adapt to the aim of the present work.

Generally, results obtained from MD simulations depend on the following factors:

-The interaction potential between the gas/wall atoms as mentioned in Chapter 3.

-The dimension of the simulation models. In general, 3D models are better than 2D since it accounts for interactions of the gas atom with all its neighbors.

-The potential between the solid atoms must be good enough to reproduce the free surface eect. It is well known that the distance between the atomic layers near the free surface are much smaller than in the bulk.

-The temperature eect must be considered as gas molecules are adsorbed easier at cold walls than at hot walls, which can results in a higher σ t .

-The surfaces are not always ideally smooth and can have dierent morphology (e.g randomly rough or textured surfaces).

This chapter aims at including these features in simulations of molecular beam experiments. The gas/wall couple under consideration is Argon and Platinum but the methodology of the present work can be used to obtain σ t for any gas/wall couple provided that an appropriate potential is used. The chapter is organized as follows. After the Introduction, Section 2 is devoted to the description of the computational method. It discusses briey the choice of potentials, the method to prepare surface samples and MD simulation of gas/wall interaction. We remark that a part of surface sample preparation requires a separate molecular dynamics simulation of lm deposition processes in order to create a realistic random roughness surface. The σ t results issued from the calculations are then shown in Section 4. Finally, conclusions and perspectives are discussed.

Computation model 4.2.1 Interatomic potential

The interatomic potentials play an important part in the molecular dynamics simulations since they govern the dynamics of the system and, thus the accuracy of the results. In this work, the following van der Waals type pair potential between At and Pt derived by Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF] is used ϕ Ar-P t (r Ar-P t ) = V 0 e -α(r Ar-P t -r 0 ) -C 6 r 6

Ar-P t , r Ar-P t = |r Ar -

r P t | , (4.1)
where r Ar-P t is the distance between an Ar atom at location r Ar and a Pt atom at location r P t , or vice versa. As mentioned in Chapter 1, contrary to usual Lennard-Jones potentials, the repulsive part of this pair potential has a Born-Mayer form and provides a better description of the strong repulsion of the electrons. The pairwise potential parameters have been empirically adjusted such that the laterally average potential reproduces the measured properties of an Ar atom adsorbed on a slab of Pt atoms, i.e. a well depth of about 80 meV [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] and, a vibrational frequency of the adsorbed atom of about 5 meV [START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF]. The van der Waals interaction of an Ar atom with a platinum surface can be evaluated from the Ar-polarizability and the Pt-dielectric function. The values of the potential parameters are given in Table 4.1 and were shown to be in good agreement with an ab-initio based calculation [START_REF] Léonard | Inuence of the pairwise potential on the tangential momentum accommodation coecient: a multiscale study applied to the argon on Pt(111) system[END_REF].

Table 4.1: Parameters of the Pt-Ar pairwise potential [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF].

V 0 (eV) α (Å -1 ) R 0 (Å) C 6 (eV.Å 6 ) 20000 3.3 -0.75 68.
15

In terms of the potential between the Pt atoms, the multi-body Quantum-Sutton Chen (QSC) potential is used [START_REF] Sutton | Long-range FinnisSinclair potentials[END_REF]. As a particular Finnis-Sinclair potential type, the QSC potential includes quantum corrections and predict better temperature dependent properties.

Table 4.2: Quantum Sutton-Chen parameters for Pt [START_REF] Sankaranarayanan | Molecular dynamics simulation study of the melting and structural evolution of bimetallic Pd-Pt nanowires[END_REF].

n m ϵ (eV) c a (Å) 11 7 9.7894e-3 71.336 3.9163

For a system of N Pt-atoms, the potential is given by the following expression

V pot,P t = ϵ     1 2 N ∑ i=1 N ∑ j=1 j̸ =i ( a R ij ) n -c N ∑ i=1 ρ 1/2 i     , ρ i = N ∑ j=1 j̸ =i ( a R ij ) m , (4.2)
where a is the lattice constant, R ij the distance between atom i and j and the local density ρ i of atom i. The parameters ϵ and a determine the scales of energy and length, respectively and, n and m the range and shape of the potential. These potential parameters are given in Table 4.2. Combining the Ar-Pt and Pt potentials, we can compute the total potential of the system

V pot = N ∑ i=1 ϕ Ar-P t (R Ar-i ) + V pot,P t , (4.3)
and the force f i acting on atom i at position r i by

f i = - ∂V pot ∂r i . (4.4)
Since we only consider the interaction of one Ar atom with a Pt surface, there is no contribution of the Ar-Ar term in the total potential formula V pot . The accuracy of the QSC potential for Pt has been justied in Ref. [START_REF] Sankaranarayanan | Molecular dynamics simulation study of the melting and structural evolution of bimetallic Pd-Pt nanowires[END_REF] as it reproduces accurately the melting temperature and the specic heat of the material. Although its implementation is more costly than the harmonic (spring) potential, it should better reproduce the surface eects, since atoms near the free surfaces are dierent from the bulk. Our tests on the QSC potential show that in fully relaxed equilibrium system, the interatomic distance near the free surfaces is much smaller than in the bulk (see Fig. In this work, three types of surfaces are considered: smooth surfaces, periodic nanotextured surfaces and randomly rough surfaces. The orientation of their free surfaces is (100) according to the Miller index. Initially, the Pt atoms are arranged in layers and the two lowest ones (phantom atoms) are used to x the system and for the thermostat purpose. The remaining Pt atoms are free to interact with other solid atoms and gas atoms. The random arrangement of these atoms denes the "rough" state of the surface and will be detailed later on.

A smooth surface model is a system composed of 768 atoms arranged in 6 layers, all of which are in perfect crystal order. The nanotextured models are constructed from the smooth surface model by adding successively atom layers to create pyramids with the slope angle 45 • . The slope is necessary to assure the stability of the system since perfectly vertical blocks (slope angle 90 • ) are less stable: in many cases atoms migrate to lower positions and the blocks evolve into step-like structures with smaller potential energy. The base of the pyramid can be a square (type A, Fig. 2) or an innite strip (type B, Fig. 3), so that both isotropic/anisotropic eect can be considered. Although these pyramids are simplied models of a real rough surface, it can show the dependence of σ t on the roughness. The latter in MEMS/NEMS is reported to be several Å [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF]. In this work, the highest peak, varying with the number of atoms layers added on the surfaces, ranges from 2 to 6 Å.

Randomly rough surface models are also constructed by adding atoms on the smooth surfaces in a random way. In the available literature, there are several mathematical mod- els [112115] that describe random roughness. However, these models are not suitable at atomic scale: it is dicult to force atoms to be at given positions and, structure parameters such as orientation (100) and lattice constant must be respected. Furthermore such atomistic systems might not be appropriate in terms of potential energy. In our opinion, a randomly rough surface which is consistent with the internal atomistic structure, should be built from molecular dynamics simulations. Rapid cooling of thin lms from the liquid state [START_REF] Hoang | Free surface eects on thermodynamics and glass formation in simple monatomic supercooled liquids[END_REF] can create rough surfaces but the nal systems could contain many defects (e.g. pores, dislocations) and non crystalline structure (e.g. polycrystal, amorphous). As the present work focuses on Pt (100), the rough surfaces are constructed by deposing atoms randomly on the existing smooth Platinum surface. Since this procedure is quite similar to the lm vapor-deposition processes, it is asssumed that the created surface is quite close to real MEMS/NEMS surfaces. The procedure of the material deposition is described as follows.

The initial system is a Pt plate made of four layers of 512 solid atoms, arranged in (100) fcc order. First, the system is relaxed towards the minimal potential energy conguration. Then, after 2000 time steps of 1 fs, a Pt atom is inserted randomly from a height of 10 Å with the initial thermal velocity corresponding to 1000 K. Under the attraction force (QSC potential) from the Pt plate, the deposed Pt atoms move downwards until they reach the plate which is maintained at 50 K. (see a snapshot of the deposition process in Fig. 4.4). Finally, when all inserted Pt atoms are attached rmly into the Pt plate, the whole system undergoes the anneal process at the ambient temperature T a = 300 K with a time step equal to 2 fs. During the whole simulation, the Leap-Frog Verlet integration scheme is employed and the temperature is kept constant by simple velocity scaling method. In what follows, we describe the MD method used to simulate the gas-wall collision and to calculate the σ t coecient. The simulations are three dimensional: an Ar atom is projected into a Pt(100) surface with dierent incident angles θ and with dierent approaching φ-planes. In spherical coordinate system, (θ, φ) are the polar and azimuthal angles, respectively (see Fig 4.6). As mentioned in Chapter 2, the coecient σ dir t (θ, φ) is given by

σ dir t (θ, φ) = ⟨v in ⟩ -⟨v rn ⟩ ⟨v in ⟩ , (4.5)
where v in and v rn are the projections of the incident and the reected velocity on the vector n, respectively. The latter is the intersection of the plane xOy and the φ-plane, i.e it lies on xOy and makes an angle φ with respect to Ox. Only one gas-wall collision is treated per simulation and the averages ⟨v in ⟩, ⟨v rn ⟩ in Eq. (4.5) are taken over a large number of simulations (or collisions). The denition (4.5) is the most accurate description of gas-wall interaction since it is associated to each direction. We also calculate the eective anisotropic σ an t (φ) coecients using the same equation (4.5) but with gas atoms arriving from all directions: the direction of v i is randomly uniform with v in > 0. In the special case where the surface is isotropic, σ an t varies little with φ and a single eective isotropic σ iso t constant is sucient for modelling gas-wall interaction as in Maxwell's model. The latter eective isotropic coecient is obtained by the similar method but ⟨v in ⟩, ⟨v rn ⟩ in Eq. (4.5) are further averaged over n (or φ).

We assume rst that an Ar atom only interacts with the Pt wall within a cuto distance r c = 10 Å. Since this distance is much smaller than the typical mean free path at atmospheric pressure or in high vacuum (λ > 69 nm), it can justify the choice of such a small region to calculate the σ t coecients. At a distance of 10 Å, an Ar atom can be considered as non-interacting with the Pt-wall atoms since the potential value at that distance (-0.0580736meV) is negligibly small when compared with the potential well depth (10.21 meV). At the beginning of each simulation, an Ar atom is inserted randomly at the height r c above the wall surface with initial incident velocity v i . The norm of v i is equal to the thermal speed corresponding to the gas beam temperature T g . Although the results of this work are obtained using a constant incident velocity corresponding to the gas temperature, we have done separate simulations using the Maxwell-Boltzmann velocity distribution and nd that σ t is insensitive to this modication. A collision is considered as nished when the atom bounces back beyond the cuto distance. Then the reected velocity v r is recorded for the statistical purpose and another Ar atom is reinserted randomly to continue the process. After approximatively 10000 collisions (simulations), converged values of σ t values were obtained. Numerical tests show that the statistical error of a typical 10000-collision average is within 1.0%.

Throughout the simulations, periodic boundary conditions were applied along the x, y directions. The velocities and positions of gas atoms and the solid atoms at each time step are calculated by the usual Leap-Frog Verlet integration scheme. To control the temperature T w of the system, the phantom technique is used: the Langevin thermostat [START_REF] Schlick | Molecular Modeling and Simulation: an Interdisciplinary Guide[END_REF] is applied to the atom layer above the xed layers. The motion of an atom i belonging to this layer is governed by the equation

m i dv i (t) dt = -ξv i (t) + f i (t) + R i (t) . (4.6)
In Eq. (4.6), v i is the velocity of the atom i, f i is the resulting force acting on it by the surrounding ones, m i is the atomic mass and ξ is the damping coecient. The third term R i in the right hand side of Eq. (4.6) is the random force applied on the atom. In the simulation, it is sampled after every time step δt from a Gaussian distribution with zero average and mean deviation of

√ 6ξk B T w /δt.
The simulations were carried out by setting time step and damping factor at the following values:

δt = 2 fs , ξ = 5.184 × 10 -12 kg/s . ( 4.7)
The wall temperature T w was kept at 200 K, 300 K and 400 K and the gas beam temperature T g was kept at a slightly higher value than T w , here T g = 1.1T w . Such choice of T g was made arbitrarily and the procedure of the present work can be applied to any gas temperature. Generally, to obtain the best statistical results, a typical run requires 4 × 10 7 time steps of 2 fs. All simulations were run on 9 processors, using a domain decomposition and the Message Passing Interface. The longest simulation takes about 20 CPU hours. We have carried out computations with dierent time steps from 1fs to 3fs and we have found that the results are insensitive to this factor.

Eects of temperature and roughness height

From the description of the models in Section 2, the coecient σ dir t can depend on the several input parameters: temperature, surface morphology, incident direction (θ, φ). The variation of σ dir t in terms of these parameters is investigated in the following subsections.

The σ dir t results at dierent temperatures are shown in Tables 4. [START_REF] Struchtrup | Macroscopic Transport Equations for Rareed Gas Flows: Approximation Methods in Kinetic Theory[END_REF] and Fig. 4.8. A general trend can be noticed here: σ dir t increases as the temperature decreases, ranging from 0.78 to 0.92 in the case of the highest roughness considered (h = 5.88 Å). This trend in σ dir t variation can be explained by the fact that the adsorption is stronger with colder walls. Gas atoms stay longer near the wall, interact more with solid atoms and, as a result, the reection is more diusive. Similar remarks have been reported in Refs. [START_REF] Cao | Temperature dependence of the tangential momentum accommodation coecient for gases[END_REF][START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF] for conned systems. For h = 5.88 Å and T w = 300 K, Table 4.4 shows that the σ dir t value varies very little with the incident angle θ and very close to the average isotropic value σ iso t = 0.85. This means that for this kind of surface, Maxwell's one parameter model is suciently accurate to model gas-wall interaction. 4.8). Computations carried out for pyramidal structures at the temperature of 300 K show that the σ dir t coecient can reach up to 0.87 for surfaces with the highest peak conguration. It is clear that the presence of peaks leads to non uniform surface potentials with local minima where gas molecules can easily be trapped: the gas atoms adhere stronger to the wall and lose their initial momentum. Moreover, the changes in local slopes produce more or less random variations in the local incident and reection angles.

Visualization of collision trajectories shows a clear dierence between a smooth surface and a rough surface. On a smooth surface, a gas molecule collides and bounces several times before nally escaping from the inuence distance r c of the wall (see Fig. 4.9). On a rough surface, it stays near the wall and moves within the valley between the peaks, a mechanism similar to surface diusion, until the wall provides enough energy to escape (see Fig. 4.4). Thus, in addition to the roughness height, the in-plane random arrangement of the atoms also plays a signicant role on the accommodation coecient.

Surface anisotropy eect

An anisotropic textured surface can obstruct or facilitate the ows dierently along dierent directions. Bazant and Vinogradova [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF] generalized Navier slip boundary conditions for anisotropic texture surface by using a tensorial slip length. In the framework of the kinetic theory, Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] proposed a new scattering kernel that accounts for surface anisotropy. Their formulation is based on three independent accommodation coecients α x , α y , α z along the three directions x, y, z. The coecients α x , α y represent the tangential accommodation coecients and α z is the normal accommodation coecient. The tangential accommodation coecient α n in direction n is then computed by the expression (see Appendix I)

σ an t (φ) = α n = α x cos 2 φ + α y sin 2 φ . (4.8)
We remark that by substituting φ = 0 • and φ = 90 • , the accommodation values α x and α y along x, y direction can be recovered. In this subsection, we study the anisotropy eect using MD and the directional σ t denition in Eq. (4.5) and examine the relation (4.8).

The anisotropy eect can be seen from Figs. 4.11 and 4.12: the σ dir t variation with φ is non uniform for rough surfaces. The accommodation process along the two directions x, y is highly dierent. The σ dir t is minimum when the atoms are projected along the longitudinal direction of the strip (φ = 90 • ), since the surface may be considered as almost smooth in that direction (see Fig. 4.3). This σ t value corresponds to α y in the model of [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF]. The maximal σ t values recorded for φ = 0 • and h > 0 can be attributed to the largest roughness eect in that direction and correspond to α x in the model [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF]. Moreover, Figs. 4.11 and 4.12 show an increase of anisotropy eect as the roughness increases: the dierence between the highest and the smallest σ t value increases with the roughness height whereas the σ t results depend very little on the beaming direction for a smooth surface. This could be explained by the fact that the smooth surface can be considered isotropic. Although Figs. [START_REF] Liou | Microuid Mechanics[END_REF].11 and 4.12 show discrepancies of σ t obtained in dierent ways, all curves can t reasonably well the analytical relation (4.8) as plotted in solid, dashed and dash-dotted lines. For anisotropic surfaces, the reected ux is not always lying in the same plane with the arriving one. Consequently, in addition to Eq. ( 4.5), we should account for the ratio of the reected ux components along two orthogonal directions m, n: ⟨v rm ⟩/⟨v rn ⟩.

According to the anisotropic model (see Appendix I), this ratio can be computed by the expression

⟨v rm ⟩/⟨v rn ⟩ = (α x -α y ) cos φ sin φ 1 -α x cos 2 φ -α y sin 2 φ .
(4.9)

By observing the surface structure, we can deduce that ⟨v rm ⟩/⟨v rn ⟩ must vanish for impinging uxes parallel to the planes of symmetry of the anisotropic surface. That remark is in good agreement with Eq. (4.9) where ⟨v rm ⟩/⟨v rn ⟩ = 0 at φ = 0 • , 90 • . Our MD simulation conrms the remark and also shows that the ratio is nonzero at φ ̸ = 0. From Table 4.5, at φ = 45 • , we nd that the ratio is signicant. It even increases with the roughness height increases, i.e the anisotropic eect is enhanced.

+D=FJAH #

Tensorial slip theory for gas ows

Introduction

Among the topics of micro-nanouidics, slip ows at solid surfaces have been studied actively from dierent aspects including experiments, modelling and computational methods [START_REF] Struchtrup | Macroscopic Transport Equations for Rareed Gas Flows: Approximation Methods in Kinetic Theory[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF][START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF][START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Bocquet | Hydrodynamic boundary conditions, correlation functions, and kubo relations for conned uids[END_REF][START_REF] Bocquet | Flow boundary conditions from nano-to micro-scales[END_REF]. In most problems under consideration, Navier slip boundary conditions are used

u s = L s ∂u ∂n , or u s = M τ . (5.1)
In (5.1) 1 , u s is the slip velocity which is assumed to be proportional to the normal velocity gradient ∂u ∂n at the wall and to a slip length L s . Equation (5.1) 2 is another formulation based on shear stress τ and surface mobility constant M . Considering a Newtonian uid of viscosity η, these two formulations are equivalent as L s = M η and τ = η ∂u ∂n . For gases, the slip length L s is shown to be related to σ t and the mean free path λ according to Maxwell's model [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF][START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF] L s = 2σ t σ t λ .

(5.2)

Micro-nano textured surfaces have drawn a lot of scientic interests recently. If suitably designed, these surfaces can become superhydrophobic and reduce signicantly friction. However, due to the strong surface anisotropy, the application of isotropic equation (5.1) is quite limited. A natural generalization of (5.1) while keeping the relations linear is the following tensorial form [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF] u s = b ∂u ∂n , or u s = Mτ .

(5.3)

In (5.3), u s , u, τ are vectors and M, b are second rank tensors having the same meanings as u s , u, τ, M, L s . Specically, the tensors M and b can capture the anisotropy eect when they are not proportional to the identity tensor I. In the gas regime, it is revealed that the accommodation coecient can be orientation dependent [START_REF] Pham | Eects of surface morphology and anisotropy on the tangential-momentum accommodation coecient between Pt(100) and Ar[END_REF] and its variation with respect to the incident direction agrees reasonably well with the anisotropic gas-wall interaction model [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF]. The latter is formulated using a scattering kernel based on three accommodation parameters. This information suggests that the tensors M and b can be obtained from the accommodation coecients in the same manner as Maxwell's isotropic model.

The present chapter which aims at addressing the aforementioned problem, is organized as follows. After the introduction, we shall revisit the concept of anisotropic scattering kernel proposed in Ref [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] and the associated slip models in Section 5.2. Explicit expressions for M and b as functions of the accommodation parameters are presented. Section 5.3 is dedicated to Molecular Dynamics (MD) simulation of Poiseuille ows where the anisotropic gas-wall interaction model is implemented. The MD results will then be compared with the available analytical solutions. Finally, some concluding remarks are given at the end of the chapter.

5.2 Anisotropic slip theory for gas ows

Scattering kernel for anisotropic surfaces

In kinetic theory, the boundary conditions for the distribution function f (v) are usually described using a scattering kernel, for example

v z f (v) = ∫ Ω ′ B(v ′ , v)|v ′ z |f (v ′ )dv ′ .
(5.4)

The scattering kernel B(v ′ , v) is the probability that a gas atom impinging the wall with velocity v ′ is reected with velocity v. In Eq. (5.4), the wall is assumed to be normal to the z direction of the Cartesian coordinate system Oxyz. The velocity components of the incoming and outgoing atoms and the velocity spaces Ω, Ω ′ are dened as follows

v = (v x , v y , v z ) ∈ Ω = {R × R × R + } , v ′ = (v ′ x , v ′ y , v ′ z ) ∈ Ω ′ = {R × R × R -} .
(

Dierent phenomenological models for gas-wall interaction were proposed based on accommodation parameters. The rst popular scattering kernel was derived from Maxwell's postulate on the reection mechanism of the gas atoms [START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF]: it can be either diusive or specular with the respective portion of σ t and 1σ t . Another wellknown model was developed by Cercignani and Lampis [START_REF] Cercignani | Kinetic models for gas-surface interactions[END_REF][START_REF] Cercignani | Rareed gas dynamics: from basic concepts to actual calculations[END_REF] and later extended by Lord [START_REF] Lord | Some extensions to the Cercignani-Lampis gas-surface scattering kernel[END_REF] using two parameters α and α Ez . The former is the accommodation coecient of the tangential momentum and the latter is the accommodation coecient of the kinetic energy due to the normal velocity component v z . Being invariant with respect to any change of basis in the plane xOy, these two models are suitable for describing the behavior of ideally or nearly isotropic surfaces, e.g smooth crystalline surfaces. For periodic nano-textured surfaces (see Fig. 5.1 ), MD simulations of beam scattering experiments in Chapter 4 showed that the coecient σ t depends on the surface orientation. Interestingly, the variation of σ t with respect to the incident azimuthal angle φ agrees quite well with results issued from the anisotropic scattering kernel proposed by Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF].

Let us call Φ - i and Φ + i the incoming and outcoming ux at the wall of the momentum i component dened mathematically by the expression

Φ - i = m g ∫ Ω ′ |v ′ z |v ′ i f -(v ′ )dv ′ , Φ + i = m g ∫ Ω |v z |v i f + (v)dv , (5.6)
where f -and f + are respectively the velocity distribution associated with the incident molecules and reected molecules. Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] proved the following relation

Φ - i -Φ + i Φ - i = α i , i = x, y , Φ - z -Φ + z Φ - z -Φ e z = α z .
(5.7)

In (5.7), the three model parameters α x , α y , α z are the momentum accommodation coecients along x, y, z axes and Φ e z is the reected ux associated to complete accommodation situation: α x = α y = α z = 1.

Since (Φ +

x , Φ + y ) and (Φ - x , Φ - y ) are components of tangential momentum vector uxes Φ + t and Φ - t expressed in the same basis, the rst linear relation in (5.7) can be written in tensorial form

Φ + t = (I -A)Φ - t , (5.8)
where A is the accommodation tensor admitting the diagonal form in the current basis

A = [ α x 0 0 α y ] .
(5.9)

Finally, we note the following relation which results from (5.8)

(Φ + t + Φ - t ) = N(Φ - t -Φ + t ) , N = 2A -1 -I .
(5.10)

In the current basis Oxy, tensor N is reduced to the matrix

N = [ (2 -α x )/α x 0 0 (2 -α y )/α y ] .
(5.11)

By applying the rotation transformation Q(φ)

Q(φ) = [ cos φ sin φ -sin φ cos φ ] ,
(5.12) to (5.8,5.9), we nd the relations between the tangential uxes Φ ± n , Φ ± m in the new basis (n, m) (see Appendix I)

Simplied anisotropic slip model for gas ows based on tangential accommodation coecients

From now, we consider gas ows over a periodic nano-textured surface (see Fig. 5.1) and assume that it is possible to replace the latter with a nominally at surface and an equivalent gas-wall interaction model. It is noteworthy to say that given the boundary conditions for the distribution f , there is no unique derivation of boundary conditions for macroscopic quantities [START_REF] Struchtrup | Macroscopic Transport Equations for Rareed Gas Flows: Approximation Methods in Kinetic Theory[END_REF][START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF][START_REF] Maxwell | On stresses in raried gases arising from inequalities of temperature[END_REF][START_REF] Dadzie | Temperature jump and slip velocity calculations from an anisotropic scattering kernel[END_REF]. Temperature jump and velocity slip equations are usually based on near-wall approximation of f and the assumption that the Navier-Stokes-Fourier (NSF) equations are valid up to the wall. In many cases, the obtained boundary conditions are strongly coupled (i.e stress, heat ux, temperature, velocity, etc.. are all present in the same non linear equation) and dicult to use in practice. In what follows, we adopt the isothermal ow hypothesis and present a simplied slip model for anisotropic surfaces. The derivation of the model relies essentially on the momentum equation and its associated properties (5.7-5.11) while relaxing the other conditions of the anisotropic kernel (A.1 -A.3, Appendix I). The advantage is that it can describe reasonably well the surface anisotropy eect for a broader class of problems and avoid unphysical issues due to the kernels. The shear stress vector τ that the wall exerts on the gas atoms can be computed directly from the dierence of tangential momentum uxes due to collisions

τ = (Φ - t -Φ + t ) , (5.13)
Figure 5.1: Equivalent gas-wall interaction model. In modelling and simulations, the structured surface is replaced by a smooth surface with a scattering kernel R(v ′ , v). In Ref. [START_REF] Pham | Eects of surface morphology and anisotropy on the tangential-momentum accommodation coecient between Pt(100) and Ar[END_REF], the authors considered a (100) crystalline Pt textured surface at 300 K with the following dimensions l 1 = 11.76Å, l 2 = 19.6Å, l 3 = 7.84Å, h = 5.88Å and used MD method to determine the accommodation parameters with Ar atoms associated to model (A.1 -A.3, Appendix I) α x = 0.96, α y = 0.83.

Assuming that the uid near the wall is not too far from equilibrium, we can make the following approximation

∫ Ω ′ f -(v ′ )dv ′ ≃ ∫ Ω f + (v)dv = ρ/2 , (5.14)
where ρ is the total number density of the wall density function f . The latter is dened from f -and f + as follows

f = f -in Ω ′ , f = f + in Ω .
(

Next, we consider that, to some extent, the normal velocity component v z in f + is independent of the tangential velocity

v t = (v x , v y ) and v ′ z in f -is independent of v ′ t = (v ′ x , v ′ y ).
Accounting for (5.14), the tangential momentum vector uxes (5.6) can be estimated by

Φ - t ≃ m g cz ∫ Ω ′ v ′ t f -(v ′ )dv ′ , Φ + t ≃ mc z ∫ Ω ′ v t f + (v)dv , ρ 2 cz = ∫ Ω ′ |v ′ z |f -(v ′ )dv ′ = ∫ Ω ′ |v z |f + (v)dv .
(5.16)

The last equality of (5.16) implies that the wall does not accumulate. Using the denition of gas velocity at the wall

ρu s = ∫ Ω f + v t dv + ∫ Ω ′ f -v ′ t dv ′ , (5.17)
we can deduce that

Φ - t + Φ + t = m g ρc z u s .
(5.18)

Combining (5.10) and (5.18,5.13) leads to the following expression of the slip velocity

u s = Mτ , M = 1 ρm g cz N .
(5.19)

It is interesting to note that the second rank tensor M, is similar to the interfacial mobility tensor in the tensorial slip theory [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF]. In the case of Newtonian uid of viscosity η, we obtain the generalized slip equation for anisotropic surface using the concept of slip length tensor b b = ηM .

(5.20)

We note that the viscosity η is related to the mean thermal speed c and the mean free path

λ [4] η = 1 2 ρm g cλ , c = √ 8kT πm g , ( 5
= λN = λ [ (2 -α x )/α x 0 0 (2 -α y )/α y ] .
(5.23)

It is interesting to remark that (5.23) is consistent with a rather heuristic and less rigorous approach presented in Appendix II. These results are the generalization of earlier works [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF][START_REF] To | A slip model for micro/nano gas ows induced by body forces[END_REF] concerning isotropic surfaces. Another interesting feature of b is that it is independent of temperature and the choice of the kernel form (A.1 -A.3, Appendix I). Finally, it is reduced to Maxwell's slip model when the behavior is isotropic, i.e α x = α y .

Tensor b can be made dimensionless with a characteristic length, for example the channel height (5.25)

In Ref [START_REF] Dadzie | Temperature jump and slip velocity calculations from an anisotropic scattering kernel[END_REF], a slip model associated directly to the kernel form (A.1 -A.3, Appendix I) was proposed. The obtained expression is still rather complex and still depend on the temperature T even if the variations along the tangential directions x, y are neglected. The gas temperature T must also be solved from the NSF equations combined with the separate temperature jump model. As shown by previous works [START_REF] Zheng | Comparison of kinetic theory and hydrodynamics for poiseuille ow[END_REF][START_REF] Xu | Super-burnett solutions for poiseuille ow[END_REF][START_REF] To | Molecular dynamics simulations of pressure-driven ows and comparison with acceleration-driven ows[END_REF], the gas temperature prole exhibits some anomalies that can not be described by available analytical solutions and even by compressible NSF equations. Furthermore, when the surface is isotropic α x = α y , Maxwell's model can not be recovered from the model proposed in Ref. [START_REF] Dadzie | Temperature jump and slip velocity calculations from an anisotropic scattering kernel[END_REF].

We must also say that the slip models are based on the validity of NSF equations at the wall which is not true. However, as shown in previous works and our MD results in the later section, these models can describe reasonably well the slip eects near the wall. Furthermore, errors can be dealt with using empirical coecients or higher order models, etc...

Poiseuille ow analytical solutions

We consider conned gas ows between two immobile parallel textured walls. The (dimensionless) slip tensors of each wall are denoted respectively b+ (upper wall) and b-(lower

The temperature dependence of the slip model can be expressed via the temperature dependence of αx, αy wall). The ows can be driven by either pressure dierence or external force eld. As shown in previous works [START_REF] Bhattacharya | Molecular dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases[END_REF][START_REF] To | A slip model for micro/nano gas ows induced by body forces[END_REF][START_REF] To | Molecular dynamics simulations of pressure-driven ows and comparison with acceleration-driven ows[END_REF][START_REF] Morris | Slip length in a dilute gas[END_REF] relevant to isotropic surfaces, Poiseuille parabolic solution can describe well the velocity prole even in nano-size channel and at high Knudsen number. These results suggest that we can adopt the following equation for the case under consideration [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF] u

= H 2 4η { 1 2 [ 1 - ( 2z
H ) 2 ] I + J + ( 2z
H ) K } g .
(5.26)

with g being the driving force, equivalent to the pressure gradient or the external body force. For MD simulation cases studied in the later section, the uid is subject to a gravity like eld γ and the driving force g can be evaluated as average mass density m g ρ 0 times γ g = m g ρ 0 γ

(5.27)

Tensors J, K are dened by

J = C -D(I + C) -1 D , K = (I + C) -1 D , C = b+ + b-, D = b+ -b-.
(5.28)

The rst term of (5.26) is the usual no-slip Poiseuille prole, the second constant term corresponds to slip-driven plug ow and the third term corresponds to a linear shear stress ow [START_REF] Bazant | Tensorial hydrodynamic slip[END_REF].

The general formulation (5.26,5.28) does not require that the two wall surfaces are aligned and their slip tensors b+ and b-are not necessarily diagonal in the current basis.

In the case where the slip tensors b+ , b-are diagonal in the basis Oxy, J, K are also diagonal in the same base and their components admit the expressions

J x = b + x + b - x + 4b + x b - x 1 + b + x + b - x , J y = b + y + b - y + 4b + y b - y 1 + b + y + b - y , K x = b + x -b - x 1 + b + x + b - x , K y = b + y -b - y 1 + b + y + b - y .
(5.29)

If the relation (5.23) holds, we can rewrite the previous expression as follows

J x = Kn(N + x + N - x ) + 4Kn 2 N + x N - x 1 + Kn(N + x + N - x ) , J y = Kn(N + y + N - y ) + 4Kn 2 N + y N - y 1 + Kn(N + y + N - y ) . K x = Kn(N + x -N - x ) 1 + Kn(N + x + N - x ) , K y = Kn(N + y -N - y ) 1 + Kn(N + y + N - y )
.

(5.30)

These analytical expressions will be examined in the later section with the Molecular Dynamics method.

Molecular dynamics models

In this section, we described the Molecular Dynamics method to simulate gas ows in a nanochannel. The gas atoms are placed in a rectangular box of dimension L × B × H along x, y, z directions and subject to a uniform force eld γ parallel to the plane xOy. We assume that the interaction between the gas atoms and the walls (normal to the z direction) can be described by the anisotropic kernel (A.1 -A.3, Appendix I). The following parameters associated to the kernel are used α x = 0.96 , α y = 0.83 , α z = 0.9 , T w = 300K .

(5.31)

It is noted that the values of α x , α y correspond to the tangential accommodations coecients for Ar and a nanotextured Pt wall at T w = 300K with dimensions described in Fig. 5.1 and in Ref. [START_REF] Pham | Eects of surface morphology and anisotropy on the tangential-momentum accommodation coecient between Pt(100) and Ar[END_REF]. In Ref. [START_REF] Pham | Eects of surface morphology and anisotropy on the tangential-momentum accommodation coecient between Pt(100) and Ar[END_REF], these values are obtained from MD simulations of beam scattering experiment.

Although (A.1 -A.3) are the constraint equations for the phase density f at the wall, we shall implement an analogous boundary conditions for MD methods that can mimic the eects of the former. From the composition of the scattering kernel B, we can say that µ i is the percentage of the molecules whose reection mechanism is governed by the elementary kernel B i . During MD simulations, when a gas atom crosses the wall boundary, we consider that it collides with the wall. Then, the same molecule is reinserted at the wall boundary with the same x, y coordinates and its velocity is reassigned on the basis of dierent elementary processes. At each collision event, a random number X with uniform distribution between [0, 1] is generated. Depending on the value of X and the arriving velocity v ′ , we shall decide the reection mechanism B i and generate the post collision velocity v accordingly, for example: -If 0 < X < µ 0 then the reection is specular (kernel B 0 ):

v x = v ′ x , v y = v ′ y , v z = -v ′ z . -If µ 0 < X < µ 0 + µ xy then the reection is specular along z, diusive along x, y direction (kernel B xy ): v x , v y ∼ N (0, √ kT w /m g ), v z = -v ′ z . -If µ 0 + µ xy < X < µ 0 + µ xy + µ xz then the reection is specular along y, diusive along x, z direction (kernel B xz ) v x ∼ N (0, √ kT w /m g ), v z ∼ R( √ kT w /m g ) and v y = v ′ y .
-etc... considering that the total sum of µ i equals 1 (see A.2, Appendix I).

Here N (0, √ kT w /m g ) is the normal distribution with zero mean and standard deviation √ kT w /m g and R( √ kT w /m g ) is the Rayleigh distribution with parameter √ kT w /m g . In addition to the previously described wall boundary conditions along the z direction, the simulation box is periodic along the x, y directions.

It is clear that the described anisotropic collision mechanism is a generalization of Maxwell's original mechanism. Here the specular and diusion reections are applied to three velocity components for dierent proportions of atom µ i in order to reproduce the momentum relative changes (5.7). Furthermore, due to a considerable molecular portion reecting diusively at the wall temperature T w , these boundary conditions correspond to a thermal wall allowing energy exchange. As the results, thermostats are no longer needed during MD simulation. A rather similar thermal wall approach can be found in Ref [START_REF] Bhattacharya | Molecular dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases[END_REF].

In our simulations, the interaction force between two gas molecules at distance r is governed by the Lennard-Jones potential :

V (r) = 4ϵ [ ( σ r ) 12 - ( σ r ) 6
] ,

(5.32)

with cuto distance r c = 2.5σ. For Argon, the reference energy ϵ and length σ are respectively ϵ = 1.67 × 10 -21 J and σ = 0.341 nm. The global number density is kept constant, ρ 0 = 0.035σ -3 while the number of the molecules N , the acceleration γ and others geometric parameters like L, B, H are changed as shown in Table 5.1. During the simulations, the Leapfrog-Verlet integration algorithm is employed and the averaging procedure starts only when the ow is stabilized, i.e. after about 10 6 time steps of 2 fs. The total computational time is 2 × 10 8 time steps. The height of the channel is divided into 100 layers to determine accurately the distribution of local stream velocity u. The ow direction determined by the angle φ made between the ow direction n (also the direction of γ) and the x axis is also varied from 0 • to 90 • to examine the anisotropy eect of the channel. The global mean free path λ used in the analysis is calculated by the formula

λ = 1 √ 2πρ 0 σ 2 .
(5.33)

Results

Walls with identical slip tensors

We consider rst the case where the two walls are identical and aligned, b

+ = b -= b.
This simplication leads to the following expression of J and K, J = 2b , K = 0 .

(5.34)

The MD simulation results conrm the parabolic velocity prole along the ow direction for most part of the channel, as predicted by (5.26). Another interesting aspect that agrees with (5.26) and some previous studies [START_REF] Belyaev | Eective slip in pressure-driven ow past super-hydrophobic stripes[END_REF][START_REF] Priezjev | Molecular diusion and slip boundary conditions at smooth surfaces with periodic and random nanoscale textures[END_REF] is that we observe the small transversal velocity prole. The latter seems to vanish at principal directions φ = 0 • , 90 • and maximal at φ = 45 • (see Fig. 5.3). We will show that these properties cannot be recovered by the isotropic slip theory.

Next let us look further into the quantitative aspect of the results. The velocity proles along the ow direction and the transversal direction are tted respectively with parabolic and constant equation to determine the tted slip velocities at the walls. These numerically computed values will be compared with analytical solution issued from (5.26) while accounting for (5.23,5.25,5.30). We remark that some deviations from the analytical solution is observed in a small region near the wall (the Knudsen layer) and the real slip velocities can be dierent from the tted velocities (see Fig. 5.2). The former are computed directly by averaging molecular velocity before and after collision with the walls, also considered in this work. However, for the sake of convenience, the terminology "slip velocities" without prex (tted or real) is used to refer to the tted velocities when no comparison is made between these two quantities.

According to (5.26), the dependence of dimensionless slip lengths on the ow orientation can be expressed by the relation All MD data can be tted perfectly with (5.35) showing excellent agreement with tensorial slip theory. On the other hand, expression (5.36) based on the quantitative estimation of b x , b y shows good agreement with MD results in slip ow regime Kn < 0.1 with average errors less than 5%. At higher Kn, considerable discrepancy is observed. The former have been reported in numerous works concerning the original Maxwell slip model and dierent correction coecients (either empirically or theoretically based) have been proposed to deal with these issues [START_REF] Karniadakis | Microows and Nanoows: Fundamentals and Simulation[END_REF][START_REF] Albertoni | Numerical evaluation of the slip coecient[END_REF][START_REF] Bhattacharya | Molecular dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases[END_REF].

L s (φ) H = b x cos 2 φ + b y sin 2 φ (5.35) or L s (φ) H = Kn ( 2 -α x α x cos 2 φ + 2 -α y α y sin 2 φ ) ( 5
Regarding the ratio between the slip velocities in transverse and longitudinal directions u sm /u sn , its dependency on the ow direction is given by (5.26):

u sm u sn = (b y -b x ) sin φ cos φ b x cos 2 φ + b y sin 2 φ
(5.37) (5.38) using the approximation (5.25). It is clear that this ratio always vanishes by assuming isotropy. For anisotropic surfaces, it is a function of φ and only vanishes when φ = 0 • , 90 • . Fig. (5.5) shows tted slip velocities obtained from MD simulations in comparison with analytical predictions (5.37) and (5.38). Again, MD results agree very well with the analytical expression (5.37) and less well with the quantitative expression (5.38). More precisely, (5.38) overestimate the tted transverse ow velocity. However, when we look into the real slip velocities in Fig. (5.6), there is a very good agreement and more interestingly, the ratio is quasi independent of Kn as predicted by (5.38).

According to previous works [START_REF] Vinogradova | Wetting, roughness and ow boundary conditions[END_REF][START_REF] Zhou | Anisotropic ow in striped superhydrophobic channels[END_REF], the ratio of averaged transverse and longitudinal ow rates is given by

⟨Q⟩ m ⟨Q⟩ n = 6(b y -b x ) sin φ cos φ 1 + 6(b x cos 2 φ + b y sin 2 φ) , (5.39)
where ⟨Q⟩ m and ⟨Q⟩ n are the averaged ow rate components along the directions m and n. This ratio can be obtained by integrating the velocity prole across the channel :

⟨Q⟩ m = ∫ H 2 -H 2 ⟨u m (z)⟩dz , ⟨Q⟩ n = ∫ H 2 -H 2
⟨u n (z)⟩dz .

(

Eects of Knudsen layer and normal tangential coecient α z on anisotropic ows

In the micro/nanoscale, the behavior of a gas ow near the wall is dominated by the eect of gas-surface interactions. This area is called the Knudsen layer : a thermodynamically nonequilibrium region extending a distance λ from the wall. Although the Navier-Stokes equations can accurately describe the ow behavior outside the Kn layer, they fail to predict the ow eld of the Kn layer, especially for high Kn number [START_REF] Dongari | Modeling of knudsen layer eects in micro/nanoscale gas ows[END_REF].

Figures (5.11) and (5.12) show the ratio of real and tted slip velocities for the case B -walls with dierent slip tensors. As shown in the gures, the real slip velocities are about (0.7 -0.8) of tted slip velocities at both two walls. These results agree with previous studies [START_REF] Lockerby | A wall-function approach to incorporating Knudsen-layer eects in gas micro ow simulations[END_REF], which indicated that the ratio of these two slip velocities is approximately equal to 0.8.

Finally, we discuss the validity of slip models based on the two tangential parameters α x , α y by considering the eect of α z on anisotropic ows and on slip behavior of gas at the walls. In previous sections, the normal accommodation coecient α z was chosen to be 0.9 for modelling anisotropic surfaces. The case A with identical walls was considered and only α z is varied while others parameters like α x , α y remained unchanged. The variation of slip length with α z is plotted in Fig. (5.13) for Kn = 0.05 and Kn = 0.104. The slip length increases slightly as α z increases and the variation is stronger for high Kn. However, the dierence between the cases is small (about 2%). We can conclude that slip models based on the two tangential parameters α x , α y adequately describe the phenomenon.

Conclusion

A multiscale study of boundary conditions at gas-solid interfaces have been conducted by using the Molecular Dynamics method to represent dynamic slip phenomena.

In Chapter 2, both pressure and acceleration driven methods are used to simulate gas ows in slip and transitional regimes. In the former method, a constant gravity -like force is applied to the gas particles. The latter method controls the kinetic pressure dierence between the inlet and the outlet. The main novelty lies in the modication of the periodic velocity conditions: the dierence in squared velocity between the ingoing and outgoing particles of the simulation domain is maintained constant. Plane Poiseuille ows were modeled and the results compared with approximate analytical solutions of the Navier-Stokes and energy equations reported in the literature.

In Chapter 3, in order to determine which empirical Pt-Ar pairwise potentials proposed in the archival literature may be used in MD simulations, a new Pt-Ar potential has been determined from periodic DFT calculations of the interaction between an argon atom and the solid Pt(111) surface using the CRYSTAL09 code. The energies have then been tted by a sum of Pt-Ar pairwise inter-atomic van der Waals type terms. The potential parameters, i.e. the interaction energy and the vertical Pt(111)-Ar harmonic vibrational frequency, are in good agreement with experimental data [START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF][START_REF] Zeppenfeld | Van Hove anomaly in the phonon dispersion of monolayer Ar/Pt (111)[END_REF]. This potential is very similar to the empirical one determined by Rameseyer et al. [START_REF] Han | Thermally-driven nanoscale pump by molecular dynamics simulation[END_REF] and Kulginov et al. [START_REF] Kulginov | An empirical interaction potential for the Ar/Pt (111) system[END_REF]. However this last one has the advantage of considering a physically correct repulsive term. The present and the Kulginov Pt-Ar potentials could be retained for future MD studies involving the Pt(111)-Ar couple. The tangential momentum accommodation coecient of argon gas colliding with a Pt surface has been determined by MD simulations. Argon atoms impinged on atomistic Pt surface and the post collision information have been collected in order to compute the directional σ dir t (θ, ϕ). We observed that σ dir t depends on the wall model, and on the Pt-Ar pairwise potential. We found that the well depth and the equilibrium position are the main potential characteristics that inuence the value of σ dir t . The surface is almost isotropic and the average of σ dir t over θ corresponds to the the gas-solid constant, σ t , of the Maxwell model. When averaging over the incident angle and using the new potential, a 0.42σ t value has been predicted for a wall temperature T w =300 K. The present multi-scale approach combines the advantages of the two computational methods: quantum electronic structure and Molecular Dynamics calculations. Since the accommodation coecient depends on the interaction strength, accurate determinations of the interaction parameters by accurate approaches are crucial. The Molecular Dynamics approach can reproduce the dynamics of the system and the collision trajectories, which allows to accurately compute the TMAC.

In Chapter 4, we have studied the eects of temperature, surface textures and anisotropy on the σ t coecient of Ar gas and Pt(100) surface. The computational model is based on the molecular beam experiments and constructed with the most accurate available potentials and interaction models. Unlike the previous chapter, the multi-body Quantum-Sutton Chen potential is used for the solid walls. Although its implementation is more costly than the harmonic spring potential used in Chapter 3, it better reproduces the surface eects. The MD results show that σ t is not a simply gas/wall constant, it decreases as the temperature increases and it increases with the roughness of the wall surface. The randomly rough surface obtained from the atomic deposition simulation is also investigated and we observed that, in addition to the roughness height, the in-plane random arrangement of the atoms also plays a signicant role on σ t values. Concerning the anisotropy eect, results on systems with anisotropic surfaces show that σ t varies signicantly with orientation. Eective σ t coecients are obtained and compared with models recently reported in the literature.

The last chapter of the thesis examines the gaseous slip ows over anisotropic surfaces from theoretical and numerical points of view. A simplied tensorial slip model based on tangential accommodation tensor is proposed and confronted with MD simulation results. In MD model, a special gas wall interaction based on anisotropic scattering kernels is implemented to mimic the surface anisotropy. The numerical results are shown to be in good agreement with theoretical predictions for the slip ow regime.

2AHIFA?JELAI

The results presented in this thesis could be completed in several directions

• The new pressure driven method using in Chapter 2 proved to be realistic. Hence, it could be applied to many extended problems with non constant axial pressure gradient: ows around obstacles or in rough channels, for example. The method appears thus relevant for modeling compressibility eects as well as temperature variations along the ow direction, a domain still unfullled by using MD methods.

• With the data obtained from the molecular beam simulations, we will easily determine the other important coecients characterizing the exchange between gas and surface such as the normal momentum accommodation coecient and the energy accommodation coecient (σ T ). The eects of the gas/wall interaction potential and the morphology of surface on these parameters could be also deduced. On the other hand, knowing the σ T coecient facilitate the quantication of the thermal slip in microchannels that not yet considered in this thesis.

• The current beam simulations model monoatomic gases and xed solid walls. It is also interesting to extend for polyatomic gases and moving walls.

• Although molecular dynamics is a very powerful technique, it is limited to small domains and small time frames. Nowadays, a MD simulation can be performed on systems containing a billion atoms at best, with dimensions up to micron. A possible solution to this limitation is to combine molecular dynamics with continuum mechanics. The results presented in this thesis such as the accommodation coecients, slip velocities could be used as boundary conditions at gas/solid interface for the continuum mechanics simulation in remainder region where the continuum description is appropriate.

Since the velocity of a gas atom is unchanged between collisions, it is possible to think that ⟨v ′ t ⟩ is equal to the velocity average at one mean free path λ from the wall, hence the stream velocity u λ at that location. By assuming that N -= N + = N/2, the slip velocity u s can be deduced as follows Considering the rst order λ approximation, i.e u λ = u s + λ ∂u ∂z , the nal slip equation becomes The term standing before ∂u ∂z is equivalent to the slip length tensor b and identical to (5.23).

List of Figures

NomenclatureB

  (v ′ , v) scattering kernel [d.u.] b slip length tensor [m] C p , C v heat specic at constant pressure and volume [J.K -1 .Kg -1 ] C 1 , C 2 rst and second order slip coecients [d.u.] c thermal speed [m.s -1 ] Greek symbols α x , α y , α z momentum accommodation coecients along x, y, z axes[d.u.] δt time step [s] dΩ solid angle [sr] ϵ depth of the potential well [J] η viscosity [Kg.m -1 .s -1 ] γ external force eld [Kg.m.s -2 ] κ thermal conductivity [W.m -1 .K -1 ] λ mean free path [m] µ i percentage of the elementary kernel B i [d.u.] Φ -, Φ + incoming and outcoming ux at the wall ψ 1 , ψ 2 , ... functions in Chapman and Enskog model ω vertical harmonic vibrational frequency [J] ρ number density [m -3 ] σ distance at which the inter-particle potential is zero in Lennard Jones potential [m] σ t tangential momentum accommodation coecient [d.u.] σ n normal momentum accommodation coecient [d.u.]

Figure 2 :

 2 Figure 2: Examples of micro/nanouidic systems.

Figure 3 :

 3 Figure 3: This chip is commercialized by Agilent Technologies and permits the identication of specic genetic sequences in a 1-µl sample of roughly puried DNA. The process time is just about 10 min [1].

Figure 1 . 1 :

 11 Figure 1.1: Elementary volume in physical and velocity space.

Figure 1 . 3 :

 13 Figure 1.3: Classication of the gas ow regimes and governing equations over the range of Knudsen numbers (Zhang et al., 2012)[13].

Figure 1 . 4 :

 14 Figure 1.4: Specular reection and diusive reection.

Figure 1 . 5 :

 15 Figure 1.5: Interpretation of the Navier slip length L s for Couette velocity proles

Figure 1 . 7 :

 17 Figure 1.7: Typical initial congurations for the solid atoms

Figure 1 . 11 :

 111 Figure 1.11: Diagram of an MPI application.

Figure 1 . 13 :

 113 Figure 1.13: Communication between two neighboring processors.

Figure 2 . 2 :

 22 Figure 2.2: Pressure distribution (p(x)) along the channel for pressure driven ow case and acceleration driven ow case (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ, γ = 4.3 × 10 -4 ϵ 0 /mσ). The full line represents the linear t and is used to determine the pressure gradient ∆p/∆x. The pressure unit is [ϵ 0 /σ 3 ].

Figure 2 . 4 : 1 Figure 2 . 5 :

 24125 Figure 2.4: Distribution of pressure Pxx along the channel for pressure driven ow case (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ) computed by Eq. 2.20 and 2.21. The pressure unit is [ϵ 0 /σ 3 ].

Figure 2 .Figure 2 . 6 :

 226 Figure 2.6: Dimensionless velocity prole in half of the channel cross-section for Kn=0.05 and Kn = 0.1, (r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The full lines represent equation v(ẑ) = v(0)(1 -

  driven, Kn = 0.1 Acceleration driven, Kn = 0.1 Pressure driven, Kn = 0.05 Acceleration driven, Kn = 0.05

Figure 2 . 7 :

 27 Figure 2.7: Dimensionless temperature prole ( T (z)) in half of the channel cross-section for two Kn-values according to the ow model (r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ).

Figure 2 . 8 :

 28 Figure 2.8: Variation of temperature prole when increasing the pressure gradient or parameter δ (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The dashed lines are based on equation c 0c 2 ẑ2c 4 ẑ4 that best ts the numerical results. c 0 = 3.07, c 2 = 4.16, c 4 = 12.8 for γ = 0.1ϵ 0 /mσ, c 0 = 1.74, c 2 = 2.5, c 4 = 1.0 for γ = 0.05ϵ 0 /mσ.

Figure 2 . 9 :

 29 Figure 2.9: Changes in the temperature prole for various values of the acceleration parameter γ (Kn = 0.05,r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The solid and the dashed lines correspond to equations b 0b 4 ẑ4 and d 0d 2 ẑ2d 4 ẑ4d 6 ẑ6 that best t the numerical results. b 0 = 2.53, b 4 = 13.0, d 0 = 2.46, d 2 = 0.54, d 4 = 4.58, d 6 = -97.8 for γ = 0.1ϵ 0 /mσ and, b 0 = 1.48, b 4 = 4.57 for γ = 0.05ϵ 0 /mσ.

Figure 2 . 10 :

 210 Figure 2.10: Axial variations of pressures p(x) = p K + p V and p K for pressure driven ow (n = 0.1σ -3 , r c = 2.5σ, ϵ = ϵ 0 ). The pressure unit is [ϵ 0 /σ 3 ].

Figure 2 . 11 :

 211 Figure 2.11: Axial variations of pressures p(x) = p K + p V and p V for pressure driven ow (n = 0.8σ -3 , r c = 2 1/6 σ, ϵ = ϵ 0 ). The pressure unit is [ϵ 0 /σ 3 ].

Figure 2 . 14 :

 214 Figure 2.14: Streamlines in a rib-roughened channel. The color represents the uid density eld.

Figure 3 . 1 :

 31 Figure 3.1: Scheme of the Pt cluster used to decompose the global Pt(111)-Ar interaction potential on a sum of two bodies Pt-Ar terms used for the molecular dynamics. The circle represents the argon atom adsorbed on the fcc-hollow site, numbers 1 represent platinum atoms of the rst layer, numbers 2 Pt atoms of the second layer, and numbers 3 Pt atoms of the third layer, respectively.

Figure 3 . 2 :

 32 Figure 3.2: Molecular dynamics scheme. The incident argon atoms are with v i velocities with no contribution along y. θ is the incident angle. The Pt wall has a fcc structure with a (111) surface. Each Pt atom is connected to the nearby ones by harmonic springs of rigidity k w (see the text).

  8 N/m, δt = 2 fs, ξ = 5.184 × 10 -12 kg/s.

Figure 3 .

 3 Figure 3.3 presents the Pt(111)-Ar potential computed together with those discussed in Section 2.1. Four types of potentials can be distinguished. The rst ones are the Lennard-Jones potentials established by Maruyama and Kimura[START_REF] Kimura | A molecular dynamics simulation of a bubble nucleation on solid surface[END_REF][START_REF] Maruyama | A study on thermal resistance over a solid-liquid interface by the molecular dynamics method[END_REF] and by Spijker et al.[START_REF] Spijker | Computation of accommodation coecients and the use of velocity correlation proles in molecular dynamics simulations[END_REF]. They depict a relatively weak value of V e , of ≈ -40 meV, and a short equilibrium z e < 3 Å. The second group of potentials includes those of Head-Gordon et al.[START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] and Svanberg

Figure 3 . 3 :

 33 Figure 3.3: Pt(111)-Ar interaction potentials. z is the distance between the Ar atom at the vertical of the fcc-hollow site and the rst layer of the Pt[START_REF] Asproulis | Boundary slip dependency on surface stiness[END_REF] surface. z = 0 is dened as the position of the Pt nuclei of this rst layer for T = 0 K.

Figure 3 . 4 :

 34 Figure 3.4: Pt-Ar pairwise potentials.

Figure 3 . 5 :

 35 Figure 3.5: Pt-Ar pairwise forces.

Figure 3 .

 3 Figure 3.6: σ dir t (θ, ϕ) values computed for T w = 300 K, θ = 45 • using the present potential, for ϕ varying between 0 and 90 • .

Figure 3 . 7 :(

 37 Figure 3.7: Distribution of the v rx velocity with v ix = 262 m.s -1 after 12000 collisions using the present pair potential for T w = 300 K, θ = 45 • and ϕ = 0 • . The red curve corresponds to the t of the results by a gaussian function f (v rx ) = 7.6 × 10 -4 +

Figure 3 . 8 :

 38 Figure 3.8: Norm of the reected velocity (in m.s -1) recorded for dierent values of the reected angle θ r (in degrees) with respect to the normal direction of the surface. The present pair potential has been used as well as T g = 600 K, T w = 500 K, θ = 45 • and ϕ = 0 • . The dotted curve corresponds to the experimental values recorded by Head-Gordon et al.[START_REF] Head-Gordon | On the nature of trapping and desorption at high surface temperatures. Theory and experiments for the ArPt (111) system[END_REF] and the markers to the simulated data from the same work.

  4.1). As shown by previous works[109111], the lattice constant, wall mass and stiness can have signicant impacts on σ t and the slip eects.

Figure 4 . 1 :

 41 Figure 4.1: Surface eects: the fully relaxed conguration (right) is dierent from initial conguration (left).The solid lm system is composed of xed atoms (bottom layer), thermostat atom (upper bottom layer) and normal atoms (remaining layers).

Figure 4 . 2 :

 42 Figure 4.2: Nanotextured surface of type A (square)

Figure 4 . 3 :

 43 Figure 4.3: Nanotextured surface of type B (strip)

Figure 4 .

 4 [START_REF] Liou | Microuid Mechanics[END_REF] shows a snapshot of the nal system whose total number of Pt atoms have reached 733. To improve the statistical results, 5 samples obtained thanks to the above described procedure are collected, as shown in Fig.4.5.

Figure 4 . 4 :Figure 4 . 5 :Figure 4 . 6 :

 444546 Figure 4.4: Snapshot of deposition process (left) and nal thin lm system (right)

Figure 4 . 7 :

 47 Figure 4.7: Molecular dynamics scheme. The incident argon atoms are with v i velocities. θ is the incident angle. The Pt wall has a fcc structure with a (100) surface. The Pt atoms are controlled by Sutton-Chen potential.

Figure 4

 4 Figure 4.8: σ dir t computed for the wall of type A (square) at T w = 200 K, 250 K, 300 K, 350 K and 400 K for three roughness heights h with θ = 45 • and φ = 0 • .

  .3 and Fig.

  4.10). The real behaviors are mixed: we sometimes observe the colliding and bouncing mechanisms on rough surfaces (not shown in Fig.4.10), but they are not typical.

Figure 4 . 9 :

 49 Figure 4.9: Typical collision trajectories (solid and dashed lines) on a smooth surface. Gas molecules collide and bounce several times before escaping.

Figure 4 . 10 :

 410 Figure 4.10: Typical collision trajectories (solid and dashed lines) on a rough surface. Gas molecules move within the valley between the peaks.

Figure 4 .

 4 Figure 4.11: σ dir t computed for type-B walls (strip) versus azimuth angle φ for dierent roughnesses (T w = 300K, θ = 45 • ). The solid, dashed, dash-dotted lines are the analytical expressions (4.8) used to t the present numerical results.

Figure 4 .

 4 Figure 4.12: σ an t computed for type-B walls (strip) versus azimuth angle for dierent roughnesses (T w = 300K). The solid, dashed, dash-dotted lines are the analytical expressions (4.8) used to t the present numerical results.

Figure 5 . 2 :

 52 Figure 5.2: Schematic of the Knudsen layer and real and tted slip velocities.

. 36 )

 36 if the simplied quantitative relation (5.25) is used. From Fig.5.4, both analytical solutions and numerical solutions at dierent Kn show the same trend of L s in function of φ.

Figure 5 . 3 :Figure 5 . 4 : 19 Figure 5 . 5 :

 53541955 Figure 5.3: Longitudinal and transverse velocity proles u n , u m for dierent values of φ and Kn = 0.104. The velocities are normalized with u max -the velocity at z = 0 for case φ = 90 • .

Figure 5 . 10 :

 510 Figure 5.10: Normalized slip velocities at the lower wall as a function of φ for dierent Kn numbers. The lines represent the analytical expression (5.45) and the points are MD data.

Figure 5 . 13 :

 513 Figure 5.13: Dimensionless slip length L s /H in term of α z with φ = 0 • for Kn = 0.05 and Kn = 0.104.

+ λ 2 2 ∂ 2 u ∂z 2 +

 22 2u s = (2I -A)u λ .(A.[START_REF] Cao | Molecular momentum transport at uid-solid interfaces in MEMS/NEMS: A Review[END_REF] We remark that the Taylor development of u at the wall yields the expressionu λ = u s + λ ∂u ∂z ... . (A.15)

  u s = λ(2A -1 -I)

1

  An inkjet printhead (left) and the visualization of droplets of ink projected onto a target (right)[START_REF] Tabeling | Introduction to microuidics[END_REF],[START_REF] Ho | ASME Int. Mechanical Engineering Congress and Exposition[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Examples of micro/nanouidic systems. . . . . . . . . . . . . . . . . . . . . 3 This chip is commercialized by Agilent Technologies and permits the identication of specic genetic sequences in a 1-µl sample of roughly puried DNA. The process time is just about 10 min [1]. . . . . . . . . . . . . . . . 1.1 Elementary volume in physical and velocity space. . . . . . . . . . . . . . . 1.2 Phase space element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Classication of the gas ow regimes and governing equations over the range of Knudsen numbers (Zhang et al., 2012)[13]. . . . . . . . . . . . . . . . . . 1.4 Specular reection and diusive reection. . . . . . . . . . . . . . . . . . . . 1.5 Interpretation of the Navier slip length L s for Couette velocity proles . . . 1.6 Periodic boundary conditions for a simulation box in two dimensions . . . . 1.7 Typical initial congurations for the solid atoms . . . . . . . . . . . . . . . . 1.8 All-pairs method for interaction computations : in a simulation box, all pairs of particles should be examined. . . . . . . . . . . . . . . . . . . . . . 1.9 Representation of a simulation box in cell subdivision method in two dimensions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Truncated sphere interactions related to r c and neighbor list related to r v , of a given particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Diagram of an MPI application. . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Spatial decomposition scheme in two dimensions : each cell is handled by a dierent processor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Communication between two neighboring processors. . . . . . . . . . . . . . 2.1 Modied periodic boundary conditions. The velocity v between the outgoing molecules are dierent from the ingoing molecules (Eq. 2.16 and 2.17) . . . 2.2 Pressure distribution (p(x)) along the channel for pressure driven ow case and acceleration driven ow case (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ, γ = 4.3 × 10 -4 ϵ 0 /mσ). The full line represents the linear t and is used to determine the pressure gradient ∆p/∆x. The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Distribution of pressure p(x) and pressure components Pxx (x), Pyy (x), Pzz (x) along the channel for pressure driven ow case (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ). The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . . . 2.4 Distribution of pressure Pxx along the channel for pressure driven ow case (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ) computed by Eq. 2.20 and 2.21. The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . . . . . . . . . . . . . 2.5 Dimensionless density prole (n(z)) in half of the channel cross-section for pressure driven simulations (Kn=0.05, 0.1, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 , δ = 10 -3 ϵ 0 /mσ, 2 × 10 -3 ϵ 0 /mσ). . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Dimensionless velocity prole in half of the channel cross-section for Kn=0.05 and Kn = 0.1, (r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The full lines represent equation v(ẑ) = v(0)(1a 2 ẑ2 ) that best t the numerical results (a 2 = 0.18 for Kn=0.05 and a 2 = 0.29 for Kn=0.1). . . . . . . . . . . . . . . . . . . . . . . 2.7 Dimensionless temperature prole ( T (z)) in half of the channel cross-section for two Kn-values according to the ow model (r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). . . . 2.8 Variation of temperature prole when increasing the pressure gradient or parameter δ (Kn = 0.05, r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The dashed lines are based on equation c 0c 2 ẑ2c 4 ẑ4 that best ts the numerical results. c 0 = 3.07, c 2 = 4.16, c 4 = 12.8 for γ = 0.1ϵ 0 /mσ, c 0 = 1.74, c 2 = 2.5, c 4 = 1.0 for γ = 0.05ϵ 0 /mσ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Changes in the temperature prole for various values of the acceleration parameter γ (Kn = 0.05,r c = 2 1/6 σ, ϵ = 0.1ϵ 0 ). The solid and the dashed lines correspond to equations b 0b 4 ẑ4 and d 0d 2 ẑ2d 4 ẑ4d 6 ẑ6 that best t the numerical results. b 0 = 2.53, b 4 = 13.0, d 0 = 2.46, d 2 = 0.54, d 4 = 4.58, d 6 = -97.8 for γ = 0.1ϵ 0 /mσ and, b 0 = 1.48, b 4 = 4.57 for γ = 0.05ϵ 0 /mσ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Axial variations of pressures p(x) = p K + p V and p K for pressure driven ow(n = 0.1σ -3 , r c = 2.5σ, ϵ = ϵ 0 ). The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . 2.11 Axial variations of pressures p(x) = p K + p V and p V for pressure driven ow (n = 0.8σ -3 , r c = 2 1/6 σ, ϵ = ϵ 0 ). The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . 2.12 Dimensionless velocity prole in half of the channel cross-section for n = 0.1σ -3 and n = 0.8σ -3 . The dashed lines represent the equations a 0a 2 z 2 that best t the numerical results. . . . . . . . . . . . . . . . . . . . . . . . 2.13 Axial variations of pressure component Pxx computed by Eq. 2.20 and 2.21. The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Streamlines in a rib-roughened channel. The color represents the uid density eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Distribution of pressure component Pxx in a rib-roughened channel computed by Eq. 2.20 and 2.21. The pressure unit is [ϵ 0 /σ 3 ]. . . . . . . . . . . 3.1 Scheme of the Pt cluster used to decompose the global Pt(111)-Ar interaction potential on a sum of two bodies Pt-Ar terms used for the molecular dynamics. The circle represents the argon atom adsorbed on the fcc-hollow site, numbers 1 represent platinum atoms of the rst layer, numbers 2 Pt atoms of the second layer, and numbers 3 Pt atoms of the third layer, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Molecular dynamics scheme. The incident argon atoms are withv i velocities with no contribution along y. θ is the incident angle. The Pt wall has a fcc structure with a (111) surface. Each Pt atom is connected to the nearby ones by harmonic springs of rigidity k w (see the text). . . . . . . . . . . . . 3.3 Pt(111)-Ar interaction potentials. z is the distance between the Ar atom at the vertical of the fcc-hollow site and the rst layer of the Pt(111) surface. z = 0 is dened as the position of the Pt nuclei of this rst layer for T = 0 K. 3.4 Pt-Ar pairwise potentials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Pt-Ar pairwise forces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 1: Some relations between macroscopic physical quantities Q and molecular quantities Q.

Table 1 .

 1 2: Approach equations of CE model for dierent order k of Kn.

	Order k	Model equations
	Kn 0	Euler
	Kn 1

Table 1

 1 .3 presents their values when σ t = 1. There is no general agreement

Table 1 .

 1

		3: Values of rst and second order slip coecients in literature
	References	Approach	C 1	C 2

  C 1 is around unity while C 2 varies from -0.5 to 9/8. In Table1.4, the expressions of two slip coecients as function of σ t are shown. These values could be derived from the Navier-Stokes equation as well as the Boltzmann equation. The rst-order slip coecients C 1 usually relate to σ t , while the second-order slip coecients C 2 are rather dierent.

			2/3	0.25
	Hadjiconstantinou (2003) [28]	HS model	1.11	0.61
	Cercignani and Lorenzani (2010) [29]	HS model 1.1209 0.2347
		BGK-based model 1.1366 0.69261
	Li et al. (2011) [30]	Lattice Boltzmann model	1	0.8
	on the results :

Table 1 .

 1

	4: Expressions of rst and second order slip coecients in literature
	References			C 1	C 2
	Maxwell (1879) [16]			2-σt σt	0
	Loyalka (1971) [21]	2-σt σt	√ π 2 (1 -0.1871σ t )	0
	Cercignani (1988) [31]	2 √ π	2-σt σt (1 + 0.1621σ t )	2 π (0.5 + C 2 1 )
	Jie et al. (2000) [32]			2-σt σt	-2-σt σt	Re Kn
	Aubert and Colin (2001) [33]			2-σt σt	9/8
	Karniadakis and Beskok (2002) [26]			2-σt σt	-2-σt 2σt
	Lockerby et al. (2004) [34]

Table 2 .

 2

		Boundary conditions
	Faces x	Eqs. (2.16-2.17)
	Faces y	PBC
	Faces z

1: Summary of the boundary conditions applied on the simulation domain

Table 2 .

 2 .3. 3: Input data for MD simulations in Section 2.5.2

	n

  Figure 2.13: Axial variations of pressure component Pxx computed by Eq. 2.20 and 2.21. The pressure unit is [ϵ 0 /σ 3 ].

	IK1 (Eq. 2.20)
	MOP (Eq. 2.21)

Table 3 .

 3 1: Pt(111)-Ar interaction potential parameters. -1 ) C vdW ( eV.Å 3 ) z 0 (Å)

	V min (eV)

∑ Rep + vdW V 0 (eV) α (Å r = |r Arr P t |.

Table 3 .

 3 2: Equilibrium interaction energy (V e in meV), equilibrium distance (z e in Å), and estimated vertical harmonic vibrational frequency (ω e in meV) of the Pt(111)-Ar interaction potentials.

	Potential Present	V e / meV z e / Å ω e / meV -81.3 3.35 4.8
	Kulginov [71]	-89.7	3.24	4.8
	Ramseyer [82]	-72.8	3.22	5.0
	Maruyama [73, 74]	-38.4	2.93	3.9
	Spijker [76]	-43.2	2.76	8.1
	Head-Gordon [68]	-76.0	4.11	5.7
	Svanberg [81]	-77.2	3.89	4.4
	Yamamoto [75]	-41.8	4.28	8.8
	Smith [69]	-104.0	3.00	5.6
	Lahaye [70]	-103.4	2.98	8.7
	Exp. [68, 79]	-80.0 a		5.0 b

a Estimated V 0 value. b Anharmonic vibrational frequency.

Table 3 . 3 :

 33 Comparison of the σ dir t (θ = 45 • , 0 • ) given by two models for T = 300 K and an incident angle of 45 • .

	Potential	D e / meV R e / Å Model 1 Model 2
	Present work	9.5	3.93	0.40	0.14
	Yamamoto [75]	11.6	4.60	0.28	0.14
	sive potential with respect to the attractive one.

R 0 has been derived from local-density approximation DFT calculations.

Table 3 .

 3 4: Dissociation energy (D e in meV), equilibrium distance (R e in Å), and σ dir t (θ = 45 • , 0 • ) computed at T w = 300 K with an incident angle of 45 • from dierent Pt-Ar pairwise potentials.This value is comparable to the σ t value of 0.19 found by Yamamoto[START_REF] Yamamoto | Slip ow over a smooth platinum surface[END_REF] using a Couette ow simulation approach and a simplied wall description. In his simulation, the Pt atoms are connected through Lennard-Jones pairwise potentials with no possible relaxation. We nd that Model 1 gives a much higher σ dir t (45

	Potential Present work	D e / meV R e / Å σ dir t (θ = 45 • ) 9.5 3.93 0.40
	Kulginov [71]	10.2	3.85	0.46
	Ramseyer [82]	9.1	3.76	0.49
	Maruyama [73, 74]	5.6	3.46	0.34
	Spijker [76]	6.8	3.30	0.28
	Head-Gordon [68]	11.6	4.60	0.24
	Svanberg [81]	6.8	4.49	0.23
	Yamamoto [75]	11.6	4.60	0.28

Table 3 .

 3 5: σ dir t (θ, 0 • ) and corresponding standard deviation σ T M AC computed at T w = 300 K with the Pt-Ar pairwise potential of the current work for dierent incident angles.

	Incident angle σ dir t (θ, 0 • ) σ T M AC 0.43 2.30 10 • 0.41 0.85 30 • 0.40 0.63 45 • 60 • 0.41 0.55

  • , 0 • ) = 0.34 for 300 K of Table 3.4. At 200 K, the present model would leads to σ dir

t

[START_REF] Todd | Pressure tensor for inhomogeneous uids[END_REF]

Table 4 .

 4 3: σ dir t (θ, φ) computed at T w = 200 K, 300 K, 400 K for three roughness heights h with θ = 10 • , 45 • , 80 • and φ = 0 • .

		Surface type		θ T w = 200 K T w = 300 K T w = 400 K
		A (h =5.88 Å) 10 • 45 • 80 •	0.96 0.92 0.90	0.87 0.85 0.83	0.79 0.77 0.74
		A (h =3.92 Å) 10 • 45 • 80 •	0.94 0.90 0.88	0.84 0.79 0.78	0.75 0.74 0.72
		Smooth	10 • 45 • 80 •	0.85 0.82 0.80	0.72 0.70 0.69	0.61 0.60 0.59
		Table 4.4: σ iso t	and σ dir
			Surface type A (h =5.88 Å) 0 • 10 • 0.87 φ θ σ dir t 0 • 30 • 0.86 0 • 45 • 0.85 0 • 60 • 0.85 0 • 80 • 0.83	σ iso t
						0.85
			Random		0.92
			(Fig. 4.5)
			Smooth		0.70
	4.3	MD simulation results

t (θ, φ) computed at T w = 300 K.

Table 4 .

 4 5: Ratio ⟨v rm ⟩/⟨v rn ⟩ computed for type-B walls (strip) with dierent roughness heights h at T w = 300K, θ = 45 • and φ = 45 • .

	h 0	⟨v rm ⟩/⟨v rn ⟩ 0
	1.96 Å	0.15
	3.92 Å	0.39
	5.88 Å	0.67

  .[START_REF] Loyalka | Approximate method in the kinetic theory[END_REF] while cz can be estimated from(5.16) by considering f (hence f + , f -) close to Maxwellian

	√
	cz =	2kT πm g	.	(5.22)
	As a result of (5.20), (5.19) and (5.11), tensor b becomes simplied into the simple form
	b

Table 5 .

 5

		1: Input data of MD simulations
	N (atoms) L(σ) B(σ) H(σ) Kn	γ(ϵ/mσ)
	74088	128.4 128.4 128.4 0.05 0.001
	27000	91.7 91.7 91.7 0.07 0.001
	8000	61.1 61.1 61.1 0.104 0.001
	4400	61.1 61.1 33.6 0.19 0.001

v ′ = αv, α = √ 1 + 3δL/v 2 .(2.16)

Remerciements

List of gures

In the case of simplied model, equation (5.39) can be expressed in terms of accommodation coecients : (5.41)

The ratio ⟨Q⟩m ⟨Q⟩n as a function of φ for dierent values of Kn is plotted in Fig. (5.7). The results show that the MD data and analytical solutions are closely similar : the ratio is largest at φ = 45 • and it approaches zero when φ tends to 0 • or 90 • . However, there are some visible discrepancies between the MD data and analytical results especially, at high Kn number, as it has been explained earlier.

Walls with dierent slip tensors

In this section, we consider the second case where the two walls are misaligned. The lower wall position is maintained while the upper wall is rotated with an angle 90 • in its plane. In this case we have

which leads to the following expression of J, K

(5.43)

If the estimation (5.25) of b from the accommodation coecients is adopted, one obtains

.

(5.44)

Normalizing the slip velocities at two walls u + sn , u - sn (at z = ±H/2) with the mid stream velocity u n0 (at z = 0) yields

(5.45)

Figure (5.8) shows a good agreement between the analytical solutions and MD results: the longitudinal velocity prole is a shifted parabola. As expected, near the wall where the Knudsen layer dominates, the velocity distribution tends to deviate from (5.26). The normalized slip velocities at the upper and lower walls are plotted in Fig.

(5.9, 5.10)as a function of φ . As can be seen, the agreement between MD results of tted slip velocities and predicted solutions (5.45) is also observed at small Kn numbers for both two walls. For gas-wall interaction, Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] proposed an anisotropic scattering kernel B(v ′ , v) dened by

in which

The vectors v ′ , v are respectively the arriving velocity and the reected one and the constants α x , α y , α z are the accommodation coecients along the directions x, y, z. The elementary kernels B k are given by the following expressions

with δ being the delta Dirac function, T w the wall temperature, k B the Boltzmann constant and m g the molecular mass of gaseous particles. The boundary conditions for particle distribution function f (v) is then dened by

We use Φ - j and Φ + j to denote the incoming ux at the wall of the momentum j component. Then

with f -and f + being the velocity distribution associated with the incident molecules and reected molecules. Dadzie and Meolans [START_REF] Dadzie | Anisotropic scattering kernel: Generalized and modied Maxwell boundary conditions[END_REF] proved the following relation

Their model is based on three parameters α x , α y , α z dened along given directions of a system of coordinate. We are interested in the accommodation coecients in an arbitrary direction. Hence, we consider a family of orthogonal directions (n, m) obtained by rotating xOy around Oz by an angle φ. Consequently, the n, m-components are related to x, y components by

(A.7)

Since α x and α y are accommodation coecients, we can deduce the the relation

with

and compute the accommodation coecient along any direction n. For example, by setting the component Φ - m = 0 (e.g we beam atoms along direction n only), we can recover the expression for α n in Eq. (4.8). The ratio between the reected components m, n can also be computed by the expression

It is clear that for isotropic model α x = α y , this ratio is always zero for all φ. Thus for anisotropic surface α x ̸ = α y , the ratio Φ + m /Φ + n is a function of φ, only vanishes at φ = 0 • , 90 • . For example, at φ = 45 • , we obtain

(A.11)

Appendix II. Slip length tensor

We consider a control surface near and parallel to the immobile wall and study the collisions at this surface. During a unit period of time, there are N gas atoms crossing the surface: Natoms go downward and N + atoms go upward with tangential velocity

and v t = (v x , v y ). The gas velocity u s at the wall can be obtained by the average

(A.12)

The notation ⟨ ⟩ indicates the average over the atoms. It is clear that mN + ⟨v t ⟩ and mN -⟨v ′ t ⟩ are equivalent to the tangential momentum vector uxes Φ + t and Φ - t . Using the linear relation (5.8), we can write θ is the incident angle. The Pt wall has a fcc structure with a (100) surface. The Pt atoms are controlled by Sutton-Chen potential. . . . . . . . . . . . . 4.8 σ dir t computed for the wall of type A (square) at T w = 200 K, 250 K, 300 K, 350 K and 400 K for three roughness heights h with θ = 45 • and φ = 0 • . . 4.9 Typical collision trajectories (solid and dashed lines) on a smooth surface.

Gas molecules collide and bounce several times before escaping. . . . . . . . 4.10 Typical collision trajectories (solid and dashed lines) on a rough surface.

Gas molecules move within the valley between the peaks. . . . . . . . . . . . 4.11 σ dir t computed for type-B walls (strip) versus azimuth angle φ for dierent roughnesses (T w = 300K, θ = 45 • ). The solid, dashed, dash-dotted lines are the analytical expressions (4.8) used to t the present numerical results. 4.12 σ an t computed for type-B walls (strip) versus azimuth angle for dierent roughnesses (T w = 300K). The solid, dashed, dash-dotted lines are the analytical expressions (4.8) used to t the present numerical results. . . . . 5.1 Equivalent gas-wall interaction model. In modelling and simulations, the structured surface is replaced by a smooth surface with a scattering kernel R(v ′ , v). In Ref. [START_REF] Pham | Eects of surface morphology and anisotropy on the tangential-momentum accommodation coecient between Pt(100) and Ar[END_REF], the authors considered a (100) crystalline Pt textured surface at 300 K with the following dimensions l 1 = 11.76Å, l 2 = 19.6Å, l 3 = 7.84Å, h = 5.88Å and used MD method to determine the accommodation parameters with Ar atoms associated to model (A.

List of publications