H D Bui

A Ehrlacher

Developments of fracture mechanics in France in the last decades

ABSTRACI': 1be purpose of the present paper is to give a survey on some developments of fracture mechanics in France in the last decades. Until about 1970, the sut.. ject was practically unknown at the universities. The topics were growing up rapidly, because of French industrial projects or programs in the field~ of aeronautics. aerospace industries, nuclear reactor technology. In France. fracture mechanics became an important subject for engineers and researchers starting witb the first Summer School ( 1970) devoted to fracture mechanics held at La Celie Saint Loup (France) (Fran~ois & Joly 1972). At this meeting, engineers from the French Atomic Centre (CEA), Electricite de France, the Nuclear Vessel Builders (Framatome, Creusot-Loire) and from other industries gave lectures on linear elastic fracture mechanics, Griffirh' s energy balance, Irwin's energy release rate G and Kies' stress intensity factors. Obviously~ fr.acture mechanics is strongly motivated by industrial applications. Ten years after the first meeting at La Celie Saint Loup~ scientists and engineers in France organised the ICF-S conference in Cannes, in April 1981. Contributions to the experimental determination of toughness of different materi-al~ the study of the influence of hydrostatic pressure on the toughness of aluminium alloys (Auger & F~is 1974), microscopic observations of cracks in meaa.Is. fatigue crack propagation [START_REF] Lemaitre | Aspect pheno~nologique de Ia rupture parendommagement[END_REF] date back to the 1970's. But most of the papers in this period are devoted to computational aspects in linear and non-linear fracture mechanics with the emphasis on stress intensity factors calculations. 3D calculatioo of SIF by the weight function method [START_REF] Labbens | Weight functions for 3D synaetrical crack problems[END_REF]), path-independent integrals [START_REF] Bonnet | Regularization of the displacement and traction BJE for 30 ela~ todynamics using indi~ct methods[END_REF], fatigue crack propagation and fatigue criterioo (Dang Van 1973), crack tip plastic zones (Nguyen & Bui 1974).

In the 1980's, studies in non-linear fracture mechanics were concerned with ductile fracture, with applications to nuclear energy technology (Rousselier 1981 ), crack initiation using damage mechanics [START_REF] Bui | Propagation of damage in elastic and plastic solid~[END_REF]Lemaitte & Cbaboche 1985), lhcnnodynamical aspects of fracture (Nguyen 1980;Ehrlacher 1985), bomogenisation of periodic structures containing cavities or cracks [START_REF] Marigo | Fonnulation d•Uile loi d•endommagement d'un materiau elastique. Comp1es Rendu[END_REF][START_REF] Marigo | Modelling of brittle and fatigue damage for elastic material by growth of micro[END_REF]Sanchez-Palencia 1980;[START_REF] Suquet | Une approche simplili~ pour 1c calcul des proprietes elastique.<~ de materiaux heterogenes a structure periodique[END_REF]). The contributions of Beremin, a pseudonym for a group of researchers from the Ecole des Mines de Paris. working in ductile fracture and the phy~ics of metals [START_REF] Mudry | Etude de Ia rupture ductile et de Ia rupture par clivage d'acier faiblement aiU~[END_REF]Pineau & Joly 1991;[START_REF] Pineau | Global and local approaches to fracture-transferability of laboralory tests results to component\[END_REF] sball be mentioned in particular.

'Then, contributions on new subjects emerged since the mid 1980's: shear band lo-calisatioo (Fressengeas & Molinari 1985. 1987;[START_REF] Benallal | Conditions de bifurcation a l' interieur et aux. frontims pour une dasse de ma~riaux non-standards[END_REF]. 1991[START_REF] Doghri | Investigation of localization due to damage in elasto.-plastic materials[END_REF]. stability analyses of systems of cracks [START_REF] Nguyen | Energy Jnethods in fracture mechanics, bifurcation and second variation[END_REF], scale effects in micro-macro analyses (Leroy & Molinari 1993; Lalai & Sab 1996). delamination of composites (Ailix & Ladeveze 1992, 1996;Allix ct al. 1995), asymptotic analyses (Geymonat et al. 1986; Leguillon & Sanchez-Palencia 1987). inverse problems and identification of cracks [START_REF] Bonnet | Regularization of the displacement and traction BJE for 30 ela~ todynamics using indi~ct methods[END_REF][START_REF] Andrieux | ldeatifJC&tion de fissure.~ planes pal' une donnee de bord unique: un proc6d6 direct de localisatiOA ct d'idenlification[END_REF][START_REF] Andrieux | Sur I' identification de fiS~rurc planes via le concept d'kart Ala n!ciprocite en 61af[END_REF], weak singularities of stresses at interfaces such as the "Epine' singularities (French word for the •thorn' like shaped function characterised by its boundedness, discontinuity and unbounded gradient) (Bui & Taheri 1989; Deperrois & Dang Van 1990; LeguiUon 1992; Santacreu 1994). experimental approaches in dynamic fracture (Bui & Maigre 1988; Bui et al. 1992; Maigre & Riu.el 1993) etc. In thesequel, some of these contributions will be discussed in more detail.

KEYWORDS: Linear elastic fracture mechanics, stres~ intensity factors. pathindependent integral, integral equation method, plasticity. damage mechanics. scale effects, inverse problems, fluid flow.

I CONSERVATION LAWS

The ]-integral, developed independendy by J.R. [START_REF] Rice | Mathematical analysis in lbe mechanics of fracture[END_REF] and G.P. [START_REF] Cherepanov | Cracks in solid~[END_REF], has motivated and stimulated work in static as well as in dynamic fracture mechanics, in small strain as well as in finite defonnation. The first paper on this subject is about the dual cooservatioo law in classical elasticity (dual to the E.~helby law). using the dual law of elasticity 2 = oU(a)loG instead of the primal one a= oW~ [START_REF] Bui | DuaJ palb-independent integral in the boundary value problems of cracks[END_REF]). The dual /-integral is then introduced as a path-independent integral equivalent to the J -integral for exoct solutions to the boundary value problems of cracked body. Considerations of dual approaches in computational mechanics provide bounds for the exact potential energy and. assuming some monotonic conditions for the meshes u.~ed in finite element approximal. ions. also provide approximate bounds for the derivative of the potential energy with re,;pect to the crack length.

The problem of separation and resolution of mixed modes by using the pathindependent integrals J 1 and J 11 is dealt with in [START_REF] Bui | A.~iated path-independent inregral~ lor separating mixed modes[END_REF]. The use of adjoint elastic fields, symmetrical and antisymmetric fields, leads to tbe same results [START_REF] Bonnet | Regularization of the displacement and traction BJE for 30 ela~ todynamics using indi~ct methods[END_REF].

The same idea can be found in the paper by H.D. Bui and H. Maigre ( 1988) who consider the adjoint dynamic fields for introducing the dynamic H-integral, which is a true path-independent integral in elastodynamics, while cla~sical extensions of the J-integral to take account of dynamic effects, as known in the literature, include area integral tenns over the crack tip region.

The analysis of fracture in thenno-elasticicy using adjoint elastic fields and adjoint thermal fields, leads to a new class of conservation laws which do not contain source terms. Thus, the true path-independent T-iotegrals in thermo-elasticity were developed (Bui &. ProiA 1984).

Area integrals over ring-shaped domains, the boundarie.~ of which consist of two non intersecting closed curves Giaa and Ge,t not containing the crack tip, were proposed for the static !-integral and tbe thermo-elastic T-integral. Independently. Ph. Destuynder and[START_REF] Destuynder | Sur une interpretation de l'integrc&le de Rke en tht1orie de Ia rupture ~lastique[END_REF]K. Ohtsuka (1981) introduced the La-grangian derivative of energy and defined the Gq-integral (the so-called q-method).

where q(x) is tbe velocity field of the cracked solid. parallel to the crack. tip propagation velocity inside the curve Gint' and equal to zero outside the curve Gexa• The ringshaped area integral corresponding to the Lagrangian version of the T-integral is called the A-integral (Bui 1993).

Stability analyses of crack propagation can be done either by the Eulerian second derivative of energy [START_REF] Nguyen | Energy Jnethods in fracture mechanics, bifurcation and second variation[END_REF] or by the Lagrangian second derivative of energy [START_REF] Suo | Sur une fonnulation matematique de Ia ~riv~ seconde de J'energie potentietle en theoric de Ia rupture fragile[END_REF]. In addition to the first field q(x). there is a second velocity field p(x) corresponding to a second crack, which justifies the naming of the so-called qp-method.

Stability analysis of crack propagation has been generalised to include propagation of damage in brittle elastic solids (Fedelich & Ehrlacher 1989;[START_REF] Francfort | Stable damage evolucion in a brittle continuum medium[END_REF]) and delamination of thin membranes (B6rest 1 987).

