A V Salsac
email: a.salsac@utc.fr

J Walter

D Barthès-Biesel

Presqu'île de Giens (Var) Coupling of finite element and boundary integral methods: Motion of ellipsoidal capsules in simple shear flow

Keywords: fluid-structure interaction, finite element method, boundary integral method

We introduce a novel numerical method to model the fluid-structure interaction between a microcapsule and an external flow. An explicit finite element method is used to model the large deformation of the capsule wall, which is treated as a bidimensional hyperelastic membrane. It is coupled with a boundary integral method to solve for the internal and external Stokes flows. The study aims at characterizing the motion and large deformation of an initially-ellipsoidal capsule subjected to a simple shear flow. Both oblate and prolate spheroids are considered. Two regimes are found depending on the value of the capillary number, ratio of the viscous to elastic forces. At low capillary number, the capsule tumbles, behaving mostly like a solid particle. At higher capillary numbers, the capsule has a fluid-like behavior. It oscillates in the shear plane, while its membrane continuously rotates around its deformed shape. During the tumbling-to-swinging transition, the capsule transits through an almost circular profile in the shear plane, for which a long axis can no longer be defined. Qualitatively, oblate and prolate capsules are found to behave similarly, particularly at large capillary numbers, when the influence of the initial state fades out. However, the capillary number at which the transition occurs is significantly lower for oblate spheroids.

Introduction