Conservation laws in finite strain analysis have been derived for elasticity. thermoelasticity and electroelasticity [START_REF] Maugin | The elecuoe1~tic cnergy~momentum rensor[END_REF][START_REF] Maugin | Pseudo-momentum and material forces in non-linear elastic-ity: Variational fonnulations and application to brittle &acture[END_REF][START_REF] Dascalu | Energy release rate and path-independent integrals in electro• elastic crack propagation[END_REF]. Also, the Lagrangian approach has been extended to take into account geometrical and physical non-linearities.

THEORY AND EXPERIMENT

New theoretical methods of analyses and new experimental techniques have to be developed in parallel for studying crack deviation in mixed modes, temperature rise during crack propagation, fatigue crack propagation under microscopic observations and dynamic fracture. Developments of new mathematical and numerical methods are necessary.

For example, the deviation of a crack in linear fracture mechanics can be studied by using an uncoupled theory, where the angle of deviation quite well predicted when the criterion sqq max is used. However, the uncoupled theory cannot predict the critical applied force. which is required for the determination of the true stressintensity factor at the tip of a kinked crack (in Mode 1). Therefore a coupled theory such as developed in [START_REF] Amestoy | Deviation infinit6imale d'une fissure dans une direction arbitrairc[END_REF](Amestoy et al. , 1985;;[START_REF] Leblond | Crack paths in plane situations. I. General form of the expan.~ion of the stress intensity factors./ru[END_REF]Amestoy & Leblond 1992;Leguillon 1993), using asymptotic matching of solutions, is necessary. A particular example is given by the exact solution for the infinitesimal kinked crack problem. The crack deviation angle. predicted by the criterion ~ = 0 of a pure local Mode I, does not differ significantly from the angle obtained from applying the criterion sqq max, with the stress field sq of the uncoupled theory. H the Irwin criterion holds true in any circumstance, then 1or the consistency reason, the critical value of applied force can be derived from the same formu1a k;= K~c-No such criterion can be found in the uncoupled theory.

Another example is given by the experimental determination of the dynamic stress-intensity factors associated with crack initiation. Based on theoretical analyses (Bui & Maigre 1988), a convolution equation is derived for the dynamic stressintensity factors K 1 (t) or K 11 (t). using only data available on the boundary of the cracked solid. the traction vector T(x, t) and the associated velocity vector v(x. t). For Mode I. the time-convolution equation is of the following fonn K.,(t)* K 1 (I) = H(T, v; T, v') where the prime stands for a known dynamic adjoint field in Mode I (velocity v'(x) and traction field T'(x) on the boundary) and the H-integral is explicitly known from the boundary data T and v.

Experiments on dynamic fracture of metals, polymers, and glass using the convolution method and the Hopkinson bar technique, which provide both the force and the velocity of the impact point. have been perfonned by Maigrc and Riucl (1993).

In the 1980's, attention has been paid to the problem of a dynamic crack acting as a moving heat source. Some experimental evidern:e hus been given by previous researchers, W cichen and Schooert ( 197 4 ). Based on tbermod ynam.ical considerations, Nguyen (1980), Bui et al. (1980, 1981), Ehrlacher (198.5), have proved that, in 20, the moving crack tip in an clastic material, due to the dissipation of energy at the crack tip, induces a logarithmic singularity of the temperature field. ExperimenlS were perfonned in thin sheets of steel, using a thennograpbical set-up. Indeed, we observed a hot point before the occurrence of diffuse plastic dissipation and before the development of damage, which is the result of shear band localisation throughout the thickness, as shown io Figure I. The final damage process induces the greatest dissipation rate.

Today, a large body of experimental work. is devoted to the study of microme~ chanisms of damage, fracture and fatigue in metals, as well as fracture and delami• nation in composites. [START_REF] Pinna | Mode II fatigue crack propagation in a maraging !lteel[END_REF] have concenlrated on Mode II fatigue micro-crack propagarion in metals, and based their work on observations made during testing in a microscope.

INTEGRAL EQUATION METHODS FOR CRACK PROBLEMS

Computational analyses of cracked bodies in 2D as well as in 3D have motivated work oo integral equations methods. The first papers on this topic by [START_REF] Bui | Application des potcntiels. elasliquc.<> ll'~de des fissures plane~t de forme arbitraire en milieu tridimensionnel[END_REF][START_REF] Bui | An integral equations method for solving the problem of plane crad of arbiuary sllape[END_REF] treated a plane crack of arbitrary shape in a 30 infinite elastic medium. A set of three integral equations was established: one uncoupled equation for the crack opening displacement of Mode I, and two coupled equations for Modes II and Ill. [START_REF] Bui | Etude expertmeotaJe de Ia dissipation dans Ia propagotion de fissure par thermo3raphie infrarougc[END_REF].

These equations establish the relationship between the stress intensity factor K 1 (s) and the intensity of tbe crack discontinuity near the crack front [u 3 (s. r)) = C(s}r 1 11.~ where s. r are the local co-ordinates near the crack front. It is shown that the relationship between C(s) and the stress intensity Kt<s) is tbe same as that obtained in 20 'plane strain' problems. This relationship holds true even for an edge crack. when the co-ordinate s approaches perpendicularly the free surface. [START_REF] Bui | Singulurite des contrain~ au voisinage d'une fissure angulcuse[END_REF].

lt is worth noting dlat the usual statement of a 2D state of plane stress at a free surface is incorrect. The proof can be derived directly from the equation of elasticity in plane stre.~s for which the singular field ~3 = O(r-112) at the crack tip. in plane stress.. violates the compatibility equations which imply that £:n(x 1 , x 2 ) is an affine function of the cG-Ordinates, ~3 = ax 1 + bx 2 + c.

Integral equations for arbitrarily curved cracks in 3D infinite space were established by [START_REF] Sladek | Three-dimensional curved crack in an elastic body[END_REF]. [START_REF] Levan | 1ntegral equations for 3D problem[END_REF], Bui and Bonnet (1987). [START_REF] Levan | Equations integrales pour fissures 3D. Contribution a Ia regulari~tion[END_REF], Bonnet and Bui (1993).

Many difficulties arise when dealing with numerical aspects of integral equations for crack problems. Firstly. the kernel is bypersinguJar (in the finite part sense of the integral), for wbicb regularisati.on techniques are needed [START_REF]Equations integrates de frontiere, derivation par rapport au domaine et approche Wergetique poor les solide[END_REF][START_REF] Bonnet | EqUQii()IU inlegraJe.t d ellments de jron1iere[END_REF]). The metbod of regularisaticn consists of subtracting and adding some singularity to the kernel. The additional term can be integrated exactly. while the remaining tenns. when combined. are weakly singular and can be computed with a better accuracy. Mathen\ati('.al~y. this procedure comes from the very definition of finite part integral in the distributional sense. Secondly, the structure of equations is asymmetric, for which variatiooal approaches are used [START_REF] Nedelec | lnaegr.csl equation.~ with non integrable kernels. lntegra1 equations and operator theory[END_REF]. Here, the symmetry is re• stored by a double boundary integration, which requires large computer memories. Variational approaches preserve both the symmetry and the regularity of integrals.

The integral equations method, in its variational form. has been used for solving 2D problems too [START_REF] Yin | Variational approach of displacement discontinuity method and application ro crack problems[END_REF] where the components of the crack dis• placement discontinuity were the unknown quantities.

Dynamic cracks problems. particularly diffraction or scattering of elastic waves in infinite space, are efficiently studied by means of integral equation methods [START_REF] Coney-Dumont | Simulation numerique de problemes de diffraction d•ondes par une fissure[END_REF]Ha Duong 1987;[START_REF] Becache | R.6!Wiution par une methode integraJe d'un probleme de diffraction d'ondes elastique uansitoire par une fissure . These Do[END_REF]Becache , 1993;;[START_REF] Bonnet | Regularization of the displacement and traction BJE for 30 ela~ todynamics using indi~ct methods[END_REF].

DUCTILE FRACfURE

This subject is motivated by potential applications in nuclear reactor technology. Considerable efforts have been paid in both mechanics of materials and metallurgy. particularly for studying stainle!OS steels and duplex stainless steels. Contrary to some initial hope, ductile fracture mechanics is not a simple generalisation of linear elastic fracture mechanics to non linear problems. nor can we state that the J-integral is the right parameeer in ductile fracture in every case. Experiments in ductile fracture sbow evidence of a specimen size effect aod theoreticaJ considerations indicate that not only the crack length, but also the geometry of the solid affect the value of J (the so-called resistance curve J-liD depends on the geometry of the solid). This structural dependence is due to physical non-linearities such as plasticity. viscoplasticity or damage and hence restricts the applicability of the resistance curve concept to very special cases {small scale yielding hypotheses). The interpretation of Jdoldt as the rate of energy dissipated inside the contour of the !-integral (taken in the elastic region) leads to the decomposition Jdo/dt =DR+ DP. with the dissipation rate by fracture DR at tbe crack tip and the plastic dissipation rate DP inside the same contour. In the classical model of n crack in plasticity. the term DR vanishes (Rice's Paradox). It is proved that DR vanishes for a large class of elasto-plastic materials, including hardening materials (Nguyen 1980). For a power law of plasticity. the paradox DR= 0 still holds true (Bui & Dang Van 1981).

In the classical approach. the path-independent J-integral is defined by means of lhe work density W(e) =fad£. It has no simple thermodynamical interpretation. unlike the case where the elastic strain energy density Wela.c;t = lf2 (~-e').E.( £4) is used for defining the path-dependent !-integral. Precisely, to obtain a non zero fracture dissipation rate, DR~ 0. the contour of the J-integral must be taken as the boundary of the damage process zone.

Analytical solutions in brittle damage of ela.~tic solids (elastic law a= E e when IEf < 0 and a • 0 otherwise) provided by [START_REF] Bui | Etude expertmeotaJe de Ia dissipation dans Ia propagotion de fissure par thermo3raphie infrarougc[END_REF], lead to DR~ 0 and ~how that the damage process zone has some thickness. In [START_REF] Bui | Etude expertmeotaJe de Ia dissipation dans Ia propagotion de fissure par thermo3raphie infrarougc[END_REF] ), a particular solution in damage plasticity is given for the stationary case of damage in Mode 3.

The interl'ace between damage and plastic zones is the cycloid N (N 1 = h(cos2cp + 1)/Jt, N 2 = h(sin2t + 2cp)l'1t). The elastic-plastic interfa~ depending on the remote stress

intensity factor K 3 , is the cycloid M (M 1 = K~ (cos2' + 1 )121t't 0 2 • M 2 = Kj sin2cp/2Jrto 2 + 2hfPht).
For classical Mode 3 in plasticity, the cycloid M becomes the cla~ical Rice circle and the damage zone reduces to the crack h = 0, Figure 2.

For a more general loading in damage plasticity, unfortunately there is no exact solution. Numerical works are only available, based on the introduction of damage evolution law on a microscopic scale. The Rousselier (1981) criterion appears to be the best one for predicting the growth rate of cavities in ductile materials like stain-Jess steels of tbe French nuclear vessels. It is an extension of the Gurson law. Finite elements solutions in plasticity. including ductile fracture damage, can be found in many works (D'Escalba & Devaux 1979;[START_REF] Mudry | Etude de Ia rupture ductile et de Ia rupture par clivage d'acier faiblement aiU~[END_REF][START_REF] Beremin | A local crkerion for cleavage f01eture of a nucle8f preSSllre vessel steel[END_REF][START_REF] Wadier | The theta method applied to analysis of 3D elastoplastic structures[END_REF][START_REF] Debruyne | An arbitrary Lagrangian Description or 2D and 30 Cracked Structures[END_REF]. The macroscopic behaviour of solids containing cavities or cracks periodically distributed at microscopic scale can be analysed by the homogenisation method [START_REF] Andrieux | Un moc.leJe de mat~riau microfissure avec frottement[END_REF][START_REF] Suquet | Ptasticite et homogeneisaljon[END_REF][START_REF] Marigo | Modelling of brittle and fatigue damage for elastic material by growth of micro[END_REF]. A continuous approach to ductile fractu~ using damage mechanics. is provided by (Dragoo & Mroz 1979; Lemaitre & Chaboche 1985). The strong dependence of local damage criterion on the materials. and perhaps the dependence on the loading conditions too. prevents its applicability to different cases. However, a methodology exists for analysing ductile fracture of components. It consists of global and local approaches to ductile fracture which are based on the transferability of laboratory tests results. on specific materials. to specific compooo nents. with the use of numerical computations in non-linear mechanics [START_REF] Pineau | Global and local approaches to fracture-transferability of laboralory tests results to component\[END_REF]. This methodology justifies the development of French Computer Codes for industtial applications. such as the CASTEM 2000 of the French Atomic Centre. the Aster Code of Electricite de France (Mialon 1996) and the Systus Code of Framatome. and requires also sophisticated equipment when laboratory testing metals.

DAMAGE MECHANICS

Damage in materials occurs due to initiation of micro-cracks and microcavities in highly defonned zones. e.g. at the crack tip. at grain boundaries, at interfaces between inclusions and matrix of composites. Oo the meso-scale of several grains, the homogenised mechanical propenies are characterised by a decrea.~ of the elastic stiffness E, softening or a decrease of the yield strength, for increasing strain. A continuum theory of damage can be developed within the framework of the theory of plasticity, using the concept of generalised standard materials introduced by B. Halphen and Q.S. [START_REF] Hajphen | Surles materiaux standardli generalis~[END_REF]. Internal parameters. [START_REF] Geymonat | Singularites dans les joints collH[END_REF];. are introduced into the constitutive laws such that their evolution rates dm;fdt = ~(a. (!);~ e") and the rate of plastic deformation df!'ldt = ').J(o, 00;, f!'), with the same plastic multiplier :A. > 0 and the familiar yield function f(o, m;, £1'), are concomitant. Some assumptions are necessary for specifying the forms of function.~ g( o. c»;, £") and A((!);) (Dragon & Mroz 1979). The simplest model [START_REF] Lemaitre | Aspect pheno~nologique de Ia rupture parendommagement[END_REF] allows the detennination of damage by stiffness measurement. The damage parameter is a scalar 0 < c» < 1 and the dependence of the elastic stiffness on ro is assumed to be of the simple form A( m) = (1 -m)Ao.

In some contributions, internal damage parameters are con,;idered as microscopic variables which corre.~pond to new stress fields satisfying additional equilibrium equations [START_REF] Fremond | Damage, gradient of damage and principle of virtual power[END_REF]. Such models using different size scales are aimed at the stabilisation of localisation phenomena.

A simpJer model of elastic-brittle damage is provided by the following damage criterion: m = 0 if f. = I G 1a 0 < 0 and co = 1 if f. = I a 1a 0 = 0. 1bis allows for analytic solutions to damage propagation problems. There exist two distinct zones in the solid. the undamaged material 20ne ((I) = 0), which is elastic, and the zone associated with complete damage, (m = 1). Since the boundary between these zones is unknown, we are let to free boundary value problems. in quasi-statics as well as in dynamics. For example, the steady state dynamic propagation of damage in an infinite medium, in Mode III fracture, has been solved analytically [START_REF] Bui | Propagation of damage in elastic and plastic solid~[END_REF].

Several papers are devoted to the propagation of discontinuities in solids due to a sudden change of properties [START_REF] Stolz | Sur le probleme d•~volution thcnnomecanique de" solides a changement brutal de caracteristiques[END_REF] or due to phase transfonnations [START_REF] Buisson | Comportement global as.'\ocie ault mouvements d•interfaces entre variantes de martensites[END_REF][START_REF] Chateau | Modclisation m6:aniquc d'un milieu continu polyphasique avec 6:hange de mari~re entre constituants[END_REF]. Thennodynamical analysis shows evidence of internal entropy production, concentrated on the moving damage front which behaves like a heat source [START_REF] Stolz | Approche energ6tique de Ia propa&ation dynamique de discontinui~ mkaaiques[END_REF]. 1be familiar energy release rate in fracture mechanics can be generalised to damage propagation or to propagation of discontinuities in tenns of the jumps of free energy and velocity gradients aero~\ the damage front.

CRACK INITIATION

Prior to rupture, localisalion of defonnation and damage by instabilities may occur and these phenomena can be used to define the onset of frc~cture. Jn a general context. the deformation und damage may be localised in the bulk of a structure, at the boundary or at interfaces. For a given loading process. the localisation conditions. related and based on earlier works of Hadamard (Mandel 1963), define the instant of occurrence of these phenomena but also the orientation of an eventual crack associated with a ma"imwn tensile stress plane. Depending on the constitutive behaviour, the type of boundary conditions applied to the structure and the type of interfacial behaviour, they also allow to predict rupture to start inside, at the free surface or at an interface [START_REF] Benallal | Conditions de bifurcation a l' interieur et aux. frontims pour une dasse de ma~riaux non-standards[END_REF](Benallal et al. , 1991)).

Cracks frequently initiate at the interface of two materials (Desmorat 1996). at the junction of dissimilar elastic materials and the free surface [START_REF] Zwiers | On lhe logarithmic singularily of free-edge stress in laminaced composites under unifonn extension[END_REF]. Most developments have been devoted to such interface singularities. However. there are new kinds of singularities. related to the interface between dissimilar thennal or plastic properties and al~ at discontinuity points of the loading di5tribution. At the intersection line between the free surface and the interface where the thermal strain eT 22 or tbc plastic sttain P!' 22 are discontinuous, there are weak singularities. (Bui & Taheri 1989) coined the term 'Epine, singularity, called after the French word for the 'thorn' like shaped function. The •Epine' singularities of the stress components o 11 are characterised by their boundedness, discontinuity and unbounded gradient. 1bey occur when tbe interface normally intersects the free surface [START_REF] Deperrois | Inclu&ions de surface et ~iogu1arit6 ~pine[END_REF]Leguillon 1992;Santacreu 1994). The same kind of singularities arise at points of discontinuous normal stress applied to the boundary of elastic solids [START_REF] Barton | The circular cylinder with a band of unifonn pressure[END_REF]. They are exactly the local effects solutions in plates or shells considered as 3D bodies, defined a., the difference between a 3D solution and the Kirchhoff plate or shell solution. Thermal striping, crack initiation at the interface of a weld can be explained by the 'Epine• nonnal stress singularity. For example at a plastic interface, the jump relations at the free surface are of the form [o 1 1 ) = 2p. [d 22 ] in plane strain [START_REF] Bui | La ~ingularitc •Epine• dan.~ le5 bi-mat~riaux. en thermo-~lastoplasticitt[END_REF]. [0' 11 ]

= 2f,L([ff 22 ] + 'U[£1' 33 ])/(1-'U) in 3D (Santacreu 1994). If
the thennal and plastic interfaces are inclined (angle a between the interface and the normal to the free surface). the ~Epine' singularity exists with the factor cos2a, but it is masked by a stronger logarithmic singularity with the factor sin2a. The 'Epine' singularity appears also in the asymptotic analysis of an adhesive joint between two elastic solids, intersccting the free surface, as the joint thickness tends to zero [START_REF] Geymonat | Singularites dans les joints collH[END_REF].

Studies of delamination in composites have important applications in material science. Damage modelling of an elementary ply in laminates composites (Allix & Ladeveze 1992, 1996, Allix et al. 1995) consists of the characterisation of the layer by homogeneous internal parameters, describing both, transverse cracking of the matrix and fiber rupture~ and of the identification of interface properties by cla.C\sical experiments (DCB specimens, traction on plates with holes).

Crack initiation originates at points near the boundary of solids or inside subsurface regions of high residual stresses. when the solids are subjected to plastic strain. Classical finite elements methods in non-linear computational mechanics provide efficient mean~ for analysing monotonic loadings. However, these methods cannot be applied to shake-down analyses, plastic acconunodation and ratchetting studies. Significant progress has recently been made (Dang [START_REF] Van | Steady-state flow in classical elastoplasticity; applications to repeated rulling and sliding &[END_REF]; Dang [START_REF] Van | Elastoplastic analysis of repeated moving contact. Application to railways damage phenomena[END_REF] for the computations of re.~idual stresses associated with moving contact loads. repetitive contact loads. cyclic surface loads etc .• [START_REF] Fouvry | Relation between variou.'l Fracture Transition Tempe111tures and the K. 1[END_REF][START_REF] Stalin-Muller | La mod~Jisation de,c; phenomenes de l'usurc, eoouche d'un modele quandtatif[END_REF]. These new computational methods have very important applications in fatigue and wear of materials and structures.

Computational work is required for detennining the macroscopic cyclic stress and strain, as well as the cumulative plastic strain. For a given load path. it is particularly important to determine the endurance limit. This is provided by the well-known Dang Van criterion in multiaxial loadings (1973) which is now widely used in French industrial applications. The Dang Van criterion expressed in the fonn ft:t, G,J = 0, where tis the maximum shear stress and a,. the hydrostatic stress. is based on a micro-macro analysis of polycrystals with various localisation fonnulae (Lin-Taylor assumption. self-consistent model etc.). Microscopic models developed for metal plasticity [START_REF] Zarka | Contribution a Ia theorie de Ia pla[END_REF][START_REF] Zaoui | Etude de l'iAfluem."e de Ia desorientation propre des grains sur le comportement 'Yisco-plastique de metaux polycristaJiins (CFC)[END_REF]) are used only to justify the theoretical form of the fatigue law f(t. om) = 0. The constants in the fatigue laws may be determined directly by experiments. Recent work in fatigue damage of metals using 'double scale models' renew the ideas of Dang Van's model developed in the 1970's.

SCALE EFFECfS IN FRACfURE MECHANICS

Fracture of quasi-brittle materials such as concrete, rock and ceramics presenlS a very complex feature because of the heterogeneity at microscopic scale and of different properties of the constituents (matrix~ inclusions. cavities, etc.). Experiments on a macroscopic scale show evidence of softening behaviour and dependence of the fracture stress, o-•. on the specimen si1.e, L.

There are two kinds of scale effects. The first one corresponds to the volume scale effect of a body with si7..e L, having cnlcks at the microscopic scale. The fraclure stress is of the fonn d-AL-a (where A and <X> 0 are characteristics of the material). This scale effect can be explained by the heterogeneity at the microscopic scale.

1be second one. called structural scale effect. is due to the non-linearity of macroscopic behaviour a~ already mentioned in ductile fracture of a body with a macroscopic crack.

Tbe simplest models of non•lineac fracture mechanics are associated with small scale yielding hypotheses. They are based on classical elasto-plasticity on the macroscopic scale (resistance curve model. Hillerborg's fictitious crack model). More complex models take into account the evolution of damage in the vicinity of the process zone which exhibits softening behaviour. To circumvent difficulties arising in numerical instabilities. such as the mesh size dependence or the localisation in one finite clement etc .• one introduces a certain characteristic width. It results in localisa-tion limiters as can be provided by higher-order gradients theories to account for micro-structural effects [START_REF] Leroy | Spatial pauems and size effects in shear zones: a hyperela~tic model with higher-order gradients[END_REF] and a non-local damage theory (Pijaudier-Cabot & Bazaot 1987; Mazars et al 1991). These models can be analysed either numerically by means of finite element methods or, in some special cases (hyperelastic models), analytically (Leroy & Molinari 1993).

Stochastic models were developed in order to take into account the heterogeneity of materials and the random chamcter of mechanical constants on a micro-scale. H.J. Hermann ct al. ( 1973) describe micro-cracking and the process of fracturing in heterogeneous solids by a lattice model of beams. Stochastic finite methods have been developed [START_REF] Rossi | The usefulness of statistical models to describe damage fracture in concrete[END_REF] for modelling the volume scale effect observed in experiments, and to study the influence o f heterogeneities. In the latter model, micro-cracking is the result of decohesion between interface elements at some critical stress. Both, volume and structural scale effects have been accounted for in the stochastic non-local model of quasi-brittle fracture [START_REF] Lalai | A stochastic non-local damage model[END_REF], as shown in Figure 3. Three size scales have been introduced in the so-called m3 model (micromeso-macro). On the microf;Cale of the grain size d. the non-local clastic-brittle law is assumed. On the meso-scale of the representative volume of size L'>> d, the assumption of heterogeneity of the grains and randomness of critical stresses arc assumed. The sizes d and L' are small in comparison with the size L of the structure (macro-scale) L >> L">> d. The model by Lalai and Sab is based on Weibull's h'talisticallaw (with constant 'm') and thus depends on tbree constants (d, V , m).

Beyond some critical size C, the volume size effect can no longer be noticed. This corresponds to the minimum volume size at the me.~scale for whicb a macrohomogeneity throughout the body is observed.

[t&pe 2

INVERSE PROBLEMS IN FRACTI.JRE MECHANICS

lnve~e problems in fracture mechanics are essentially concerned with nondesuuctive-tcsting of materials and structures for the detection of cracks, cavities or inclusions in solids [START_REF] Cielo | Thermographic non-desuuctive evaluation of industrial materials and structures[END_REF]: Zorgat et al. 1991 ). Classical approaches to solve such an inverse problem include ultrasonics. Eddy current metbod~ X-ray or r-ray tomography. and thermography. Recently, some theoretical methods for crack detection have been worked out for inverse boundary value problems.

Inverse boundary value problems address the identification of cracks in three dimensions, where, similar to situations in potential theory as well as in the theory of elasticity. a surplus of data along the boundary is available. Such problems have been widely studied in 20. We shall report here some new results in 3D elasticity.

In elasticity theory, the data consist of the pairs {u, T(uJ}. which denote the displacement and the traction vector fields on the exterior boundary ~ respectively. These pairs correspond to some plane crack to be detennined. In the classical solution method, a family of solutions. v, of the elasticity equation for the uncracked body is considered~ and the functional RG(v): =l: IJT(v)•u -T(u) •v}dS is defined, with the integral taken over the exterior boundary of the cracked solid, ~. If tbere is no crack, the Betti reciprocal theorem states that RG(v J = 0 . Therefore the functional RG(v) provides a mechanical measure of the 'signature' of the crack and is called 'Reciprocily gap' functional. The RG(v) functional being known from the pair of data, can be expressed in tenns of the displacement discontinuity [u] as follows Jp{ T(v) • (u]dS = RG(v), where the integral is taken over the unknown crack surface F and the right hand side RG(v) is the known linear functional. This equation, for appropriate choices of basis functions v, delermines completely the crack plane F (normal vector Nand position in 3D space) [START_REF] Andrieux | ldeatifJC&tion de fissure.~ planes pal' une donnee de bord unique: un proc6d6 direct de localisatiOA ct d'idenlification[END_REF]. Recently. the complete identification of the crack bas been obtained by determining the support of the displacement discontinuity, with the fields, v, chosen so a.s to obtain the Fourier transform of (u] (Andrieux et aJ. 1997). Mathematically, the crack identification problem is solved explicitly by determining the projections of the singularity of u. or its discontinuity [u], on appropriate basis vector fields v of an Hilbert space.

FLUID-FILLED CRACKS

Fluid-filled cracks have been tbe subject of many papers. Previous studies of a fluidfilled crack in an elastic medium have been based on approximate global equations of continuity and motjon. and on an assumption of one-dimensional flow in the crack. Such investigations have led to the existence of singularities in the fluid pressure and hence the conclusion that a crack cannot be completely fiJied with fluid. The lal\t assumption can be later introduced in another natural way by considering the surface force (no complete filling is pos.\ible when the radius of curvature of the crack tip is smaller than some water molecules). However, at a macroscopic scale. without considering the surface force. the assumption of partial filling of fluid is not necessary if we incroduce the two-dimensional viscous flow near the crack tip. Such an extension from one-dimensional to two-dimensional fluid flow bas been given [START_REF] Bui | A reexamination of the pressure at the lip of a fluid-filled cra<:k, 1111[END_REF]. It is shown that the exact equations for 20 fluid-filled cracks do not give rise to a singular pressure.

The case of crack propagation with a fluid-filled crack is most interesting, from both, a physical and mechanical point of view. As tbe fluid is assumed to be viscous, the physical adherence condition at the crack surface is not compatible with the assumption of a fluid completely filling a Griffith crack, because of the discontinuhy and singularity of the vertical components of velocity, u 2 (0, ± 0) = ± oo, of the solid points lying at the moving crack tip D (x 1 • 0). The elastic solid point at the crack tip D breaks up into two parts, moving in opposite directions. Therefore. in order to avoid flow singularities the viscous fluid is assumed to partially occupy the crack zone, with a convex (or concave) meniscus in equihorium with a 1\urface tension force. A model of steady state flow of a viscous fluid-filled crack has been provided (Bui et at. 1991 ). Using the N3S Codes (Navier-Stokes computer codes of Electricite de France), the coupled problem of viscous fluid flow in the moving frame at the crack tip has been solved numerically as shown in Figure 4.

The case of a stationary crack has been studied by H.D. [START_REF] Bui | The interaction between the Griffith crack and fluid: Theory of Rebbinder's ef-Fects[END_REF] who gives the analytical solution of the problem of a partially fluid-filled crack in an infinite medium, with the unifonn trctction at infinity. It is shown that for a convex meniscus, a stress intensity factor exists even when tbe tradion at infinity vanishes. Consequently. the existence of a 'capiltary• stress intensity factor lowers the applied effective stress inten.~ity factor at rupture. This model provides u simple explanation of the well known Rehbinder effect (Bui 1996). The toughness of wet material is generally smaller than that of a dry material, as shown experimentally (Rossi & Boulay 1990).

Most interpretations involre the chemical change of material properties by water or moisture. Here. the interpretation of the Rehbinder effect is ba.<ied on a purely mechanical argument and on the assumption of capillary force.

Viscous fluid flow through a longitudinal crack in a pressurised cylindrical shell has been studied with applications to the LBB problem (Leak-before-break) [START_REF] Bui | Modeles analytiqucr. de debit de fuilc dans une coque cylindrique n~uree Joogi• 1udinafement[END_REF]. lt is shown that there are two kinds of fluid ftow corresponding to the values of Folias• parameter A = a(12(1 -v 2 ))14 (hRr'h~ where •2a• represents the crack length, R is the mean radius of the cylindrical shell, h the thickness of the shell. vis Poisson ratio. for l> 4. 78, the semi-axis of the elliptical crack opening displacement decreases from the inner radius to the outer radius of the tube (convergent flow) and for A-<4.78, the semi-axis of the elliptical crack opening increases (divergent flow). Theoretical predictions of the fluid tlow rate agrees quite well with experiments on c surface tension force along AB, no fluid-solid contact along AB; adherence along BC wilb relative velocity -V (Bui etal. 1991 ).

pressuriscd cracked tubes. The LBB studies of a longitudinally cracked cylinder are based on the use of crack singularity analyses provided by D. Bergez and D. Radenkovic ( 1972). The singularities of cracked plates or shells are characterised by 5 coefficients: k 1 , lc 2 • k 3 (the clas.~ical three modes in the middle surface of the shell) and the rates of change through the thickness k' I ;::; aK.(xl)taxJ, "2 = oK2(x'J)I()x3 of Modes I and IT stress-inten.~ity factors.

Figure I .

 I Figure I. Thermographic observarion of crack propagation in a [hin &heet of steel[START_REF] Bui | Etude expertmeotaJe de Ia dissipation dans Ia propagotion de fissure par thermo3raphie infrarougc[END_REF]).

  Figure 2. (a) Damage plasticity; (b) Damage and plastic zones; (c) 1be limiting case of a crack as ER-+-(Rice's circular pla~tic zone) (Bui & Ebrlacher 1981).

Figwe 3 .

 3 Figwe 3. Numerical simulation of quasi-briUle fracture wilh 3600 S-Nodes square elements, with raadom critical stress in lhe sides and diagonals of elemeats. Steps 1. 2. 3: d itluse dama~ with increasing toad. Steps 4 and S: localisation in a crack-like ~one. with decreasing load[START_REF] Lalai | A stochastic non-local damage model[END_REF].

Figure 4 .

 4 Figure 4 . Steady swe viscous fluid flow in the movine frame at the crack tip D. wiah velocity V;


Maxime Rondeau

L Isnard

R Arès

Point-to-source path tracing Monte Carlo to compute the Clausing and distribution functions in high-vacuum systems

This paper presents an approach to simulate the free molecular flow in vacuum systems by using a Monte Carlo method for solving the Boltzmann particle transport equation with no intermolecular collisions. Sometimes referred to as a point-to-source Monte Carlo path tracing in image rendering, in this paper the name is borrowed from the thermal radiation and heat transfer field, reverse Monte Carlo path tracing. It is shown that this method provides better accuracy and stability when computing the Clausing function when compared to the standard test particle Monte Carlo method used for free molecular flow. The Clausing function leads to the distribution function of positions and velocities from which the particle density map, pressure gradient, energy flux, and other local quantities can be computed. Using reverse path tracing, the particle concentration in a conical segment is computed, and the maximal flow input is determined by calculating the mean free path at the maximum density position.

I. INTRODUCTION

Many applications, such as spectroscopy, microscopy, epitaxy, lithography, distillation, arc remelting, and particle accelerators, benefit from or require high vacuum when operating. In those environments knowledge of the particle distribution anywhere on the walls may, for example, help predict coating areas due to gases being injected or particles etching off surfaces during processes, optimize the location of pumps to better maintain vacuum, and therefore improve the overall design of systems. Simulating the particle distribution in the whole or key sections of such systems requires the solution to the collisionless Boltzmann transport equation. However, solving this equation is a computationally expensive task in environments with complex geometries. This is why efficient solving methods are interesting when it comes to optimizing the design of vacuum systems.

Such systems operate in a regime called free molecular flow which was first clearly demonstrated by Fleming [START_REF] Fleming | On a phenomenon of molecular radiation in incandescence lamps[END_REF] when he observed a shadowing effect made by a horseshoe carbon wire blocking the evaporation of a thin copper film in a light bulb. Such a phenomenon can be explained when the mean free path (λ) of molecules is much higher than the characteristic dimensions of the system (L), meaning the molecules radiate from the copper layer and directly hit the bulb. The Knudsen number (K n ), the ratio between λ and L, is thus very large (K n ≫ 1) for highly rarefied gases. A more mathematical description of this flow regime was introduced by Knudsen [START_REF] Knudsen | The laws of molecular flow and of inner friction flow of gases through tubes[END_REF], who studied the flow of gases through narrow tubes and estimated the energy throughput to be proportional to the pressure gradient, with a proportionality coefficient being purely geometrical. Only the shape introduces flow resistance as expected for a regime with a large K n . However, it was Clausing [START_REF] Clausing | The flow of highly rarefied gases through tubes of arbitrary length[END_REF] who developed an integral equation, with a line-of-sight approach, to compute flow resistance induced by a cylinder of length L and radius R. Solving the integral equation actually gives access to the Clausing function, the probability w(x) for a molecule that leaves any position x along the length of a tube to exit at the end:

w(x) = w rs (L -x) + L 0 dξw(ξ )w rr (ξ -x), (1)
where w rs (Lx) is the probability for a molecule to directly leave a ring section of the tube at position x and exit at position L through a disk. On the other hand w rr (ξx) is the probability for a molecule to go directly from a small ring section at position x to another one at position ξ . Solving Eq. ( 1) therefore gives access to w(x), which helps to compute the Clausing factor (W), defined as the probability for a molecule that diffusely enters the tube to exit at the other end, which is finally used to compute the conductance of pipes. The flow conductance helps to estimate pumping times for vessels using pipes of various aspect ratios. The Clausing function w(x) also gives access to the distribution function and therefore to the local density, mean free path, mass transport, and other quantities that help characterizing vacuum systems. Using local particle conservation, one can derive the transport equation for free molecular flow, which can be seen as a generalized Clausing equation (not only for pipes). Solving this equation is no easy task. Simple geometries, such as shapes with axial symmetry, reduces the Clausing equation to a one-dimensional integral, allowing numerical integration to be carried out with the use of orthogonal polynomials. This effectively transforms the integral into a matrix equation. Such a method was used by Iczkowski [START_REF] Iczkowski | Effusion of gases through conical orifices[END_REF] for the effusion through conical orifices. When no symmetry can be exploited it is usually necessary to use numerical methods to carry out the integrals. Combining the numerical integral with an iteration process and a cutoff procedure, the solution to the Clausing equation was computed by various authors for applications that use gas injection or evaporation, molecular beam epitaxy by Curless [START_REF] Curless | Molecular beam epitaxy beam flux modeling[END_REF], chemical beam epitaxy by Isnard [START_REF] Isnard | Theoretical and experimental molecular beam angular distribution studies for gas injection in ultra-high vacuum[END_REF], and physical vapor deposition in small trenches by Cale [7]. This last computing method is very fast, but is harder to implement when the number of objects in the system increases. Another approach is to solve the equation with Markov chains to sample the integral. Again, with surface discretization, this Monte Carlo approach was used by Lobo and Gomez-Goni [START_REF] Lobo | Comparison between monte carlo and analytical calculation of Clausing functions of cylindrical and conical tubes[END_REF] for pipes and conical shapes. Essentially, most of the solutions require defining small areas or patches that make up the objects, a tedious task when faced with complex systems, and also a source of small error. Moreover, in this last Monte Carlo method, smaller area elements scale with greater computing time.

The Clausing equation can be linked to the particle flux g, defined as the number of impinging or emitted particles per unit time per unit area, on any surface in the system, and is in fact nothing more than the collisionless Boltzmann integral transport equation with diffuse surfaces, giving the so-called radiosity equation g = g 0 + Tg. In this equation g 0 is the distribution of particles on the surfaces of the system with no surface collision prior and T is the transport operator of the system. In the case of free molecular flow the transport operator is simply a matrix of geometric factors. Since the flow is described by a standard particle transport equation, many solving methods are known for heat, light, and neutron transport. One of the methods to solve the transport equation is to use the adjoint solution from the adjoint transport equation [START_REF] Carter | Particle Transport Simulation with the Monte Carlo Method[END_REF][START_REF] Lux | Monte Carlo Particle Transport Methods: Neutron and Photon Calculations[END_REF]. However, in the case of free molecular flow, the transport equation is self-adjoint, and the adjoint solution should be the solution to the transport equation.

This work is inspired from a Monte Carlo path-tracing sampling method used in optics [START_REF] Kajiya | The rendering equation[END_REF], neutron diffusion [START_REF] Lux | Monte Carlo Particle Transport Methods: Neutron and Photon Calculations[END_REF], and radiative heat transfer [START_REF] Howell | Thermal Radiation Heat Transfer[END_REF], which is applied to the free molecular flow, similarly done by Malhotra [START_REF] Malhotra | High Knudsen number physical vapor deposition: Predicting deposition rates and uniformity[END_REF]. The procedure proposes to begin the particle path at the point where one wants to compute the impinging flux and keep emitting from this location in a diffuse way until one hits an exit or a particle source area. This procedure is called reverse path tracing from now on, a name borrowed from the field of thermal radiation heat transport [START_REF] Howell | Thermal Radiation Heat Transfer[END_REF]. This is in contrast with the method used by Gomez-Goni [START_REF] Lobo | Comparison between monte carlo and analytical calculation of Clausing functions of cylindrical and conical tubes[END_REF] mentioned previously. In this Monte Carlo approach, molecules originate at the particle source and are emitted diffusely until they exit or return to the source. To compute the impinging flux one needs to count the number of collisions made in a small area of length ε centered around the desired location. Both processes are illustrated in Fig. 1. This latter FIG. 1. Illustration of the processes involved in the two path-tracing methods. The particle source is a uniformly distributed diffuse emitter represented by the bottom cap of the cylinder, the top cap surface is called the exit area, and the body of the cylinder is diffuse re-emitting surfaces called walls. Particles are generated at the black spot bounce of the walls, and their trajectories are illustrated in two colors. Line paths do not contribute to the impinging flux whereas dashed paths do. In the direct path example, the particles are captured in a small collecting area of length ε.

method will be called forward or direct path tracing from now on, but is also known as a test particle Monte Carlo method [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF][START_REF] Sharipov | Data on internal rarefied gas flows[END_REF].

The objective is to demonstrate that when reverse path tracing is used properly it offers greater accuracy and immensely more precision to calculate any quantity when compared to direct path tracing for a given computation time. This is especially true when the wanted solution is only for a subset of points in a large system. Furthermore, being a path tracer, it offers the same statistical information on any desired quantity that direct path tracing would. For example, it is possible to compute the transmission probability as in direct path tracing; see Table I.

Moreover, it does not require the use of small areas to collect particles while maintaining the advantage of using Boolean operations between simple geometrical shapes, such as cones, pipes, planes, and spheres, to construct complex systems. However, it does suffer from the same shortcoming of any Monte Carlo method, slow convergence (when compared to other numerical integration methods) proportional to the inverse square root of the number of generated paths. TABLE I. Clausing factors (W ) expressing the probability for a molecule that enters a cylinder of aspect ratio length over radius L/R to exit at the other end. The integer N represents the number of paths generated to obtain the transmission probability in this work, which is compared to Ref. [START_REF] Li | A modified numerical method for the accurate calculation of molecular flow transmission probabilities and density distributions of cylindrical tubes[END_REF].

L/R N W (this work) W ([ 16

II. THEORY AND DEFINITIONS

A. The transport equation applied to highly rarefied gases

A free molecular flow system in a steady state, with isothermal walls and sources, is composed of three surface types:

(1) Particle source areas, or inlets, are characterized by a cosine law of emission at each point of the surface, with an emission rate that may change from point to point. However, the angular distribution of emission rates on the source area remains constant. This means that incoming particles from other surfaces hitting the source area are recycled so that the emission pattern stays unchanged.

(2) Walls re-emit all incident particles with a cosine law assuming there is no sorption and particle migration on such areas. Even though the cosine law of emission does not apply [START_REF] Davies | The failure of theory to predict the density distribution of gas flowing through a tube under free molecular conditions[END_REF][START_REF] Sharipov | Influence of gas-surface interaction on gaseous transmission probability through conical and spherical ducts[END_REF] for all systems, this paper treats particles that are thermally accommodated by the system walls. A previous paper has shown that such conditions are experimentally attainable [START_REF] Isnard | Theoretical and experimental molecular beam angular distribution studies for gas injection in ultra-high vacuum[END_REF].

(3) Particle exit areas, or outlets, do not reemit anything, particles hitting such a surface are lost outside the system.

Let g(r) be the impinging or emitted flux of particles at point r in a system with a steady flow. Then it is possible to compute the flux g(r):

g(r) = ν S i dS i p(r i → r) + S ′ g(r ′ ) dS ′ p(r ′ → r), (2)
where the integral on S ′ spans over all the walls, and ν is a uniform (uniformity is not an essential criteria) and constant source of flux, S i .T h et e rm[ g(r ′ ) dS ′ ]p(r ′ → r) dS can be read as follows: the number of particles transferred to an area dS located at r that originated from an area element dS ′ around position r ′ with an emitted flux of g(r ′ ). Since g(r ′ ) dS ′ is an emission rate, the number of emitted particles per second, then p(r ′ → r)dS must be the fraction exchanged between the two areas, with Eq. ( 3) giving the exchange fraction per unit area. All area elements are scalar quantities, the vectorial components of the surfaces being incorporated,

p(a → b) = (r b -r a ) • n a (r b -r a ) • n b π |r b -r a | 4 , (3)
where (r b -r a ) is the distance vector between the two positions, and n b , n a are unit vectors normal to the area elements, always oriented so that the scalar product with the distance vector is positive.

Just by looking at Eq. ( 3), one can convince oneself that the exchange fraction per unit area is reversible:

p(a → b) = p(b → a). (4)
Equation ( 4) also works for specular reflection and may hold for other scattering laws. Therefore, reverse path tracing in free molecular flow is valid for gas-surface interaction models that respect Eq. ( 4). This paper concentrates on diffuse scattering expressed by Eq. ( 3), and therefore the particle transport Eq. ( 2) can be related to the Clausing equation [START_REF] Clausing | The flow of highly rarefied gases through tubes of arbitrary length[END_REF]. Rewriting Eq. ( 2) using a matrix-vector notation,

g = g 0 + Tg, (5)
where g 0 is the first integral on the right-hand side of Eq. ( 2) and Tg is the second one, then the solution is (see Refs. [START_REF] Lux | Monte Carlo Particle Transport Methods: Neutron and Photon Calculations[END_REF][START_REF] Courant | Methods of Mathematical Physics[END_REF][START_REF] Arfken | Mathematical Methods for Physicists International Student Edition[END_REF] for more details)

g = (1 -T) -1 g 0 . (6)
Since the fraction of particles that are transferred between any two small sections of emitting or reemitting walls is always inferior to one, given by the integral operator T when all possible section pairs are considered, Eq. ( 6) can be developed into a sum using the Neumann expansion [START_REF] Arfken | Mathematical Methods for Physicists International Student Edition[END_REF]: Basically, Eq. ( 7) means that in a stationary flow, the total impinging flux at a certain position on any surface of the system is a direct transfer of particles with 0 collisions [g 0 (r) = S i dS i νp(r i → r)] plus the transfer with one collision [Tg 0 (r) = S ′ S i (νdS i )p(r i → r ′ ) dS ′ p(r ′ → r)] and so on. Equation ( 5) can be imaged as shown in Fig. 2, where the thick arrow represents the total flux at one point and the thinner ones represent g 0 , and each collision on a wall is an integral over the area of the re-emitting walls in the system.

g = g 0 + Tg 0 + TTg 0 + TTTg 0 +••• . (7)
Since the exchange fraction per unit area is reversible, then looking closely at g 0 and using Eq. ( 4) one finds that

g 0 (r) ν = S i dS i p(r i ← r)( 8 )
is actually the fraction of particles that are emitted at r transferred directly on the source area S i .

The second term in Eq. ( 7) can also be rewritten,

[Tg 0 ](r) ν = S ′ S i dS i p(r i ← r) dS ′ p(r ′ ← r), (9)
which represents the fraction of particles that originate from r making one collision and then ending up on the source area S i . Since the next term of the series is simply the previous one with an additional T operator applied, it becomes obvious, by recurrence, that each term is the fraction of particles that originates from r and finishes on the area S i with an increasing number of collisions on the path they follow. The sum in Eq. ( 7), divided by ν, is thus simply the total fraction of particles leaving r that end up on S i where they are not re-emitted. If one generated N random paths starting at point r, then a numerical approximation of Eq. ( 2) would be the arithmetic average given by

g(r) ν = lim N→∞ 1 N N i=1 N i , (10)
with N i given by N i = 1,i th particle exits through source 0,i th particle exits elsewhere .

If the emission pattern was not uniform on the inlet surface, then Eq. ( 11) needs to be modified. Indeed, the path weight would differ from 1 when the particle would exit at the source area. This weight would depend on where on the source area the particle exited in the following manner: it is the ratio between the inlet emission flux at the point of exit, noted ν(r Exit ), and the average flux emitted by the inlet area (S source ), 1/S source Source dSν(r). This ratio gives 1 in the case of uniform emission at the inlet.

In the case discussed here, every direction generated in a path history is made with a spatial cosine distribution of angles with respect to the surface normal, and this from start to end. Of course, in a path history a collision with an exit or a source area makes the particle vanish. However, if the starting point is not a wall, the cosine distribution remains the one that needs to be generated. Indeed, if the goal is to get the flux on a imaginary surface within the volume of the system, one needs to define an orientation of the surface and generate the paths as if the surface was a wall. However, during the trajectory, the particle cannot hit this imaginary surface since it is not a system wall. Such cases include the need to compute the flux on the exit area or the returning flux on the source area. One notes that the flux computed on a imaginary surface with a certain normal unit vector is different if the normal is taken in the opposite direction. The difference between the opposite fluxes gives the net flux going through the imaginary area.

To sum up, in reverse path tracing, when generating paths from the desired point, one needs only to count the number of particles that end up on the source area to compute the impinging flux. This is in contrast to direct path tracing, where the particles are generated at the particle source area and one counts the particles that make it to a small discrete area around the desired point. Since g(r)/ν is the fraction of particles that reach the particle source, then by definition w(r) = 1 -g(r)/ν must be the Clausing function defined in Sec. I.

B. Convergence and error

Looking closely at Eq. ( 11), one sees that counting particles that exit through the source area is a binary event, a random walk process. From this point of view, the convergence of the calculation can be estimated from the fluctuations of the random walk process, which is quantified by its standard deviation. Not having the Bernoulli trial probabilities, which are actually what is being computed, one must use the sample variance of Eq. ( 10) given by

σ 2 ≈ N N -1 1 N N i=1 (N i ) 2 -( N) 2 , with N = 1 N N i=1 N i . ( 12
)
Convergence of the computation is evaluated using the definition given in the following, called precision from now on and noted δ:

δ ≡ 2 N σ 2 N × 100. ( 13
)
However, the precision of the computation does not guarantee that the result is accurate. To evaluate accuracy, noted in Eq. ( 14), the result is compared to trusted literature values, g L (r), computed with other calculation methods:

≡ |g(r) -g L (r)| g L (r) × 100. ( 14
)

C. The distribution function

The distribution function f (r,v) is defined as a density of particles located at point r that have a velocity v in a given direction. In a free molecular flow regime with uniform temperature, the solution to the transport equation gives the distribution function as

f (r,v) = 4 v m 2πk B T 3 2 exp - mv 2 2k B T g r -t v |v| S , (15)
where v = √ 8k B T/πm is the average velocity of a particle in a Maxwell-Boltzmann gas, T the temperature, k B the Boltzmann constant, m the particle mass, S the system surfaces, and (rt v |v| S) is simply the nearest intersection between the reverse path line of the particle and the surface located at a distance t [START_REF] Sparrow | Velocity profile and other local quantities in free-molecule tube flow[END_REF]. From this distribution we shall compute the particle density anywhere in a system's volume,

n(r) = f (r,v)dv, ( 16
)
and compute the local mean collision rate Ŵ(r) of a particle,

Ŵ(r) = σ T dv 1 f (r,v 1 ) dv 2 |v 1 -v 2 |f (r,v 2 ) n(r) , (17)
FIG. 3. Accuracy and precision comparison between direct and reverse path tracing of the Clausing function at the entrance point (z = 0) of a pipe of aspect ratio L/R = 10. It is compared to the number of collisions made in the pipe, which is proportional to the computation time. For the direct path tracing, the pipe was separated in 300 equal sectors. Square points represent the direct path simulation, and the circles the reverse path tracing. Empty dots are for accuracy, and the filled ones are for precision.

where σ T = πd 2 is the total collision cross section for hard spheres of diameter d. The local mean free path is defined with the help of this collision rate:

λ(r) ≡ v Ŵ(r) . (18)

III. REVERSE AND DIRECT PATH-TRACING COMPARISON

Reverse path-tracing samples only the ones that contribute to the impinging flux at one point, which is not the case in direct path tracing. Because of this, the reverse path will always converge faster than direct path tracing for each single point that is considered. This paper attempts to illustrate this with the comparison of a reverse and direct path Monte Carlo simulation in a pipe of aspect ratio length over radius L/R = 10.

The simulation was carried out for multiple points along the tube's height (see Fig. 1), identified as the z axis, and the same conclusion was made for each of those points. From this fact, it is sufficient to illustrate the comparison using only one point, in this case the entrance at z = 0. The reference value of the Clausing function at this point is taken from [START_REF] Lobo | Comparison between monte carlo and analytical calculation of Clausing functions of cylindrical and conical tubes[END_REF]a sw(0) = 0.085 062 5631 giving an impinging flux ratio of g(0)/ν = 1 -w(0). For the direct path Monte Carlo the pipe was segmented in 300 sectors of equal width to accumulate the collisions made within the area of each circular band. Such discretization is not needed in reverse path Monte Carlo, only points where the impinging flux is wanted need to be specified. To illustrate the stability and performance of reverse path tracing, accuracy and precision [Eqs. [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] and [START_REF] Malhotra | High Knudsen number physical vapor deposition: Predicting deposition rates and uniformity[END_REF]] are plotted with respect to the number of collisions made within the tube. This represents roughly the number of operations made by the computer to attain a certain degree of accuracy and precision. Now proceeding with the result of this comparison illustrated in Fig. 3. One striking feature is that reverse path is already stable and converged accurately beneath 0.1% after 1 × 10 7 collisions within the tube. On the other hand, the direct path still fluctuates considerably when compared to reverse path tracing and converged to 1.2% only. This is more than an order of magnitude in computing time gain. Moreover, since the direct path uses a set of discrete bands that extend over a small range of the length axis, the impinging flux cannot be evaluated as accurately for a given number of collisions. This is well represented in the inset of Fig. 3 where after 1 × 10 10 collisions the accuracy is off by 0.075% but fluctuations are now four times less than this offset. This is due to the numerical error introduced by the finite width of the collecting sector.

Indeed, let ε be the collecting section length on the z axis. The average number of incident particles on a circular band of area, g(z) 2πRε, that is computed with direct path tracing should be

g(z) 2πRε = lim N→∞ 2πRε N N-1 i=0 g z - ε 2 + i N -1 ε . ( 19
)
With a Taylor expansion to the second order and a few sum identities it is possible to show that the error made on the flux ratio g(z)/ν is of second order with ε:

g(z) ν - g(z) ν = ε 2 24 d 2 dz 2 g(z) ν + O(ε 3 ). ( 20
)
As stated in Eq. ( 20), the numerical error is proportional to the second derivative of the flux along the pipe (d 2 [g(z)/ν]/dz 2 ). Since the flux along the pipe is not exactly linear as reported by Clausing [START_REF] Clausing | The flow of highly rarefied gases through tubes of arbitrary length[END_REF] this derivative is nonzero, therefor producing a small offset. For this reason, while reverse path tracing gives an accurate result for any point along the pipe, the direct path will always average over the collecting patches, giving a small error when the flux is not linear. This error diminishes when the collecting sector width is reduced, but at the price of increasing the convergence time. Since the flux is about linear in a pipe, the numerical error offset is quite small.

IV. CONTROLLED GAS INJECTION IN HIGH-VACUUM SYSTEMS

With the help of reverse path tracing the flux ratio g(z)/ν along the length of a conical section was computed to a precision below δ 1% for each of the 500 points chosen. The flux ratio also gives access to the distribution function given by Eq. [START_REF] Sharipov | Data on internal rarefied gas flows[END_REF]. The starting radius is 0.95 mm, the length is 20.1 mm, and the opening half angle is 10 • . This shape was initially studied by Isnard [START_REF] Isnard | Theoretical and experimental molecular beam angular distribution studies for gas injection in ultra-high vacuum[END_REF] to experimentally show that the transport equation solution is valid for free molecular flow and that the impinging flux arriving on a wafer destined for epitaxial growth could then be controlled by optimizing the shape of a gas injector. However, too high a flux at the entrance illustrated the inadequacy of the free molecular flow hypothesis to describe the exiting flux. To find out what too much flux represents in certain gas injection designs one must compute the distribution function. From this, it is possible to quantify the density anywhere in the conical section and evaluate the mean free path of molecules. An upper limit to the input flux can then be estimated if the injector is to work in free molecular flow.

Figure 4 shows a density map for uniform temperature of the walls and input gas. The cone is illustrated in the inset with the flow direction. The color gradient represents the density ratio within the system and because of the axial symmetry it needs to be plotted only with respect to radius and height. The calculations assume a uniform input flux ν emitted at the entrance from a Maxwell-Boltzmann reservoir located upstream. The input flux is then given by ν = n i v/4, with n i the input density of the gas and v the mean velocity of a particle defined after Eq. [START_REF] Sharipov | Data on internal rarefied gas flows[END_REF]. The density in the cone varies slightly with the radius, being denser in the middle, and very strongly with the height. Indeed, the particle density drops rapidly after a few millimeters of the entrance. Such a drop helps to keep the free molecular flow for most part of the cone. Even at the very beginning, the density drops to 75% of the input flux. With these data it is possible to estimate the mean free path at the entrance and see when the free molecular flow hypothesis breaks down. Let us attempt to quantify a net flux for which the free molecular flow breaks down, the net flux being what is measured in experiments, ν net = Wν, with W = 0.56 being the cone's Clausing factor defined previously. The source area is given by A = π (9.5 × 10 -4 ) 2 m 2 from which the net particle flow can be computed

where x is the number of net standard cubic centimeter (sccm). Using Eq. ( 21), the input density may be computed from a Maxwell-Boltzmann gas,

n i = 4ν net vW , (22)
with v the average velocity of a particle according to Maxwell-Boltzmann distribution and 4.465 × 10 17 the number of particles per second, equivalent to 1 sccm. With the definition of the mean free path and Eq. ( 22), it is now possible to compute the Knudsen number at the entrance,

λ L = v Ŵ(r)L , (23)
considering an input gas of nitrogen N 2 [START_REF] Isnard | Theoretical and experimental molecular beam angular distribution studies for gas injection in ultra-high vacuum[END_REF]. Then the effective intermolecular collision cross section is σ T = 4.3 × 10 -19 m 2 [START_REF] Hirschfelder | Molecular Theory of Gases and Liquids, Structure of Matter Series[END_REF][START_REF] O'hanlon | A User's Guide to Vacuum Technology[END_REF] and v = 475 m/s at 295 K. With Eqs. ( 10), ( 15), ( 17), ( 21), [START_REF] Hirschfelder | Molecular Theory of Gases and Liquids, Structure of Matter Series[END_REF], and (23) and using a characteristic length of L = 2 mm we compute the Knudsen number for various net flows which yields Fig. 5. The plot suggests that for a net flow of about ν net A = 0.4 sccm the free molecular flow does not hold. This is in agreement with the results reported in Ref. [START_REF] Isnard | Theoretical and experimental molecular beam angular distribution studies for gas injection in ultra-high vacuum[END_REF], where at 1.2 sccm the measured pattern strongly deviates from the free molecular flow simulation. Such a flow represents a Knudsen number of about K n = 0.4 according to Fig. 5, which is in the transition regime.

V. CONCLUSION

Solving the transport equation for the free molecular flow in complex systems with a Monte Carlo method requires sampling techniques to accelerate the convergence of the calculations. It was shown that path reversal can be used when it comes to cosine emitting surfaces, and that the transport equation could be solved using this method. Thus reverse path tracing is a sampling method that is easy to implement and that removes the need for surface discretization, improving both precision and accuracy. An example of these gains was illustrated with the computation of the flux ratio at the entrance of a pipe of aspect ratio L/R = 10. Using the precision and accuracy gain of reverse path tracing, a density map was computed inside a conical section. The entrance density was then used to compute a net flow for which the free molecular flow regime would break down thus altering the performance of a controlled gas injector. Other local physical quantities, like pressure gradient, could have been computed using reverse path tracing but were left out since the goal of this paper was only to demonstrate the usefulness of this technique for vacuum system applications.

FIG. 2 .

 2 FIG. 2. Schematic of the iterative process expressed by the equation g = g 0 + Tg 0 + TTg 0 + TTTg 0 +••• Each direction change on a wall represents an integral over all the re-emitting walls.

FIG. 4 .

 4 FIG. 4. Relative particle density distribution within a cone. The density at the entrance is 0.75 times input density and decreases rapidly with the height of the frustum. The red arrows in the inset indicate the direction of the flow in the conical section.

FIG. 5 .

 5 FIG. 5. Knudsen number in function of the net flow in the conical segment.

ACKNOWLEDGMENTS

We express our gratitude for the financial support from the Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG) and the Fond de recherche Québécois sur la nature et les technologies (FRQNT).