Bioartificial capsules, which consist of an internal liquid protected by a thin hyperelastic wall, have many practical uses in pharmaceutical, cosmetic and bioengineering applications. They can also be seen as simplified models of red blood cells. When suspended in another flowing liquid, capsules undergo large deformations and strong fluid-structure viscous coupling due to the low-Reynolds-number flows of the internal and suspending liquids. The membrane mechanical behaviour has been extensively studied numerically for initially-spherical capsules in flow and different strategies have been considered to derive numerical solutions. Many studies have used a fluid solver based on the boundary integral method to solve for the Stokes flow equations (e.g. [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF][START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow : Large deformations and the effect of capsule viscosity[END_REF][START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow : effect of the membrane constitutive law and onset of buckling[END_REF][START_REF] Lac | Deformation of a capsule in simple shear flow : effect of membrane prestress[END_REF][START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF]). The velocity field at any position within the fluid domain is given by surface integrals calculated on the geometric boundaries. This method therefore has two main advantages : it reduces the geometric dimension of the problem by one, which largely decreases the total number of nodes, and avoids re-meshing the fluid domain at each time step, which is hugely costly.

The model the most used for the capsule wall is that of a 2D hyperelastic membrane : the wall is considered to be infinitely thin and to have a negligible bending stiffness. Two approaches may be considered to model the capsule membrane mechanics : the equations of the force equilibrium on the capsule wall may either be written locally at each point (strong form) or globally integrated over the capsule surface (weak form). Most capsule studies have used the strong form of the equations. Capsules in simple shear flows have been considered by Pozrikidis [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF], Ramanujan & Pozrikidis [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow : Large deformations and the effect of capsule viscosity[END_REF] and more recently by Li & Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF], who computed the membrane load as a piecewise constant function. Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow : effect of the membrane constitutive law and onset of buckling[END_REF] and Lac & Barthès-Biesel [START_REF] Lac | Deformation of a capsule in simple shear flow : effect of membrane prestress[END_REF] used instead bi-cubic B-splines as interpolation functions in order to compute the loads with high accuracy. These studies have shown that the capsule wall is dominated by in-plane membrane tensions, but compressive tensions can occur locally, creating wrinkles.

An alternative option is to write the balance equations in their weak form and to use a Finite Element (FE) method. The local equilibrium equations are converted into their variational equivalent. Eggleton & Popel [START_REF] Eggleton | Large deformation of red blood cell ghosts in a simple shear flow[END_REF] and Doddi & Bagchi [START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF] have coupled a FE method with an immersed boundary method. Ho-wever their FE method, based on the use of linear (P1) elements, lacks generality in the implementation. The results are not very precise, which is probably a consequence of the method used for the fluids.

Recently, we have proposed a new method coupling the boundary integral method and the finite element method [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a Stokes flow[END_REF]. It has the advantage to use the same discretization for the fluids and capsule wall, which allows a Lagrangian tracking of the membrane position with high accuracy. The numerical stability and accuracy of the method have been demonstrated on a spherical capsule placed in different linear shear flows. The method has been shown to be very stable in the presence of in-plane compression.

Capsules with a non-spherical initial shape have been considered very little [START_REF] Kessler | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Sui | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Bagchi | Dynamics of nonspherical capsules in shear flow[END_REF], even though they provide more realistic models of red blood cells. The extensive wrinkling that appear when such capsules deform had previously rendered most numerical simulations unstable. We presently study initiallyellipsoidal capsules placed in a simple shear flow such that their revolution axis is in the shear plane. We consider both oblate and prolate spheroids to characterize the influence of the capsule shape on the deformation. Particular attention will be paid to the regimes observed and to the transition between them.

Problem statement

We consider a spheroidal microcapsule of semi-axes denoted a along the revolution axis and b along the orthogonal directions. It consists of an infinitely thin membrane of surface shear elastic modulus G s that encloses an incompressible liquid. It is suspended in a simple shear flow in the (x 1 ,x 2 ) plane :

v ∞ = γx 2 e 1 , (1)
where γ is the shear rate. The inner and outer fluids are supposed to be Newtonian and to have the same viscosity µ and density ρ. Two main parameters influence the capsule deformation. The capillary number

Ca = µγℓ G s (2)
compares the viscous stresses in the flow to the elastic stresses in the membrane. The length scale ℓ is defined as the radius of the sphere having the same volume as the ellipsoid (ℓ

= 3 √ ab 2 )
. The other parameter is the aspect ratio a/b, ratio of the ellipsoid semi-axis lengthes. The spheroid is oblate when a/b < 1 and prolate when a/b > 1.

In the following sections, all studies are conducted for two initial values of the aspect ratio : a/b = 0.5 (oblate spheroid) and a/b = 2 (prolate spheroid). These two aspect ratios are well-suited to comparing the influence of the capsule initial shape. The two spheroids have, by definition, the same internal volume and therefore the same length scale, l. With initial surface areas differing by less than 2%, the two capsules essentially have the same initial value of the surface-area-to-volume ratio.

Membrane mechanics

Following [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF], we treat the membrane as a purely bidimensional sheet of hyperelastic material. In particular, we neglect strain variations across the thickness and therefore the bending stiffness of the material. We introduce the displacement U(X ,t) = x(X ,t) -X, where X and x are the positions of a given material point in the reference and deformed state. It is related to the velocity of the membrane through the kinematic condition :

∂ ∂t U(X ,t) = v(x,t) . (3)
The numerical procedure consists in following the position of the capsule membrane after the start of flow. At each time step, the membrane deformation is thus known : it is quantified by the principal dilation ratios λ 1 and λ 2 in its plane. Two deformation invariants are generally used :

I 1 = λ 2 1 + λ 2 2 -2, I 2 = λ 2 1 λ 2 2 -1 = J 2 s -1. (4)
The Jacobian J s = λ 1 λ 2 represents the ratio of the deformed to the undeformed surface areas. Elastic stresses in an infinitely thin membrane are replaced by the elastic tensions corresponding to forces per unit arc length measured in the plane of the membrane. When the membrane is a twodimensional isotropic material, the Cauchy tension tensor τ can be related to a strain energy function per unit area of undeformed membrane w s (I 1 , I 2 ). Presently we will consider two hyperelastic laws. The neo-Hookean law describes the behaviour of an infinitely thin sheet of a three-dimensional isotropic material

w NH s = G s 2 I 1 -1 + 1 I 2 + 1 . (5)
The area dilation modulus is then K s = 3G s . Skalak et al. [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF] derived a constitutive law for twodimensional materials with independent surface shear and area dilation moduli

w Sk s = G s 4 I 2 1 + 2I 1 -2I 2 +CI 2 2 , C > -1/2 . ( 6
)
The area dilation modulus is then K s = G s (1 + 2C) and the Poisson ratio ν s = C/(1 + C). The Sk law was initially designed to model the area-incompressible membrane of biological cells, such as red blood cells, corresponding to C ≫ 1.

The finite element method is then used to solve the equilibrium of the membrane

∇ s • τ + q = 0. ( 7
)
where q is the unknown viscous load exerted by the fluids on the membrane and ∇ s • the surface divergence operator in the deformed configuration. The local equilibrium is written in a weak form using the virtual work principle. Let V be the Sobolev space H 1 . For any virtual displacement field û ∈ V , the internal and external virtual work balance requires

S û • q dS - S ε( û) : τ dS = 0, (8)
where ε( û) = 1 2 ∇ s û + ∇ s ûT is the virtual strain tensor. The equation is solved for the load q. Note that most finite element procedures look instead for the displacement field, treating the forces as known. This is not possible when dealing with a capsule : the lack of kinematic boundary conditions implies that the displacement solution is not unique.

Internal and external flows

Once the load q is known, the velocity field can be calculated solving for the internal and external flows with the boundary integral method. The flows are indeed governed by the Stokes equations, as the Reynolds number Re = ρl 2 γ/µ is infinitely small. To determine the new position of the capsule at the following time step, one needs to solve for the velocity of the points of the capsule. The latter can be expressed as an integral equation over the deformed capsule surface S :

∀x 0 ∈ S, v(x 0 ) = v ∞ (x 0 ) + 1 8πµ S J(x 0 , x) • q(x) dS , (9)
where v ∞ is the undisturbed flow velocity, J is the Green's single layer kernel, r = x 0x and r = r . The new position of the membrane points is obtained solving (3).

Numerical method

For a given deformed state of the capsule, we first solve the solid problem [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a Stokes flow[END_REF]. V is discretized as a finite element space using an unstructured mesh, with quadratic (P 2 ) curved triangular elements. The discretized problem leads to the following matrix system :

[M]{q N } = {R} , (10)
where {q N } corresponds to the degrees of freedom of the discretized load, [M] has the structure of a mass matrix and {R} corresponds to the right-hand side of ( 8) and depends non-linearly on the displacement U. The tensions τ are computed directly using the strain energy function w s . Solving the solid problem consists in assembling {R} and [M] on the deformed state. Surface integration is performed using Hammer points on the elements. Eq. ( 10) is then solved using the sparse solver Pardiso [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. The velocity field is obtained explicitly from the boundary integral equation ( 9), which is discretized on the same mesh as for the solid problem, once again using Hammer points for the integration. Finally, the new position is obtained by convecting the nodes after integrating (3) numerically with a secondorder Runge-Kutta method. More information on the convergence and stability of the method may be found in [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a Stokes flow[END_REF].

Results

Motion modes : tumbling and swinging

Depending on the value of the capillary number Ca, a given ellipsoidal capsule may exhibit two types of motion. At low flow strength, a "solid-like" regime occurs, called tumbling. Oblate and prolate spheroids exhibit a similar behaviour. Figure 1(a) shows the time evolution of a prolate ellipsoidal capsule with initial aspect ratio a/b = 2 (a/ℓ = 1.59) and a Sk membrane (C = 1) for Ca = 0.1. It shows that, during tumbling, the capsule rotates like a quasi-rigid ellipsoid subjected to the flow vorticity, while the internal flow is almost stationary with respect to the membrane. This is illustrated by the fact that point P that was initially on the small axis of the ellipsoid, remains in the vicinity of its initial location.

Since the capsule profile may be quite difficult to characterize, we evaluate the capsule distortion by the deformation of its ellipsoid of inertia [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow : Large deformations and the effect of capsule viscosity[END_REF]. This method is widely used, but gives approximate results when the deformed particle shape is far from ellipsoidal. By symmetry, the material points initially located in the shear plane (x 1 , x 2 ) remain in it and two of the principal axes of the ellipsoid of inertia with semi-axes L 1 and L 2 (L 1 ≥ L 2 ) are also located in the shear plane. Correspondingly, it is convenient to quantify the three-dimensional capsule deformation with the deformation of the intersection of the profile with the shear plane. The Taylor deformation parameter is then defined as

D 12 = L 1 -L 2 L 1 + L 2 . ( 11
)
Note that, contrary to a spherical capsule, the initial value of D 12 is not zero and is given by D 0 12 = |a -b|/(a + b). The time-evolution of the Taylor parameter D 12 is shown in figure 2a. The large oscillations prove that the membrane can undergo large displacements without large deformation, which is due to the initial shape flaccidity (measured by the surface area to volume ratio). At higher flow strength, a "fluid-like" regime occurs, called swinging. The membrane rotates around the deformed shape of the capsule similarly to the tank-treading motion observed for spherical capsules. Because of the initial anisotropy of the reference shape, the membrane points are not equivalent : a stationary steady state is not possible and a periodic motion of the deformed capsule occurs. Like in the tumbling regime, oblate and prolate spheroids behave similarly in the swinging regime. We illustrate the salient motion features for the same ellipsoidal capsule at Ca = 0.9 (fig. 1c). The capsule assumes an elongated shape with a long axis aligned with the maximum flow strain direction, while the membrane continuously rotates around the deformed shape. The time evolution of the Taylor parameter D 12 is also shown in figure 2a. It can be noted that D 12 oscillates around a mean value that is larger than in the tumbling regime. The capsule reaches its maximum deformation, when the material points located originally on the longer axis of the ellipsoid are in the straining direction (see fig. 1c, γt = 23.5). Conversely, D 12 is minimum when the points originally on the smaller axis are in the straining direction, or equivalently, when the points originally on the longer axis are aligned with the flow compression direction (γt = 30.5).

Tumbling-swinging transition

In order to study the tumbling-to-swinging transition, we consider the same capsule (a/b = 2) at a capillary number Ca = 0.3. One characteristic difference between the tumbling and swinging regimes is the time-evolution of the capsule long axis. In the tumbling regime, the long axis rotates over time, while it oscillates around a mean value in the swinging regime. This is how the transition has been previously defined [START_REF] Kessler | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Sui | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Bagchi | Dynamics of nonspherical capsules in shear flow[END_REF]. However, we believe that such a criterion is not quite appropriate to determine the critical capillary number Ca ⋆ at which the transition occurs. Indeed, figure 1(b) shows that during each halfcycle, the capsule transits through an almost circular profile in the shear plane (here at γt = 26). At this time, the two principal axes of the ellipsoid of inertia in the shear plane have roughly the same length (L 1 ≈ L 2 ), thus there is no clearly identified "long" axis. The fact that the capsule profile becomes quasi circular in the shear plane shows that the transition corresponds to a Taylor parameter getting close to zero at its minimum (figure 2(b). In conclusion, we propose to define the critical capillary number Ca ⋆ as an interval using the criterion min D 12 < 0.05 ,

min D 12 being the minimum value of the Taylor parameter over one period. This criterion corresponds to a relative difference of 10% between the lengths L 1 and L 2 of the axes, which is the minimum difference that can be quantified with precision. Figure 1 also shows that the transition is associated with an increased wrinkling of the capsule wall. The most extensive wrinkling is seen to occur as D 12 approaches 0, i.e. as the long axis of the undeformed ellipsoid becomes shorter. However, wrinkling is transient during the cycle and the wrinkles disappear when the capsule long axis is in the direction of the viscous stretch (γt = 20.5 in figure 1b).

Effect of capillary number

We now conduct a systematic study of the motion of a spheroidal capsule in a simple shear flow as a function of the capillary number. Two material laws are used to describe the membrane, the neo-Hookean law ( 5) and Skalak's law ( 6) with C = 1, which have the same behaviour at small deformation. Two initial aspect ratios are studied : a/b = 0.5 and a/b = 2. This makes it possible to compare the behaviour of oblate and prolate spheroids, but a systematic study of the influence of the aspect ratio is outside the scope of this article.

The minimum value min D 12 is shown as a function of Ca and membrane law in figure 3a. A larger value of Ca is required to reach the same value of D 12 with the Sk law than with the NH law. If the NH and Sk (C = 1) laws behave similarly at small deformation levels, they are known to diverge at larger deformation, the Sk law exhibiting a strain-hardening behaviour and the NH law a strain-softening one.

Figure 3a also proves that min D 12 indeed goes through a global minimum in all cases. The criterion min D 12 < 0.05 provides values of Ca ⋆ confined within a small interval, because of the sharp variations of min D 12 around the global minimum. To determine the intervals for the critical capillary number Ca ⋆ , the capillary number was increased systematically by steps of 0.01 for oblate spheroids and 0.05 for prolate spheroids. The values of capillary number for which min D 12 < 0.05 are provided in table 1. Note that, for a given aspect ratio, the Ca ⋆ intervals are almost equal for the two material laws considered. This is due to the fact that the transition takes place at moderate deformation levels, for which the two laws behave similarly. However, oblate and prolate capsules have values of Ca ⋆ that differ by a factor ∼ 10. As the 2 capsules studied (a/b = 0.5, 2) share the same initial value for the Taylor parameter in the shear plane (D 0 12 = 0.33), any difference in their behaviour is due to the difference in their initial geometry in the orthogonal direction. The prolate ellipsoidal capsule therefore requires a higher energy to reach the shape for which D 12 ∼ 0 than an oblate capsule, because of its short characteristic length in the orthogonal direction.

It is shown in figure 1 that widespread in-plane compression can occur. In the absence of a physical bending stiffness in the numerical model, such compressive tensions cause numerical wrinkles. In order to study in-plane compression, figure 3b shows the maximum value over one period of τ min (t), the minimum principal tension, denoted max τ min :

max τ min = max t (τ min (t)) = max t min x,i=1,2 (τ i (x,t)) , (13)
where τ i are the principal tensions. In all the cases studied, τ min is negative through most of the period indicating that compression always occurs somewhere for spheroidal capsules. The positive values of τ min occur when the long axis of the original ellipsoid is in the straining direction.

Figure 3b shows that, for large values of Ca, even the maximum value of τ min is negative. It means are provided for an initially-spherical capsule for reference [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow : effect of the membrane constitutive law and onset of buckling[END_REF].

a/b = 0.5 a/b = 1 a/b = 2 NH Ca ⋆ ∈ [0.
that negative tensions occur even when the capsule reaches its maximum elongation. The reason is that, at large values of Ca, negative tensions and wrinkles appear at the tips of the elongated capsule. Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow : effect of the membrane constitutive law and onset of buckling[END_REF] observed this phenomenon for an initially-spherical capsule (a/b = 1) and defined Ca H as the capillary number above which negative tensions appear at steady state. In the case of spheroidal capsules, we define Ca H as the critical capillary number above which max τ min < 0. The values of Ca H found for the different cases studied are provided in table 1.

Except for a prolate NH spheroid, the largest amount of wrinkling occurs in the swinging regime, for values of Ca slightly above Ca ⋆ . During tank-treading motion, strong wrinkling tends to occur when the long axis of the initial ellipsoid has to be compressed to become the short axis of the deformed capsule. However, as the capillary number is increased and the capsule becomes more elongated, the isotropic component of the tensions (related to the Poisson ratio of the membrane) increases and compensates the negative tensions : wrinkling becomes less important. The maximum amount of wrinkling therefore occurs during transition and for capillary numbers slightly above it.

In the case of a prolate NH spheroid, the wrinkling does not subside as Ca increases. This is a consequence of the proximity of the two critical capillary numbers Ca ⋆ ∈ [0.20, 0.25] and Ca H = 0.35 in this particular case. Indeed, if we consider a material point originally on the long axis of the ellipsoid, at Ca ≈ 0.35 the point is on the short axis of the deformed capsule and strong wrinkling occurs, since Ca is only slightly above Ca ⋆ . A quarter of a period later, the material point is in the straining direction, but buckling and wrinkling occur at the tips, as Ca is around Ca H . These two phenomena then lead to a constant wrinkling of the membrane, that even seems to amplify over time, but this is probably a numerical artifact due to the lack of a proper bending stiffness in the model of the capsule wall.

Conclusion

We have modelled the behaviour and large deformation of an ellipsoidal capsule in a simple shear flow using the novel method of Walter et al. [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a Stokes flow[END_REF], that couples boundary integrals for the flows to finite membrane elements. The study has shown that the coupling method is well-suited to the simulation of non-spherical capsules and that it remains numerically stable in the presence of in-plane compression, even though the mechanical wall model does not account for bending stiffness. This is a convenient feature when the bending effects remain localised and weak. Still, to model exactly the physical behaviour of the capsule wall, a proper shell model remains to be implemented.

The coupling method allowed us to study the behaviour of oblate and prolate spheroids, with aspect ratios a/b = 0.5 and a/b = 2, and to recognise two regimes : a quasi-solid regime ('tumbling') at low capillary numbers, where the long axis of the capsule rotates in the shear plane, and a quasi-fluid regime ('swinging') at high capillary numbers, where this axis oscillates around a mean inclination and membrane rotation (tank treading) occurs. These two regimes are separated by a transition region, during which the capsule transits through a phase where the two axes of the capsule in the shear plane are approximately of the same length (D 12 ≈ 0). For the transition to occur, the critical stage is for the points initially located on the small axis to have enough energy to pass the long axis of the deformed capsule. The transition can thus be understood as the crossing of an energy barrier as proposed by Skotheim & Secomb [START_REF] Skotheim | Red Blood Cells and Other Nonspherical Capsules in Shear Flow : Oscillatory Dynamics and the Tank-Treading-to-Tumbling Transition[END_REF]. However, in their semi-analytical theory, they supposed that the shape of the capsules remained constant and postulated a periodic variation of the strain energy of the membrane. The present study shows that the capsule shape changes over time and that the energy variation is thus more complex than predicted.

We have also shown that oblate and prolate spheroids behave qualitatively similarly in most respects. Still the capsule initial shape influences the transition between the tumbling and swinging regimes, which occurs at a much lower value of Ca for oblate spheroids. The tumbling regime is therefore limited to a very small range of capillary numbers for oblate spheroids, as the value of the critical capillary number at transition is low.

We do not believe that the transition can be considered as a separate regime, distinct from tumbling and swinging. It rather corresponds to the parameter range where the two regimes behave so closely that they cannot be accurately distinguished from one another. Along with the present study, all the existing numerical studies of capsules in simple shear flow indicate that the tumbling-to-swinging transition is associated with a phase when D 12 ≈ 0 [START_REF] Kessler | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Sui | Swinging and tumbling of elastic capsules in shear flow[END_REF][START_REF] Bagchi | Dynamics of nonspherical capsules in shear flow[END_REF]. A similar behaviour has been found experimentally for lipid vesicles [START_REF] Kantsler | Transition to Tumbling and Two Regimes of Tumbling Motion of a Vesicle in Shear Flow[END_REF][START_REF] Deschamps | Phase Diagram of Single Vesicle Dynamical States in Shear Flow[END_REF]. These findings are, however, at odds with the experiments conducted by Abkarian et al. [START_REF] Abkarian | Swinging of Red Blood Cells under Shear Flow[END_REF] on red blood cells in a simple shear flow. They observe that, during the tumbling-to-swinging transition, the red blood cell maintains an almost constant shape and that the transition occurs through an intermittent regime during which the cell alternately swings and tumbles. However, our study and the other numerical studies cited above fail to find such an intermittent regime and always observe that the tumbling-to-swinging transition is associated with large variations of the capsule shape.

FIGURE 1 -

 1 FIGURE 1 -Evolution of the capsule shape in the shear plane over one half period. The initial shape is a prolate spheroid, a/b = 2, and the membrane follows the Sk law with C = 1. The grey scale corresponds to the normal component of the load, q • n. The dot shows the position of material point P, originally on the short axis. The value of the non-dimensional time step γt is given below each shape.

FIGURE 2 -

 2 FIGURE 2 -Taylor parameter D 12 as a function of time. (a) solid line : Sk law C = 1, a/b = 2, Ca = 0.1 (tumbling regime), dotted line : Sk law C = 1, a/b = 2, Ca = 0.9 (swinging regime), horizontal line : D 0 12 = 0.33 (b) Sk law C = 1, a/b = 2, Ca = 0.9 (transition).

FIGURE 3 -

 3 FIGURE 3 -(a) Minimum values of the Taylor parameter D 12 of the ellipsoid of inertia. The horizontal line corresponds to the limit min D 12 = 0.05 used to define the transition. (b) Maximum value during a period of the minimum principal tension. Results are shown for an oblate spheroid (a/b = 0.5, filled symbols) and a prolate spheroid (a/b = 2, open symbols). : NH law ; △ : Sk law with C = 1.

TABLE 1 -

 1 Values of the critical capillary numbers Ca ⋆ and Ca H for the cases studied. The values of Ca H

	02, 0.04] Ca H = 0.70	⋆ Ca H = 0.63	Ca ⋆ ∈ [0.20, 0.25] Ca H = 0.35
	Sk, C = 1 Ca ⋆ ∈ [0.02, 0.05] Ca H = 2.5	⋆ Ca H = 2.4	Ca ⋆ ∈ [0.25, 0.35] Ca H = 1.4


J Yvonnet
email: julien.yvonnet@univ-paris-est.fr

A Mitrushchenkov
email: alexander.mitrushchenkov@univ-paris-est.fr

G Chambaud

Q.-C He

Presqu'île de Giens (Var) A multiscale method for modeling size and surface effects in crystalline nanowires based on finite element and quantum mechanics

Keywords: Nanomechanics, Surface effects, Multiscale methods

A multiscale procedure for modeling crystalline nanostructures such as nanowires is proposed. The size effects exhibited by nano objects are captured by taking into account a surface energy, following the classical Gurtin Murdoch surface elasticity theory. An appropriate variational form and a finite element approach are provided to model and solve relevant problems numerically. The methodology is completed with a computational procedure based on ab initio calculations to extract elastic coefficients of general anisotropic surfaces.

Introduction

Nanosized objects have at least one dimension in the nanometer range, roughly between 1 to 100 nm. They are intermediate in size between atomistic/molecular and microscopic/continuum structures. By assembling nanostructural elements, it is possible to manufacture nanocomposites or nanodevices with desired properties and functions. Various types of nanostructural elements like nanowires, nanotubes, nanorods, nanorings, nanobelts or nanocombs have been synthesized and can serve as building blocks for nanosystems and nanostructures.

Nanowires and in general nanostructures are characterized by non zero surface energy that make their effective properties size-dependent. Due to a different local environment, atoms near a free surface or interface have different equilibrium positions than do atoms in the bulk.

Surface free energy is neglected in traditional continuum mechanics because it is associated with only a few layers of atoms. For objects with dimensions larger than tens of nanometers, the ratio between volume and surface is extremely small. However, for nano-sizes particles, wires and films, this ratio becomes significant, and so does the effect of surface free energy. The importance of surface stress has led to important theoretical modeling advances that has started with the works of Gurtin and Murdoch [START_REF] Gurtin | Continuum theory of elastic-material surfaces[END_REF].

This work provides a complete procedure for modeling size-dependent mechanical effects in nanowires. Firstly, a continuum variational framework enriched with surface energy is developed, and a three-dimensional FEM methodology is presented to solve the problem numerically. Secondly, a new technique based on ab initio computations is detailed to determine the surface parameters related to the continuum model.

Continuum modelling with surface energy

We consider an open domain Ω ∈ R 3 with a bounding surface Γ. The latter Γ is composed of two disjoint complementary parts Γ u and Γ F where the displacements and forces are prescribed, respectively. To introduce the surface energy and related size effects, the Gurtin-Murdoch model [START_REF] Gurtin | Continuum theory of elastic-material surfaces[END_REF] of surface elasticity is adopted. The equilibrium equations of the problem are then given by

∇ • σ + b = 0 in Ω, (1)
∇ s • σ s + σn = 0 on Γ (2)
where σ is the Cauchy stress tensor, b is body force and n is the outward unit normal vector to Γ. Equation (1) refers to the bulk equilibrium equation, while Eq. ( 2) refers to the surface equilibrium (see e.g. [START_REF] Chen | FDerivation of the generalized Young-Laplace equation of curved interfaces in nanoscaled solids[END_REF]). Superscript (.) s refers to surface quantities or operators. In Eq. (2), ∇ s • T denotes the surface divergence of a differentiable second-order tensor T and is defined by :

∇ s • T = ∇T : P. (3)
In the above equation, the symbol " :" denotes double contraction of indices and P(x) = 1n(x) ⊗ n(x) is an orthogonal projection operator describing the projection on the plane tangent to Γ at x ∈ Γ. The surface strain and stress tensors are defined by :

σ s = PσP, ε s = PεP, (4)
where ε = 1 2 (∇u + ∇u T ) is the linearized strain tensor and u is a displacement field. It is worth noting that the operator P does not involve any basis change : the reader has to keep in mind that T s denotes the projection of T on the surface Γ. Thus T s = PTP does not express the components relative to the basis of the tangent plane to the surface. This projection tensor nicely allows one to mix bulk and surface quantities in one equation. The surface stress σ s is related to the surface strain by

σ s = ∂γ s 0 (ε s ) ∂ε s + ∂γ s (ε s ) ∂ε s . ( 5
)
Above, γ s (ε s 0 ) and γ s (ε s ) are specified by

γ s 0 (ε s ) = τ s : ε s , γ s (ε s ) = 1 2 ε s : C s : ε s . (6)
Physically, γ s 0 (ε s ) is an energy density due to surface residual stress τ s , γ s (ε s ) is a surface strain density function with C s denoting a surface elastic tensor. We do not assume an isotropic surface, i.e. τ s ̸ = τ s P. The equations (1-2) are completed with boundary conditions prescribed on the surface Γ u and Γ F as

{ σn = F on Γ F u(x) = ū(x) on Γ u . ( 7
)
Finally, it is assumed that the surface adheres to the bulk :

[[u]] = 0 on Γ. (8)
To apply this framework to finite element analysis, we provide in the following the weak form of the equations. It can be derived straightforwardly by taking the classical weak form and taking into account a term related to the virtual internal work of the surface [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF]. Another way to obtain it is to express the potential energy of the system and minimize the energy with respect to the displacement field. Thus, the potential energy of the system is given by :

E = W b (ε) +W s (ε s ) -W ext . ( 9
)
In the above expression, W b (ε) is the bulk elastic strain energy given by

W b (ε) = ∫ Ω 1 2 ε : C bulk : εdΩ. (10)
In (9), W s (ε s ) is the surface elastic energy provided by

W s (ε s ) = ∫ Γ [γ s 0 (ε s ) + γ s (ε s )] dΓ. ( 11
)
Finally, the work of external forces is expressed by : The equilibrium is found by writing the stationarity condition :

W ext = ∫ Γ F u • FdΓ + ∫ Ω u • bdΩ. ( 12
)
D δu E = 0 (13)
where D v f is the Gâteaux derivative of the functional f in the direction v. It leads to the weak form suitable for finite element discretization : Find u ∈ H 1 (Ω) and u = ū on Γ u such as

∫ Ω ε(u) : C bulk : ε(δu)dΩ + ∫ Γ ε s (u) : C s : ε s (δu)dΓ = ∫ Γ F F • δudΓ + ∫ Ω b • δudΓ - ∫ Γ τ s : δε s (δu)dΓ (14)
for all δu ∈ H 1 0 (Ω), H 1 (Ω) and H 1 0 (Ω) being the usual Sobolev spaces.

Discretization with finite elements

In the following the Voigt's notation is adopted. The vectorial form of the stress tensor is given by [σ] = {σ 11 , σ 22 , σ 33 , σ 23 , σ 13 , σ 12 } = {σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 } and the strain tensor counterpart is given by

[ε] = {ε 11 , ε 22 , ε 33 , 2ε 23 , 2ε 13 , 2ε 12 } = {ε 1 , ε 2 , ε 3 , ε 4 , ε 5 , ε 6 }.
The indices 1, 2 and 3 correspond to the x, y and z directions, respectively (see figure 1). In this work we consider the nanowire which is a monocrystal with wurtzite structure, as for example AlN, which grows with hexagonal cross-sections (see figure 1) and whose bulk elastic tensor can be expressed by 5 independent constants C 11 ,C 33 ,C 44 ,C 12 and C 13 .

To solve the problem by FEM, it is necessary to discretize the weak form (14). While the bulk terms lead to classical matrix and vector forms, terms related to surface energy require a specific treatment. The domain is first discretized into volume elements, as shown in figure 1.

To discretize the surface terms, surface elements are required. In [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF], we proposed a framework avoiding the explicit mesh of surfaces when considering surface energy by using an XFEM/level-set approach. In this work however, we adopt the classical FEM framework for the sake of simplicity. Let (e 1 , e 2 ) a local orthonormal basis related to the surface such that e 1 = t and e 2 = z, with z and t = z × n unit vectors along the main direction of the nanowire and tangent to the surface Γ (see figure 1).

The surface residual stress can be expressed by :

τ s = τ s 1 e 1 ⊗ e 1 + τ s 3 e 2 ⊗ e 2 = τ s 1 t ⊗ t + τ s 3 z ⊗ z. ( 15
)
Due to the symmetry of the hexagonal nanowire, the six facets are identical and correspond to the ( 10 10 ) surfaces. Then the surface stress can be related to surface strain through four independent constants C s 11 , C s 13 , C s 33 and C s 55 and two residual stress components τ s 1 and τ s 3 :

  σ s 1 σ s 3 σ s 5   =   C s 11 C s 13 0 C s 13 C s 33 0 0 0 C s 55     ε s 1 ε s 3 2ε s 5   +   τ s 1 τ s 3 0   . ( 16
)
The elastic tensor is expressed in the local basis as :

C s = C s 11 t ⊗ t ⊗ t ⊗ t +C s 33 z ⊗ z ⊗ z ⊗ z +C s 13 (t ⊗ t ⊗ z ⊗ z + z ⊗ z ⊗ t ⊗ t) +C s 55 (t ⊗ z ⊗ t ⊗ z + z ⊗ t ⊗ z ⊗ t). (17)
On substituting the FEM discretization into the weak form (14), and using the arbitrariness of nodal variations δu, the following discrete system of linear equations is obtained :

(K + K s ) q = F -F s (18)
with q being the vector of unknown displacements,

K = ∫ Ω B T C bulk BdΩ ( 19
)
being the bulk rigidity matrix and

K s = ∫ Γ B T M T p C s M p BdΓ (20)
being the surface rigidity matrix. In (20) C s is the matrix form of the tensor (17) such as

[σ s ] = C s [ε s ], with [σ s ] = { σ s 11 , σ s 22 , σ s 33 , σ s 23 , σ s 13 , σ s 12 } and [ε s ] = { ε s 11 , ε s 22 , ε s 33 , 2ε s 23 , 2ε s 13 , 2ε s 12 }
. Note that [σ s ] is different from the left-hand term in (16) which express the surface stress in the local basis (tangent to the surface) while [σ s ] is expressed in the cartesian basis. As C s is symmetric (see Eq. ( 17)), the matrix K s is also symmetric. The matrix M p relates the surface strains to the bulk strains through

[ε s ] = M p [ε]. Its expression in 3D is specified by M p =      
)         . ( 21
)
We finally obtain the following approximation for K s :

K s ≃ ∑ e∈S Γ B T (x e Ω )M T p (x e Γ )C s (x e Γ )M p (x e Γ )B(x e Ω ) |Γ e | ( 22
)
where |Γ e | denotes the area of the triangular element Γ e , x e Ω is an integration point in the element and x e

Γ is an integration point on the external surface of the element (see more details in [START_REF] Yvonnet | Finite element model of ionic nanowires with sizedependent mechanical properties determined by ab initio calculations[END_REF]). Remark that integration points being located inside triangular facets, the unit normal vector n and tangent projector P are well defined. In this work we exclude the surface energies of both ends surfaces where displacement and external forces are prescribed. The generalized force vectors are obtained by

F = ∫ Γ F N T FdΓ + ∫ Ω N T bdΩ (23)
F s = ∫ Γ B T M T p τ s dΓ, ( 24
)
where N is a classical FEM shape functions matrix. In the following, body forces are neglected.

Extracting coefficients using Quantum Mechanics calculations

In the present work, bulk and surface elastic parameters are extracted through atomistic ab initio calculations. In the present work we used the CRYSTAL code [START_REF] Dovesi | CRYSTAL06 User's Manual[END_REF] for all ab initio calculations. B3LYP KS-DFT method with extra large space integration grid, was used. Basis set was 86-21G* for Al [START_REF] Montanari | Influence of the exchange-correlation functional in all-electron calculations of the vibrational frequencies of corundum (alpha-Al2O3)[END_REF], and 6-31Gd1G for N [START_REF] Gatti | rystal-field effects on the topological properties of the electron-density in molecular-crystals -the case of urea[END_REF]. We have chosen the hybrid B3LYP functional because it is known to provide quite accurate energetic properties. As we are interested in elastic behavior, which is evaluated through energy derivatives, this choice seems appropriate.

Surface coefficients

To compute surface elastic parameters, a surface (slab) model shown in figure 2 is used, consisting of n layers of atoms in the direction normal to the surface. Periodic conditions are applied along both other directions.

The possible strains for the slab system are ε s 1 , ε s 3 and ε s 5 which are prescribed along the z-( 3) and t-( 1) directions. The upper and lower surfaces normal to n are free to relax. The total energy of the system is the sum of two components, the surface energy and the bulk energy. The surface energy has to be isolated from the total energy. We assume the following model :

E slab (w) = wE s + (1 -w)E slab (w → 0) , w = 2 n ( 25
)
where w is the relative weight of the surface for a n-layer slab, as the model depicted in figure 2 contains two external atomic layers. Defining the elastic constants as derivatives of the energy with respect to strains, we obtain

C slab i j (w) = wC s i j + (1 -w)C slab i j (w → 0) (26)
where C s i j are the surface elastic properties and C slab i j (w → 0) are the limit values which can be obtained from bulk term by fully relaxing ε 2 for fixed ε 1 and ε 3 and optimizing the energy of the system.

Eq. ( 26) can be re-arranged as

C i j (w) = w(C s i j -C slab i j (w → 0)) +C slab i j (w → 0). ( 27
)
The procedure consists into computing values of C i j for different values of w (by increasing the number of layers n) and fitting the obtained curve with a linear function to identify the slope r = C s i j -C slab i j (w → 0) and thus C s i j . For the residual stress we have :

τ i (w) = wτ s i . ( 28
)
In that case τ s i are directly obtained from the slope of the linear fit (see an illustration in [START_REF] Mitrushchenkov | Towards an elastic model of wurtzite AlN model[END_REF]). We define per-area elastic constants as where S = ac is the unit cell area. The total surface energy E slab is determined using the DFT method described above for several values of strains ε 1 , ε 2 and ε 5 . As for bulk computations, a polynomial fitting was used to extract the derivatives. The resulting surface parameters were provided for the (10 10) AlN wurtzite surfaces in [START_REF] Mitrushchenkov | Towards an elastic model of wurtzite AlN model[END_REF][START_REF] Yvonnet | Finite element model of ionic nanowires with sizedependent mechanical properties determined by ab initio calculations[END_REF] (parameters A).

C slab i j = 1 S ∂ 2 E slab ∂ε i ∂ε j (29)
4 Numerical examples

Validation of the model : nanowire in traction

The aim of this first example is to validate the mechanical continuum model by comparing it with ab initio calculations. Different atomistic models fully solved with an ab initio method are constructed with increasing diameters. We classify the nanowires according to the number of atomic layers n in the nanowires (see figure 1 b). The biggest calculated nanowire corresponds to n = 7 and contains 588 atoms in the unit cell. As both periodic conditions along z-axis and hexagonal symmetries are taken into account, a reduced model containing 98 atoms is employed.

We compute two quantities for both models : the Young modulus Ēzz and the axial strain of the nanowire in absence of loading due to surface residual stress τ s and (c) the bulk radial residual strain ε0 rr (see more results in [START_REF] Yvonnet | Finite element model of ionic nanowires with sizedependent mechanical properties determined by ab initio calculations[END_REF]).

Results are presented in figures 3 (a)-(b). A mesh containing 8905 nodes was used to plot the continuous curves.

We can observe from figure 3 that the agreement between the continuum and ab initio models is good in regards to Young's modulus and axial relaxation strain, when employing the coefficients obtained from the slab procedure without any modification. The use of optimized coefficients [START_REF] Yvonnet | Finite element model of ionic nanowires with sizedependent mechanical properties determined by ab initio calculations[END_REF] gives rise to a slight improvement.

From the above results, we conclude that the model has a good accuracy regarding the effective properties of the nanowire. In what follows, we examine the local fields. A continuous displacement field is constructed from the ab initio discrete atomic displacements.

We can notice from figure 4 that in the bulk region both displacement fields are in very good agreement, especially when the set of optimized surface coefficients B is used. However, near the surface, the local fields are poorly reproduced. This can be partly explained by the fact that both Al and N atoms have very different kinematics [START_REF] Yvonnet | Finite element model of ionic nanowires with sizedependent mechanical properties determined by ab initio calculations[END_REF]. On the surface layer, Al atoms possess a relatively large inward radial displacement as compared to the Al atoms in the next interior layer and with opposite sign to the N atoms which have an outward displacement on the surface. This highly heterogeneous field cannot be captured accurately by the proposed continuum model.

Bending of a nanowire

In this example, the constructed continuum model is employed to study size effects in bending analysis of a nanowire. This problem is not tractable via ab initio calculations, as no periodicity can be considered. This test is of practical importance, as this type of loading can be prescribed for example using an atomic force microscope (AFM). The nanowire is clamped at its both ends and subjected to a force distribution F = σ y e y ⊗ e y n, σ y = 1e -8 MPa on a width l = L/20. As in previous example, the length of the nanowire is taken as L = 4d.

To evaluate the influence of surface effects, we compute the transverse stiffness as a function of the nanowire diameter. The results are presented in figure 6. The size effects can be clearly observed.

FIGURE 1 -

 1 FIGURE 1 -(a) Local basis and FEM model of nanowire with hexagonal cross-section ; (b) ab initio model of nanowire with n layers.

FIGURE 2 -

 2 FIGURE 2 -Slab models with number of layers n.

FIGURE 3 -

 3 FIGURE 3 -(a) Effective Young modulus of the nanowire : comparison between continuum FEM model and ab initio model (prescribed stress) ; (b) Nanowire Bulk radial strain in relaxation : comparison between continuum FEM model and ab initio model

FIGURE 4 -

 4 FIGURE 4 -Radial displacement field along nanowire radial direction (d = 3.11 nm). Comparison between continuum model and ab initio interpolated atomic displacements.

FIGURE 5 -

 5 FIGURE 5 -Deformed configuration of the nanowire in free relaxation (magnified 10 times) for (a) d = 2 nm ; (b) d = 5 nm ; (c) d = 100 nm. Red and black meshes correspond to undeformed and deformed configuration, respectively.

FFIGURE 6 -

 6 FIGURE 6 -(a) ε yy strain field and deformed configuration of the nanowire in bending ; (b) Transverse stiffness of the nanowire versus diameter.

  23 P 13 + P 12 P 33 ) (P 22 P 13 + P 12 P 23 ) (P 13 P 23 + P 12 P 33 ) (P 2 13 + P 11 P 33 ) (P 12 P 13 + P 11 P 23 ) (P 13 P 22 + P 12 P 23 ) (P 11 P 23 + P 13 P 12 ) (P 11 P 22 + P 2 12

		P 2 11 P 2 12 P 2 13	P 2 12 P 2 22 P 2 23	P 2 13 P 2 23 P 2 33
	 	2P 12 P 13 2P 22 P 23 2P 23 P 33 2P 11 P 13 2P 12 P 23 2P 13 P 33
		2P 11 P 12 2P 12 P 22 2P 13 P 23
	P 12 P 13	P 11 P 13		P 11 P 12
	P 22 P 23	P 12 P 23		P 12 P 22
	P 23 P 33	P 13 P 33		P 13 P 23
	(P 2 23 + P 22 P 33 )	(P


V A Yastrebov

G Cailletaud
email: georges.cailletaud@mines-paristech.fr

F Feyel

9-13 Mai 2011, Presqu'île de Giens (Var) Enrichment of the contact geometry within the Finite Element Method

Keywords: Contact, Node-to-Segment contact discretization, Finite Element Method, multiscale approach, enrichment of shape functions, smoothing techniques

A new approach for sub-mesh enrichment of the contact geometry has been proposed in the framework of the Finite Element Method and the Node-to-Segment contact discretization. The method is very general and multipurpose : it allows, for example, to extend the limits of the contact modeling of thin-walled structures or to account for a change of the contact geometry due to loading, for example, to simulate a shallow wear.

Introduction

Contact occurs at the interface between bodies. In the Finite Element Method this interface is discretized by nodes and segments. In the classical master-slave approach and the Node-to-Segment (NTS) discretization, this interface is represented by so-called "contact elements", which consist of one slave node and one or several master segments. In general, the contact geometry of such an element is fully described by the interpolation function of the master segment or by a smoothing function aimed to represent a C 1 geometry over several master segments. However the contact geometry can be enriched independently of the finite element mesh and the related interpolation functions.

In the second half of the 90s several approaches based on enriching of the element interpolation functions have been proposed for different problems [START_REF] Heyliger | Stress intensity factors by enriched mixed finite elements[END_REF], [START_REF] Babuška | A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution[END_REF]. The entire class of enriched or extended finite element methods got the name XFEM (extended finite element method) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] or GFEM (generalized finite element method) or PUM (partition of unity method) [START_REF] Melenk | The partition of unity finite element method : Basic theory and applications[END_REF]. This method is used for modeling of dislocations, solidification, two-fluid flows, cracks and cohesive cracks.

We propose to enrich the geometry of the master surface in a rather similar manner (Fig. 1), as the solution field is enriched within the XFEM, so that the master surface ρ within the contact element is described by the following equation :

ρ e = ρ + h e n ⇔ ρ e ( ξ ∼ ) = ρ( ξ ∼ ) + h e ( ξ ∼ , Θ ∼ )n( ξ ∼ ), (1)
where ρ and ρ e are vectors describing the original and enriched master surfaces respectively, n is a unit normal vector to the original master surface and h e ( ξ ∼ , Θ ∼ ) is an enriching function which depends on the convective coordinate of the master segment(s) ξ ∼ and, in general, may depend on the local strain-stress state and its history. So this approach can be used, for example, to account for geometry change due to wear, deformation of asperities, dislocation escape, relaxation, etc. This dependence is taken into account by means of an array of variables θ ∼ , which can also include time. Remark that the single and double line underlined quantities are vectors and second-order tensors respectively, quantities underlined by a wave represent an ordered array of scalar quantities (e.g., two convective coordinates ξ ∼ = [ξ 1 , ξ 2 ] T , array of parameters θ ∼ , etc.), quantities underlined by a double wave represent matrices of scalars (e.g., a surface metric matrix), the quantities underlined by a line and a wave designate ordered arrays of vectors (e.g., a set of basis vectors), see [START_REF] Yastrebov | Computational contact mechanics : geometry, detection and numerical techniques[END_REF] and [START_REF] Yastrebov | Tensors with tensor components : a new algebra and its applications to geometry and mechanics[END_REF].

There are several motivations to enrich the contact geometry : 1) a multiscale modeling of contact with thin-walled structures with a complex geometry of the surface as well as modeling of contact with a soft bodies with hard coatings ; 2) implicit modeling of anisotropic friction ; 3) account of a complex change of the surface geometry due to loading, etc. The main advantage of the method is an enrichment of the geometry on the sub-mesh scale without increasing the computational costs. The main drawback is a challenging analytical formulation of the geometrical quantities entering in formulae for the residual vector and the tangent matrix of contact elements.

First of all, a general weak form for contact problems will be presented. Next some remarks on the derivation of main equations will be given. In conclusion, some perspectives and applications will be proposed.

F. 1 -Enriched geometry ρ e of the master surface ρ.

Weak formulation of the contact problem

The Finite Element Method is based on the integral (weak) form of the equilibrium equation. A rigorous construction of a weak form for the contact problems results in a variational inequality [START_REF] Kikuchi | Contact Problems in Elasticity : a Study of Variational Inequalities and Finite Element Methods[END_REF]. However this formulation is hard to apply for the case of nonlinear material, large deformation and large sliding contact. That is why for the general case it is advantageous to use the approach of variational equalities [START_REF] Wriggers | Computational Contact Mechanics[END_REF], assuming known in advance the active contact zone Γ c :

Ω σ = •• δε = dΩ + Γ c δW c dΓ c = Γ f σ 0 • δu dΓ + Ω f v • δu dΩ, V = δu ∈ H 1 (Ω) δu = 0 on Γ u , (2)
where σ = is the Cauchy stress tensor, ε = is the Green-Lagrangian strain tensor, Ω is the interior of contacting bodies in the current configuration, σ 0 is a prescribed traction on the surface Γ f and f v is a vector field of volume force density. The displacement vector field u belongs to the Hilbert space H 1 (Ω) and fulfills the Dirichlet boundary conditions u = u 0 on the boundary Γ u . The form of the term δW c , which has to be integrated over the contact area Γ c of any of contacting surfaces, depends upon the choice of the resolution method (penalty, Lagrange multiplier, augmented Lagrangian methods and others) and in general can be written as a linear function of four variations :

δW c = Aδg n + B ∼ T δξ ∼ + Cδλ n + D ∼ T δ λ ∼ t , (3)
where g n is the local normal gap between contacting surfaces (g n > 0 denotes that there is no contact, g n = 0 means that bodies are in contact and g n < 0 designates a local penetration), ξ ∼ is the convective parameter describing the projection point of the slave node on the master segment, λ n is the Lagrange multiplier representing the contact pressure and λ ∼ t is a set of Lagrange multipliers representing the components of the tangential contact stress vector in the local basis 2) needed in an implicit resolution scheme, requires the variation of the virtual work term δW c :

∆δW c = ∂δW c ∂g n ∆g n + A∆δg n + ∂δW c ∂ ξ ∼ T ∆ ξ ∼ + B ∼ T ∆δ ξ ∼ + ∂δW c ∂λ n ∆λ n + ∂δW c ∂ λ ∼ t T ∆ λ ∼ t (4)
The closed form expressions for variations δg n , δξ

∼ , ∆δg n , ∆δ ξ ∼
for the continuous geometry independent on discretization were first obtained in [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems[END_REF]. The terms ∂δW c ∂• depend on the resolution method and remain unchangeable in case of enrichment of the contact geometry. However the variations of the geometrical quantities have to be recalculated.

Main equations for the enriched contact element

To derive a consistent framework one has to impose two conditions on the enriching function h e ( ξ ∼ , Θ ∼ ) :

- ) .

h e ( ξ ∼ , Θ ∼ ) ∈ C 2 ( ξ ∼ ) C 1 ( θ ∼ ), i.e. h e ( ξ ∼ , Θ ∼ ) is a C 2 -
If the master geometry is enriched locally (within each master segments), in order to avoid gaps in the surface these conditions have to be complemented by a third one : the value of the enriching function has to be zero at the boundary of the contact element h e (∂Ω e ) = 0. The latter condition is very restrictive and renders a mesh dependent method. This difficulty can be overcome if instead of local enriching, one enriches the smoothing surface constructed over several master segments, see Fig. 2.

F. 2 -Enriched geometry of the smoothed master surface.

If the enriching function depends only on the surface parameter (i.e. it does not change in time and does not depend on the deformation state) or if it changes very slowly, then all the variations of the kinematic quantities remain the same as in [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems[END_REF] if one replaces all quantities related to the master surface by enriched ones

+ n ∂h e ∂ ξ ∼ + h e ∂n ∂ ξ ∼ , A ≈ e = ∂ρ e ∂ ξ ∼ • ∂ρ e ∂ ξ ∼ T , H ≈ e = n e • ∂ 2 ρ e ∂ ξ ∼ 2 ∂ 2 ρ e ∂ ξ ∼ 2 = ∂ 2 ρ ∂ ξ ∼ 2 + ∂n ∂ ξ ∼ ∂h e ∂ ξ ∼ T + ∂h e ∂ ξ ∼ ∂n ∂ ξ ∼ T + h e ∂ 2 n ∂ ξ ∼ 2 + n ∂ 2 h e ∂ ξ ∼ 2
The explicit forms of these equations can be found in [START_REF] Yastrebov | Computational contact mechanics : geometry, detection and numerical techniques[END_REF]. However to incorporate these modified quantities in a Finite Element framework one has to derive the equations for the basic variations of the enriched master vector and its derivatives

δρ e , δ ∂ρ e ∂ ξ ∼ , δ ∂ 2 ρ e ∂ ξ ∼ 2 .
If the first expression is easily derivable, however to obtain the two latter, a significant efforts have to be undertaken [START_REF] Yastrebov | Computational contact mechanics : geometry, detection and numerical techniques[END_REF]. However if one confies himself to the case of infinitely small gaps g n = 0, as, for example, in [START_REF] Konyukhov | Covariant description for frictional contact problems[END_REF], the formulation becomes much simpler, since in this case there is no need to derive δ

∂ 2 ρ e ∂ ξ ∼
δ ∂ρ e ∂ ξ ∼ → n • δ ∂ρ e ∂ ξ ∼ , ∂ρ ∂ξ ∼ • δ ∂ρ e ∂ ξ ∼ T ; δ ∂ 2 ρ e ∂ ξ ∼ 2 → n • δ ∂ 2 ρ e ∂ ξ ∼ 2 , ∂ρ ∂ξ ∼ • δ ∂ 2 ρ e ∂ ξ ∼ 2 T
, the latter quantities are significantly easier to derive than the original ones.

In general case, in order to avoid these complicated derivations and their programming, one can incorporate in the employed finite element software, an automatic procedure which derives the needed quantities and constructs the residual vector and the tangent matrix for contact elements. This approach has been proposed in [START_REF] Wriggers | Development of 2d smooth polynomial frictional contact element based on a symbolic approach[END_REF] and [START_REF] Korelc | Symbolic approach in computational mechanics[END_REF], where the AceGen software, based on the package of symbolic mathematical calculations Mathematica, has been used to construct the needed quantities for frictional contact in case of smoothed master surface.

It is worth mentioning that the detection procedure has to be also changed. It is often based on the detection of the normal projection of a slave node on the master surface. Obviously in case of enriched geometry the projection point has to be determined on the enriched surface (Fig 3).

As already mentioned, the enriching function h e has to be rather smooth (C 2 ) and its value |h e | must remain smaller than the minimal local curvature radius. Moreover if one uses the normal projection, it is necessary to keep in mind that Newton's method, often used for this purpose, allows to determine only one projection point (closest to the starting point), that is why the function h e has to be chosen properly in order to avoid multiple solutions of min F(r s , ξ ∼ ) (Fig. 4 a), where

F(rs, ξ ∼ ) = 1 2 (r s -ρ( ξ ∼ ) 2
is a distance function. Contrary to the normal projection, the shadow projection proposed in [START_REF] Yastrebov | Computational contact mechanics : geometry, detection and numerical techniques[END_REF] is unique if there is no "shadows" from the master surface on its own, so one has to pay attention to avoid shadows due to an enrichment of the master surface (Fig. 4 b).

F. 3 -The procedure of the normal projection in case of an enriched master geometry (r s is a slave node position).

As already mentioned, in order to preserve the continuity of the master surface in case of local enrichment, enriching functions have to be zero at edges of each master segment h e (∂Ω e ) = 0 (Fig. 4 c). It has to be also mentioned that there is a possibility of intersection of enriched geometries of adjacent master segments, which has to be avoided as well (Fig. 4 d).

Perspectives and applications

The enrichment of the contact geometry by an arbitrary function permits :

1. to take into account a complicated geometry within one contact element ; 2. to account for a change of the local geometry due to loading conditions. As mentioned, if the enrichment is chosen to be localized within NTS contact elements, the choice of the enrichment function is limited : its value must be zero at the edges of the master segments. It implies a strong connection between the finite element mesh and the enrichment. A possible application of this approach is the modeling of periodic structures using a regular mesh, Fig. 5. Enrichment of thin-walled or beam structures by a constant enriching function seems to be meaningfull, since the predominant deformation of such structures does not affect the geometry of the surface (Fig. 5,a-b). A possible applications is a modeling of contact with grid structures, micro contact with fiber, etc. The proposed approach is also valid for the case, when the deformation of the surface geometry is small in comparison to the deformation of the bulk material (Fig. 5,c-d), for instance, hard coating on a soft substrate. An anisotropic friction can be simulated implicitly by a special enrichment of the master surface (Fig. 6). If one couples the enrichment function with a global smoothing procedure (Fig. 6, c-d), the mentioned shortcoming of the enrichment approach vanishes. Since the master surface is globally smooth, it is not anymore required that the enriching function is zero along the edges of the master segments. However it becomes a real challenge to obtain the needed variations of the geometrical quantities for the resulting enriched surface. On the other hand this coupling makes possible to simulate a shallow wear and to enrich the master geometry independently on the mesh.

The proposed method presents a new general and multipurpose framework in computational contact mechanics and requires further developments and applications.

;

  A, C are scalar functions of geometrical quantities and the local stress state and B ∼ , D ∼ are components of vector functions in the local basis. The linearization of the contact term (3) from Eq. (

smooth function by convective coordinate ξ ∼ and C 1 -)| < min i 1 /

 11 smooth by θ ∼ ; -in order to avoid self-intersection of the enriched surface, we require that the value of enriching function remains smaller than the minimal local curvature radius of the surface |h e ( ξ ∼ κ i ( ξ ∼

2 ,

 2 • e : ρ e = ρ + h e n, n e = δρ e = δρ e + ∂ρ e

  equations for the variations of the geometrical quantities as dot products with the normal vector to the enriched surface n e and its local basis ∂ρ e ∂ ξ ∼ , which are in turn can be expressed by the original normal vector n and the local basis ∂ρ ∂ξ ∼ :

F. 4 -

 4 Incorrect choices of enriching function h e : a -multiple normal projection within one segment ; b -presence of self-shadows for shadow projection ; c -discontinuity of enriched master surface ; d -self intersection of enriched master surface.

F. 5 -

 5 Enriched geometry of the master finite element mesh : a-b -periodic thin-walled structure, c-d -mesh of a screw with 4 turns, represented by enrichment of 16 segments with a screw function ρ e of the master surface ρ. F. 6a-b -An example of contact geometry enrichment for modeling of anisotropic friction.