E Sauvage

A Gagnoud

Y Fautrelle

P Brun

J Lacombe

Thermoconvective instabilities of molten glass heated by direct induction in a cold crucible

The aim of this study is to perform numerical simulation of a vitrification process developed by the French atomic commission. This process uses direct induction technology in a cold crucible to melt specific glass at high temperature (~1250°C).

The numerical simulation of this process needs a coupled approach of the different phenomena: induction, thermal and hydrodynamic. Indeed, those three phenomena are strongly coupled because of the temperature dependence of the glass properties. For example, the hotter the molten glass, the higher the electrical conductivity.

Two softwares are used to achieve the purpose of the study. Thermal and hydrodynamic aspects are modeled with Fluent® whereas induction is modeled with Flux®. A home-made iterative coupling is performed between the two softwares which leads to a full 3D simulation of the molten glass heated by induction.

Buoyancy and Marangoni convection are taken into account. Thermoconvective instabilities appear within the glass bath when the total Joule power injected reaches a specific threshold. Qualitative comparison of the aspect of free surface is performed. The direct measurement of the temperature in the glass bath gives a quantitative comparison between experimental and numerical results. Reasonably good agreement is found in both cases.

Introduction

A direct induction process using the cold crucible technology has been developed in the French Atomic Commissariat for 35 years. In this process, walls of the crucible are cooled by internal water circulation which protect them from the corrosion and allow higher temperature of the molten glass. Moreover, a mechanical stirrer and air bubbling system promote a good homogeneity of the melt. In results, the lifetime of the process as the incorporation rate of the nuclear waste is increased.

This study is a continuation of previous works modeling this process (Jacoutot [START_REF] Jacoutot | Strategy of coupling to model physical phenomena within molten glass bath heated by direct induction[END_REF][START_REF] Jacoutot | Modélisation numérique de phénomènes couplés dans des bains de verre brassés mécaniquement et élaborés en creuset froid inductif[END_REF]). The aim is to model the molten glass heated by direct induction. Thus, hydrodynamic, thermal and induction phenomena are taken into account. As the physical properties of the glass are functions of the temperature, these three phenomena are strongly coupled. Consequently, a coupling between two softwares is achieved. Fluent® is used to solve hydrodynamic and thermal equations whereas Flux® computes the Maxwell equations. All computations are three-dimensional. The coupling is based on files data transfers. Each software interpolates the field on the computing nodes of the other software.

Results are compared with those of Jacoutot [START_REF] Jacoutot | Modélisation numérique de phénomènes couplés dans des bains de verre brassés mécaniquement et élaborés en creuset froid inductif[END_REF] obtained with simulation in which induction heating is solved under axisymmetric assumption. This former coupling will be called 'pseudo-3D coupling'.

Modelisation

Geometry

The molten glass bath is a cylinder of 185mm high and 250mm radius. In this study the mechanical stirrer is not modeled, only four thermocouples within a sheath made of alumina (electric insulated) are taken into account.

Electromagnetic model

The commercial software Flux® is used to solve induction equations. The crucible as well the inductor are not modeled, they are approximated by a current sheet around the glass. In this coil, an alternative current (I eff = 1500 A and 280 kHz) is imposed. Owing to the high value of the frequency, the quasi-steady approximation is made. Different formulations for induction equations are available. The A-V formulation gives best result in a material with high gradient of electrical conductivity.

0 1 0 0 A A V t A V t V P V § • § • w u u ¨w © ¹ © ¹ § • § • w ¨w © ¹ © ¹ )& )& )& )&
The coil is not meshed so a reduced scalar potential formulation is used. The resolution is achieved with iterative methods such as a conjugate gradient one. The domain is discretised with approximately 200 000 first-order elements.

Thermal-Hydrodynamic model

Fluent® software solves the Navier-Stokes and thermal equations. The flow is assumed to be laminar due to the high viscosity of the glass.

0 0 0 u u u u p T T 0 g t U W U w § • ¨w © ¹ 0 P P t c T u c T T Q t U O w § • ¨w © ¹ & & & & h (4)
The source terms in the right hand side of equation ( 4) is the Joule power density dissipated in the glass. This term is calculated by Flux® as a function of induced current:

2 th j Q and j T V & & & ( 5
)
The convection is driven by two phenomena. First, buoyancy forces are modeled in the Boussinesq approximation (last term of right hand side of equation ( 2)). Second, thermocapillary convection at the surface is taken into account. This effect is due to the variation of the surface tension with the temperature. This dependence is well described with the law:

0 s s T T V V J 0 T where s T J V w w .
For almost all liquid, J is constant and positive and for the glass its value is 4 1 10 . . N m K 1 [START_REF] Parmelee | The surface tensions of molten glass[END_REF]. The boundary condition at the free surface relates the viscous strain with the thermal strain in the direction of the surface, for example, in the radial direction:

surf surf u T z r P J w w w w (6)
Thermal boundary conditions are handled via a global exchange coefficient on the crucible wall and bottom. At the free surface a mixed condition convection-radiation is considered with an emissivity of 0.9. The glass is supposed to be a Newtonian fluid but all the physical properties are complex functions of the temperature. These laws are confidential but Table 1 summarizes order of magnitude for two temperatures. The confinement glass of this study is opaque so there is no need of internal radiation model. In fact, the internal radiation is naturally included in the variation of thermal conductivity with temperature. Tab. 1: Order of magnitude of physical properties of the glass

Coupling Strategy

Mesh to mesh interpolation

The main issue is to couple a finite element and a volume finite based software. Mesh refinement requirements are different for the induction and hydrodynamic phenomena. Consequently, using a unique mesh is not possible and interpolations between the two meshes are done to minimize losses of precision. These methods have to be robust and fast in order to not penalize computation time. Flux® software natively includes these methods but specific functions have been created to allow Fluent® to do these interpolations using the calculated temperature gradient in each cell.

Coupling strategy

The Joule power distribution is calculated every 20s second of flow computation. This time was chosen in regard of the diffusion time which is 10 3 s. As Jacoutot shown, this does not change the results. The time step for Fluent® is 1s. Initialization is performed with the results obtained with the pseudo-3D coupling due to Jacoutot [START_REF] Jacoutot | Modélisation numérique de phénomènes couplés dans des bains de verre brassés mécaniquement et élaborés en creuset froid inductif[END_REF]. In this coupling, the electromagnetic computation was axisymmetric.

The convergence criterion for the iterative coupling is based on the square error between the old and new map of temperature. It is calculated by the software Flux® during importation of the new field of temperature, the convergence threshold is set at 10 -8 .

Results

Total Joule power injected of 45kW

In this first case, a 45kW Joule power heating is imposed in the glass container. Figure 2 shows iso-contours of temperature and Joule heating in a vertical plane including thermocouples. Cold glass in contact with the cold crucible appears clearly. The maximum Joule power density is located in the hot parts of the bath which have a better electrical conductivity. For such a Joule heating power, glass flow due to natural convection and thermocapillarity is stable and stationary. Without thermocouples, an axisymmetric flow appears. Such a result has already been reported by Jacoutot [START_REF] Jacoutot | Strategy of coupling to model physical phenomena within molten glass bath heated by direct induction[END_REF]. 3D simulations only exhibit a small perturbation of the Joule heating field near the thermocouples.

Total Joule power injected of 55kW

For this total Joule power rate, the flow obtained with the pseudo-3D coupling exhibits thermoconvectives instabilities [START_REF] Jacoutot | Modélisation numérique de phénomènes couplés dans des bains de verre brassés mécaniquement et élaborés en creuset froid inductif[END_REF]. These instabilities are convection cells of rising hot glass similar to Bénard-Marangoni cells [START_REF] Sauvage | Comparison between numerical and experimental results on themoconvective instabilities of a high-Prandtl-number liquid[END_REF] (cf. figure 4).

When the 3D-3D coupling is performed, the Joule power concentrates in the hotter glass zone increasing temperature differences. Thus, the instability is reinforced, and the convections cell growths and their number decreased. Moreover, this instability was initially stationary (non time dependent). The figure 5 shows the evolution of the number of convective cells with physical time. The initialization exhibits 22 cells. During the first 1000s, this number decreased to 11 because cells growths, then cells number decrease a second time because of loss of periodicity, and cells appears and disappears in the glass. Cell number oscillates between 5 and 9. Unfortunately, the construction of the classical non-dimensional numbers as Rayleigh or Marangoni is very difficult because physical properties of the glass are strongly dependent of the temperature.

Characteristic time of evolution of convective cell is about ~15min which is in qualitative agreement with experimental visualization. Number of instabilities in the glass could not be determined experimentally because of the limited view angle of the camera, but size of cells (~10cm) is in agreement (cf. figure 6).

Figure 7 shows temperature measurement with the thermocouples compared with calculated temperature with the pseudo-3D and 3D-3D coupling. Results are in good agreement too.

Conclusion

The 3D coupling enlightens the flow instability, that may occur in the molten glass heated by induction. It has been demonstrated that three-dimensional effect of induction heating could have a significant impact on the flow configuration. If a certain threshold of the total injected power is reached, then 3D thermoconvective cells appear. This instability becomes unsteady and presents chaotic behavior. The wave numbers and also characteristic times are in good agreement with the experimental data, obtained in CEA Marcoule. However, due to the difficulty inherent to the choice of reference values in term of thermophysical properties, it was not been possible, yet, to derive relevant non-dimensional number.

Authors (1) SIMAP -EPM laboratory (2)

Fig. 1 :

 1 Fig. 1: Process overview (left) and Fluent® mesh (right). The four black elements are sheath for thermocouples made of alumina.

Fig. 2 :

 2 Temperature isocontours (a), Joule power density isocontours (b) and velocity vectors (c) on crossing vertical plane including thermocouples. Scales are respectively [300;1670]K, [0 ;7,76.10 6 ]W.m í3 and [0 ;4,68.10 -4 ] m.s í1

Figure 3

 3 Figure3shows how induced current in the molten glass avoid the sheath of thermocouple made of alumina (electrical insulated). The Joule power then concentrates on both side of the sheath creating two local overheating zones.

Fig. 3 :

 3 Fig. 3: Joule power isocontours and induced current vectors at the free surface near thermocouple TC3.

Fig. 4 :

 4 Fig. 4: Temperature isocontours (left) and Joule power density isocontours (right) at the free surface and in a crossing vertical plane. Scales are respectively [310;1450]K and [0 ;4.10 6 ]W.m í3 .

Fig. 5 :

 5 Fig. 5: Time evolution of the number of convection cells.

Fig. 6 :

 6 Fig. 6: Numerical (left) and experimental (right) temperature isocontours at the free surface of the molten glass.

Fig. 7 :

 7 Fig. 7: Comparison of measured temperature and calculated ones at the location of the thermocouples.

& & & & &(1)

& & & & & & & &(2)T u u

& & & &(3)2

Acknowledgment

We wish to thank the Cedrat company for the loan of their software Flux®. Web: www.cedrat.com Email: cedrat@cedrat.com


A M Korsunsky

A Constantinescu

The influence of indenter bluntness on the apparent contact stiffness of thin coatings

In the present paper, the influence of punch tip sharpness on the interpretation of indentation measurements is considered. Firstly, in order to obtain analytical insight into the nature of the problem, closed form solutions are presented for the indentation of a homogeneous elastic half-space by an axisymmetric indenter of arbitrary shape, including Hertzian, conical, and conical indenter with a rounded tip. Next, a fast and efficient numerical implementation of a semi-analytical approach to the solution of problems about frictionless axisymmetric indentation of an elastic layer on a dissimilar substrate is described. The approach allows rapid determination of the load-displacement curves for an arbitrary punch shape. The aim of the study was to explore the implications of a finite indenter tip radius for the problem of property identification of thin films. The variation of the apparent substrate stiffness with indentation depth was established for several generic indenter shapes, namely for flat, conical and parabolic punches, and for a conical punch with a rounded tip. It is demonstrated that in each of these cases the depth variation can be described by a simple curve belonging to a family of two-parametric functions. On the basis of these findings we assess the efficiency of using different punch shapes for property determination. A procedure is proposed for this purpose which relies on the use of the depth variation profiles established here. Finally, the influence of imperfect punch shape on the accuracy of analysis is considered.

Introduction

Modern surface engineering relies on the use of thin coatings to achieve the required improvements in the stiffness, hardness, friction, scratch and wear resistance over the uncoated material. Although a variety of approaches are being developed, direct measurement of the coating's mechanical characteristics is often difficult. Frequently the key properties are assessed indirectly from indentation experiments.

In essence, the task objective of the indentation experiment is then to 'back out' the desired parameters (hardness, modulus) for the system (coating in particular) by interpreting the data obtained from continuous recording indentation testing (CRIT). However, even when sample response can be assumed purely elastic, this inverse problem (indentation interpretation) remains challenging. As in most cases of inverse problems, the direct problem has to be addressed first.

The main difficulty of modelling the indentation experiment stems from the nature of the contact boundary condition between the indenter and the sample. The contact condition is expressed in the form of an inequality that expresses mutual non-interpenetrability of the indenter and sample. Otherwise the problem may be thought as a mixed boundary value problem (part displacement, part traction boundary conditions) with unknown boundary that may move and is determined by both the applied load and the substrate deformation response. These effects result in problem non-linearity even in the case of a linearly elastic sample. Furthermore, the problem becomes more and more difficult with progressively thinner coatings.

Nevertheless, for some contact problems, complete solutions of the indentation problem have been obtained analytically, as described in general reviews of the classical results on the subject in [START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Gladwell | Contact Problems in the Classical Theory of Elasticity[END_REF]. Literature reviews indicate that closed form elastic contact solutions for coated systems (see for example [START_REF] Li | Elastic field of a thin-film/substrate system under an axisymmetric loading[END_REF][START_REF] Schwarzer | Arbitrary load distribution on a layered half space[END_REF][START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF]), or finite element calculations for elastoplastic systems using the finite element method (see for example the comparison of 2D and 3D computations by Lichinchi et al. [START_REF] Lichinchi | Simulation of Berkovich nanoindentation experiments on thin films using finite element method[END_REF]) show that they are difficult to obtain even when all material properties are known. This ensures that in practice parametric studies for the identification of material properties are outside a practical reach.

In order to overcome these difficulties, in the last decades a series of simple response functions have been proposed that relate directly the measured indentation depth and force to the material parameters. A series of papers in the case of homogenous half-spaces are based on the pioneering results by Tabor [START_REF] Tabor | The Hardness of Metals[END_REF] and more recently by Doerner and Nix [START_REF] Doerner | A method for interpreting the data from depth sensing indentation instruments[END_REF] or Oliver and Pharr [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing instruments[END_REF].

In the paper we propose a simplified response function for the indentation of an elastic coated system. The proposed response function uses relative indentation depth (i.e. indentation depth divided by the coating thickness) as the argument, and aims to compute the indentation force using the elastic properties of the coating and substrate as parameters. In fact, the model response function proposed here represents the apparent elastic modulus. This is achieved by normalising the indentation force by the indentation solution corresponding to a homogenous half-space.

In order to assess the efficiency of representation achieved by the model response function, a comparison is made with detailed numerical simulation results. Numerical solutions were obtained for frictionless indentation of an layered half-space by punches of practically relevant geometries (flat punch, Hertzian indenter, cone, and cone with a tip blunted to a sphere or parabola), representing the more realistic indenter tip shape. The mathematical implementation of the desired solutions is achieved by solving an integral equation using the method proposed by Lebedev and Uflyand [10], and previously utilised and reported by Yu et al. [START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF]. The selection of this method from a series possible solutions [START_REF] Schwarzer | Arbitrary load distribution on a layered half space[END_REF][START_REF] Li | Elastic field of a thin-film/substrate system under an axisymmetric loading[END_REF] is motivated by the ease of programming for a large class of indentor shapes, and acceptable stability of the solution even for extreme values of layer thickness, i.e. the cases of very thin or thick coating.

The experimental results for hardness have been obtained from the data on nanoindentation of thin coatings on tool steels [START_REF] Korsunsky | On the hardness of coated systems[END_REF].

The present paper introduces and validates the use of the proposed response function. The parametric study carried out here

• shows that for all considered indentor shapes the apparent elastic modulus can be represented in a compact and efficient manner by a two-parameter function; • provides estimates of the elastic moduli of the coating from a small number of measurements, and indicates the range of relative indentation depths needed for obtaining reliable assessments of coating properties; • provides estimates of measurement errors arising from the use of blunted (rounded) conical indenters when the analysis is carried out as if the indenter were a perfect cone. The blunting radius may also be estimated from the measurement results.

It follows that the model response function can be used for the identification of elastic moduli, and provides significant simplification of laborious computations required for the underlying direct and inverse problems [START_REF] Tardieu | On the determination of elastic coefficients from indentation experiments[END_REF].

It is worth noting that the present analysis is focused on the consideration of elastic deformation of the coating and substrate in response to indentation by axisymmetric rigid indenters. The limitations of such an approach are apparent, and can be identified up front. Firstly, in real indentation experiments involving almost all classes of materials, the response of the material contains the combination of elastic and plastic behaviour. Although for ideally sharp indenters the onset of plasticity is immediate, in practice some degree of indenter bluntness is always present, and, in fact, requires careful calibration [START_REF] Pharr | Understanding nanoindentation unloading curves[END_REF]. At low indentation loads (and indentation depths) it is conventional to identify the regime where the response is dominated by elasticity. At higher loads and deeper penetrations plasticity becomes the dominant deformation mechanism.

The linear elastic models considered in the present paper have their applicability clearly delimited by the onset of plasticity within the coating or the substrate. It is apparent that, whatever the material considered, plastic effects will come for the fore at certain stage in the indentation experiment; in the notation used in this paper, as the relative indentation depth (RID) becomes greater. It should therefore be understood that model predictions for large RID values (RID ≫ 1) can only be viewed as extrapolations of the linear elastic model, for the case of ideally elastic material. In practice, deviations from these predictions due to plastic flow are expected to manifest themselves.

The value of the proposed approach, nevertheless, lies in the clear identification of the trends and interrelationships between test parameters in the indentation experiment. For example, it is apparent that even in the case of purely elastic response evaluation of the coating-only properties must take into account the substrateindenter interaction, via the intermediate of the coating layer. Furthermore, limits of resolution and sensitivity of the system can be clearly identified.

A simple response function for the apparent contact modulus

The material parameter that governs the indentation response of elastic homogeneous substrates can be referred to as the contact modulus, also known as the plain strain modulus. It is defined as:

E T = E 1 -m 2 : ð1Þ
In the present study we assume that the indenter is rigid and undeformable. In cases when the indenter tip deformation needs to be taken into consideration, a combined contact modulus should be defined in an appropriate way. For the coated system, the contact modulus appearing in Eq. ( 1) can be expressed directly as a function of the applied load P, the displacement of the indenter d and the contact radius a H (d) for the equivalent homogeneous substrate problem:

E 4 = E 1 -m 2 = P 2da H d ð Þ : ð2Þ
General solution techniques using integral transform methods for frictionless elastic contact mechanics problems for homogeneous substrates were described by Sneddon. The functional dependence of the contact radius a H on the indenter displacement d for important indenter geometries, such as the cone and sphere, is found to have the form of a power law. This dependence, a=a H (d), can be used to eliminate the contact radius from the expression for the contact modulus. Consequently, E ⁎ can be defined in terms of the applied load and indenter displacement.

For a coated system the apparent contact modulus is therefore defined as:

E 4 d ð Þ = P 2fda H d ð Þ : ð3Þ
For coated systems, the apparent contact modulus varies with the indentation displacement d. It is important to stress again that the above definition contains the function a H (d) which is derived from the solution for the given indenter shape upon a homogeneous substrate, and not the real contact radius a. The dimensionless parameter f appearing in the above expression is a number calculated from the punch shape alone, in a manner defined in the following section.

The solution for the contact radius as a function of the applied force or indenter displacement is not available in closed form for the case of coated systems. Therefore, a closed form expression for the apparent contact modulus cannot be obtained. This parameter is, however, of great importance for the purposes of experimental indentation analysis using depth sensing indentation, since it provides a measure of the composite response of the substrate-coating system in terms of the parameters directly recorded by the testing system.

The purpose of this paper is to present the numerical results for the apparent contact modulus for elastic coated systems, and to give a simple representation of the apparent contact modulus as a function of the indentation depth (Fig. 1).

We show that for coated systems a good approximation for the entire family of numerical results is given by the following proposed model response function [START_REF] Korsunsky | On the hardness of coated systems[END_REF]:

E 4 d ð Þ = e 2 4 + e 1 4 -e 2 4 1 + 1 β 0 d h η : ð4Þ
Here h is the film thickness, and e 1 ⁎ , e 2 ⁎ , β 0 and η are positive constants to be determined from fitting. It is clear from the above formula that e 2 ⁎ represents the apparent contact modulus obtained in the extreme of very deep indentation. Similarly, e 1 ⁎ is the apparent contact modulus found from very shallow indentation. The significance of all parameters is discussed further in the text.

The elastic indentation problem

In order to analyse a series of indentation problems and check the validity of the proposed model function, we have chosen to compute the solution of the elastic indentation problem using a boundary integral approach proposed by Yu et al. [START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF]. To illustrate the basis of the solution procedure and provide some confidence in the validity of the proposed numerical solution in the sequel we shall sketch the main steps of the solution proposed by Yu et al. [START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF] and describe the main features of our numerical implementation. We may further note that the symbolic form of the equations can be readily obtained using the Tensor2Analysis and Mathematica package defined in Constantinescu and Korsunsky, 2007 [START_REF] Constantinescu | Elasticity with Mathematica[END_REF].

Let us consider an elastic layer of thickness h (designated by the subscript j = 1) lying on a dissimilar semi-infinite elastic substrate (designated by the subscript j = 2). The elastic solution can be expressed in the cylindrical coordinate system as:

2μ j u rj = -u j;r -zψ j;r ð5Þ
2μ j u zj = -k j ψ j -u j;z -zψ j;z ð6Þ
σ zzj = 2 1 -m j ψ j;z -u j;zz -zψ j;zz ð7Þ
τ rzj = 1-2m j ψ j;r -u j;rz -zψ j;rz ð8Þ
where u rj , u zj are the components of the displacement vector and σ zj , σ rzj are the components of the stress tensor. These elastic fields can be expressed in terms of a pair of harmonic Papkovich-Neuber potentials ψ j , φ j , j = 1,2: 1

ψ 1 r; z ð Þ = Z ∞ 0 A 1 cosh λz + A 2 sinh λz ð Þ J 0 λr ð Þ sinh λh dλ ð9Þ
u 1 r; z ð Þ = Z ∞ 0 A 4 cosh λz + A 3 sinh λz ð Þ J 0 λr ð Þ sinh λh dλ ð10Þ
ψ 2 r; z ð Þ = Z ∞ 0 A 5 exp -λ z -h ð Þ ð Þ J 0 λr ð Þdλ ð11Þ
u 2 r; z ð Þ = Z ∞ 0 A 6 exp -λ z -h ð Þ ð Þ J 0 λr ð Þ λ dλ: ð12Þ
The unknown functions appearing in the potential formulation above can be determined from the boundary conditions at the surface z = 0:

u z1 r; 0 ð Þ = d -δ r ð Þ 0 V r V a ð Þ ð 13Þ
σ zz1 r; 0 ð Þ = 0 a V r V ∞ ð Þ ð 14Þ
σ rz1 r; 0 ð Þ = 0 0 V r V ∞ ð Þ : ð15Þ
The conditions at the film-substrate interface, z = h, 0≤ r ≤ ∞, depend on the nature of bonding. For a perfectly bonded film these are

u z1 r; h ð Þ = u z2 r; h ð Þ ð16Þ
u r1 r; h ð Þ = u r2 r; h ð Þ ð17Þ
σ zz1 r; h ð Þ = σ zz2 r; h ð Þ ð18Þ
σ rz1 r; h ð Þ = σ rz2 r; h ð Þ ð19Þ
and for a freely sliding film they are:

u z1 r; h ð Þ = u z2 r; h ð Þ ð20Þ
σ zz1 r; h ð Þ = σ zz2 r; h ð Þ ð21Þ
σ rz1 r; h ð Þ = σ rz2 r; h ð Þ = 0: ð22Þ
1 For a general presentation of displacement potentials see [START_REF] Solomon | Elasticité linèaire[END_REF].

Table 1

Different combinations of Young's moduli and Poisson's ratios used in the computations, and the corresponding plane strain moduli of the substrates and layers. The determination of the elastic fields now requires the determination of the unknown functions A i (λ), i = 1, 6. This is accomplished in two steps.

v E E ⁎ E E ⁎ 0.
Firstly, the boundary conditions at the interface are used together with the condition [START_REF] Korsunsky | The influence of punch blunting on elastic indentation response[END_REF] to express all A i , (i = 1, 6) in terms of A 1 .

Secondly, the sought function A 1 is expressed in terms of another unknown function, H, as follows

A 1 λ ð Þ = 1 -g λ ð Þ ½ sinh λh Z a 0 2μ 1 d π 1 -m 1 ð Þ H t = a ð Þcos λtdt: ð23Þ
As a consequence, the surface contact conditions [START_REF] Pharr | Understanding nanoindentation unloading curves[END_REF][START_REF] Solomon | Elasticité linèaire[END_REF] are now expressed as the following single integral equation

H τ ð Þ - 1 π Z 1 0 M y; τ ð ÞH y ð Þdy = F τ ð Þ; 0 V τ V 1 ð Þ : ð24Þ
The parameter τ represents the normalised radial coordinate within the patch of contact, τ = t / a. The terms in the above equation relate to some specific aspects of the contact problem. Function F(τ) depends solely on the indenter shape and the indentation depth. The kernel M(y, τ) depends on the elastic properties of the layer and substrate, the thickness of the former, and the contact radius a. It is defined by M(y, τ) = K(y + τ) + K(y -τ), where

K u ð Þ = a h Z ∞ 0 g w ð Þcos auw h dw; ð25Þ
where h is the layer thickness, and a is the contact radius. We note for the future that the key computational step employed in the solution concerns the evaluation of integral in Eq. ( 25). It should be reported that, particularly for thin coatings, the kernel of this integral becomes oscillatory, and requires careful treatment. This is particularly important since the integral needs to be evaluated repeatedly, particularly so in the case of incomplete contacts, when the extent of the contact patch needs to be found by iteration.

The approach adopted to the evaluation of this integral involved (i) identifying the "tail" where the behaviour of the kernel is monotonic and in fact allows analytical evaluation; (ii) use of special adaptive quadrature in Mathematica [START_REF] Wolfram | The Mathematica Book[END_REF] to evaluate the oscillatory part of the integral, and (iii) implementation of the resulting module in a compiled form, allowing faster computation. All these measures together resulted in adequately good performance and efficient evaluation of the integral in question (see below).

The function g appearing in Eq. ( 23) and the expression for the kernel (25) is non-singular, and depends on the combination of elastic properties of the layer and substrate, and on the conditions at the interface. The expressions for g(w) in [START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF] for the perfectly bonded layer contain a typographical error. The authors of that paper have kindly provided with a correction upon our request. We present the correct formulae in Appendix A.

In particular, in the case of a layer perfectly bonded to a substrate of identical material it is found that g(w) = 0. Eq. ( 25) for the kernel must then be understood in the generalised function sense, and yields a constant value, which results in the integral equation kernel M(y, τ) being identically equal to zero. The homogeneous substrate solution is therefore given by H(τ) = F(τ).

Cases when incomplete contacts are considered (i.e. the extent of contact is not known a priori from the punch shape) introduce a nonlinearity in the formulation, and force the corresponding solution procedure to be iterative. However, even in this case important computational savings can be achieved if the lack of coupling between the elastic properties and the punch shape is explored.

The punch shape can be defined with respect to its tip, i.e. the point of first contact, which for the axisymmetric case will be assumed to lie on the axis of rotational symmetry. A function of the form δ(r / a), δ(0) = 0 was used by Yu et al. [START_REF] Yu | The effect of substrate on the elastic properties of films determined by the indentation testaxisymmetric Boussinesq problem[END_REF] to define the punch shape. We use a different convention here, writing δ(r), δ(0) = 0. This formulation is preferred since the punch shape is manifestly independent of the contact radius a, which is generally not known a priori. The only modification that ensues from this choice is that the function F(τ) appearing in the right hand side of Eq. ( 24) is now given by

F τ ð Þ = 1 - aτ d Z π = 2 0 δ V aτ sin θ ð Þ dθ: ð26Þ
Finally, the total indentation load is obtained by integrating the pressure distribution under the indenter, and is given by

P = 2adE 1 4 Z 1 0 H τ ð Þdτ; ð27Þ
where

E 1 ⁎ = E 1 / (1 -v 1 2
) is the plane strain modulus of the coating.

The solutions for frictionless axisymmetric indentation of a homogeneous substrate can be obtained analytically in closed form for a Fig. 3. The apparent contact modulus versus relative indentation depth for perfectly bonded and freely sliding film respectively, indented with a cone.

Table 2

Computed apparent contact moduli for the indentation of a perfectly bonded layer with a cone. The condition of smooth contact at the contact perimeter assumes the form F(τ) = 0 (note that this condition does not generally hold for a layered system). The relationship between the indenter displacement d and the contact radius a H can therefore be found explicitly by solving Eq. ( 26) with respect to a. The homogeneous substrate solution is therefore characterised by this functional dependence of contact radius on indenter displacement, a =a H (d), and the following value of the indentation load:

v 1 E 1 ⁎ v 2 E 2 ⁎ e 1 ⁎
P H d ð Þ = 2da H d ð ÞE 1 4 Z 1 0 F τ ð Þdτ: ð28Þ
This formula is the basis of the apparent contact modulus definition given in the introduction. Indeed, for indenter of any chosen shape on a layered material it is now possible to normalise the solution for the coated system with respect to the values for homogeneous material.

The application of this analysis to the interpretation of modelling and experimental results is as follows. Let the indenter shape be known, and a pair of indenter displacement and load values be obtained in some way, by measurement or calculation. Based on Eq. ( 28), the apparent contact modulus can be computed by

E 4 = P 2fda H d ð Þ = P 2da H d ð Þ Z 1 0 F τ ð Þdτ ! -1 : ð29Þ
Note that for the evaluation of this parameter, only the knowledge of the displacement-load pair of values is required, and that of the contact radius-indenter displacement function a H (d) for the homogeneous solution. The parameter f defined implicitly in the above expression is

f = Z 1 0 F τ ð Þdτ = 1 - a H d Z 1 0 τdτ Z π = 2 0 δ V aτ sin θ ð Þ dθ: ð30Þ
For example, for the flat punch δ′ = 0, and hence f = 1. Using the results of the following section for a H (d), it is easy to establish that for the cone f = 1/2, and that for the Herztian indenter f = 2/3.

The numerical implementation of the system of equations described above has been programmed in Mathematica [START_REF] Wolfram | The Mathematica Book[END_REF]. The solution was found through an iterative procedure.

A Gauss-Legendre quadrature rule and the roots of Legendre polynomials as collocation points were used for the representation of the function H, as well as for computing the kernel M. This approach allows Eq. ( 24) to be transformed into a linear algebraic system which can be easily solved for a given value of the contact radius a.

The order of approximation for Eqs. ( 25) and ( 24) was varied to establish convergence, the value of 30 being usually sufficient for both. Augmenting the number of collocation points above 30 does not necessary enhance the accuracy of the solution, for reasons associated with the deteriorating accuracy of Legendre interpolation. As a consequence, the kernel computed with a large number of points may be close to singular.

The correct value of the contact radius is found by imposing the condition of vanishing contact pressure at the edge of the contact region: σ zz1 (a / h) = 0, i.e. H(1) = 0 using our representation of the solution. This is done using an iterative procedure using e.g. the dichotomy method, starting with an initial value corresponding to the contact radius in the case of a homogenous substrate.

In order to speed up the numerical computations we have used one of the memory control options of the Mathematica programme (f[x_] := f[x] = …) to avoid repeated computations of the same expression, for example of the integral equation kernel (25). We also used the automatic compilation procedure (Compile) in order to speed up the evaluation of functional calls. The combined use of these options improved the computing time by a factor of 20.

The accuracy of the solution procedure was tested by comparison with the well known analytical solution for the homogenous substrate.

Results and discussion

In the sequel we shall present a series of numerical results for the indentation of perfectly bonded or freely sliding films on a substrate with frictionless indenters of different shapes (Fig. 2). The chosen values of the elastic properties for the substrate and the moduli are shown in Table 1, which shows that both stiff on compliant and compliant on stiff layer-substrate configurations were considered.

In the following figures we present the evolution of the apparent contact modulus with the relative indentation depth. We compare the computed values of the apparent contact modulus, denoted by solid dots, with the model response function (4), found by fitting these data, and represented by a continuous line. The least-squares fitting algorithm used was the standard Levenberg-Marquardt method implemented within Mathematica (function NonlinearFit).

Fig. 4. The apparent contact modulus versus relative indentation depth for perfectly bonded and freely sliding film respectively, indented with a sphere of radius R = 0.1.

Table 3

Computed apparent contact moduli for the indentation of a freely sliding layer with a cone. Please note that in all figures the same set of coating-substrate parameter combinations is used, as listed in Table 1. Individual legends are therefore omitted to avoid crowding the graphs.

v 1 E 1 ⁎ v 2 E 2 ⁎ e 1 ⁎

Conical indenter

The conical indenter shape presents special significance for indentation problems, since it is the only axially symmetric punch shape which does not possess an inherent length scale. As a consequence the indentation problem for a homogeneous substrate becomes self-similar. From this result it follows immediately that the load for indentation into homogeneous substrate must be proportional to the contact area. We have already mentioned in the preceding discussion the function a H (d) that describes the relationship between the contact radius and the indenter displacement into homogeneous material. This function is fixed for any given indenter shape. It is easy to demonstrate that for the cone this function is linear, and depends solely on the cone half-angle.

Indeed, for the cone of half-angle α, the indenter shape is given by

δ r ð Þ = r cot α;
and hence the right side of Eq. ( 26) is

F τ ð Þ = 1 - π 2 aτ d cot α: ð31Þ
The contact radius is found by requiring F(1) = 0, and hence for the conical indenter the contact radius is given by the linear function of indenter displacement,

a H d ð Þ = 2d π tan α: ð32Þ
It also follows from the application of Eq. ( 30) that f = 1/2 for the cone.

Therefore the apparent contact modulus is defined for the cone as:

E cone 4 d ð Þ = P 2fda H d ð Þ = π 2 P d 2 cot α: ð33Þ
For the homogeneous substrate, the indentation load is proportional to the square of the indenter displacement. The contact modulus, which is still defined by Eq. ( 29), is independent of the indentation depth and load.

For the coated system, the contact radius increases with applied load, and the region of significant deformation gradually extends across the layer and into the substrate. The dependence of the contact modulus and contact radius on the displacement deviates from the simple forms found for the homogeneous case. The apparent contact modulus undergoes a transition from the value characteristic for the coating to that for the substrate (Fig. 3). The transition zone depends on the Poisson ratios of the two materials and the boundary condition between the layer and the substrate. In all cases the transition is well represented by the model response function [START_REF] Schwarzer | Arbitrary load distribution on a layered half space[END_REF].

An inspection of the transition zone indicates that the apparent contact modulus gives a good estimate of the layer plain strain modulus up to relative indentation depth of 0.05 for a compliant layer. For a stiff layer the relative indentation depth of 0.1 already leads to an error exceeding 10%. However, the values estimated using the model response function (see Tables 2 and3) are always within 2%-10% of the real values. The large error in the case of the compliant layer comes from the nature of the problem, since the apparent contact modulus even at the relative indentation depth of 10 still shows a 10%-20% difference from the target value.

Spherical indenter

Spherical indenters are particularly important for the determination of elastic properties of materials (compared e.g. with conical Fig. 5. The apparent contact modulus versus relative indentation depth for perfect bonded and respectively a freely sliding film indented with a sphere of radius R = 1.

Fig. 6. The apparent contact modulus versus relative indentation depth for perfectly bonded and freely sliding film respectively, indented with a sphere of radius r = 10. indenters), since the former do not immediately induce plasticity at the very early stages of the indentation process.

The shape of the spherical indenter given to the first significant order of approximation by:

δ r ð Þ = r 2 2R
and consequently the right hand side in ( 24) is

F τ ð Þ = 1 -a 2 τ 2 = Rd;
where R is the indenter radius. It also follows from the application of Eq. ( 30) that f = 2/3 for the Hertzian indenter. The contact radius is given for a homogenous system by:

a H d ð Þ = Rd ð Þ 1 2 ;
and therefore the apparent contact modulus is defined as follows:

E sphere 4 d ð Þ = P 2fda H d ð Þ = 3P 4 1 Rd 3 1 2 : ð34Þ
In the case of spherical indenters with different radii, the numerical examples presented in Figs. 456show that the apparent contact modulus follows a similar evolution as in the case of the conical indenter, from the value of the layer to the value of the substrate with increasing relative indentation depth. However, the transition zone is extended over a larger range of indentation depths and depends on the sphere radius.

The model response function again describes well the evolution of the apparent contact modulus. The parameter values are given in Tables (5-7).

Blunted conical indenter

The case of blunted conical indenters presents a practical interest as experimental indenters are never perfectly conical. The analysis of these indenters is especially important in cases where the film thickness is of the same order of magnitude as the blunting radius of the indenter.

For a conical indenter blunted to the tip radius of R the shape is given by

δ r ð Þ = r 2 2R if r V b b R r -b = 2 ð Þ if r N b: 8 > > < > > :
According to the above definition, the rounded region extends as far as the blunt spot radius of b. The cone half-angle α is defined by cot α = b / R.

The right hand side of Eq. ( 24) is given by [15]

F τ ð Þ = 1 - a 2 τ 2 bd cot α if aτ V b1 - a 2 τ 2 bd cotα 1 - ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 - b 2 a 2 τ 2 s + b aτ arccos b aτ 0 @ 1 A if aτ N b: 8 > > > > > < > > > > > :
Asymptotic analysis of the above equation shows that for small contact radii function F is identical with that for the Hertzian indenter with tip radius R = cot α / b, while for large contact radii the function F of Eq. (31) emerges, characteristic of the conical indenter.

In practice it is often very difficult to determine precisely the degree of blunting of a nominally sharp conical or pyramidal indenter. Consequently it is common to ensure that the blunting radius does not exceed a certain allowed value, and interpretation is carried out by assuming that the indenter is infinitely sharp. Unfortunately, when this approach is applied to very thin coatings, it may be the source of significant error. Using our solution we can assess the error in estimating the apparent contact modulus that results from assuming that the indenter is infinitely sharp.

Following the logic of measurement interpretation described above, it is natural to define the apparent contact modulus in the same way as for the cone, i.e.

E blunt 4 d ð Þ = E cone 4 d ð Þ = P 2fda Hcone = π 2 P d 2 cot α: ð35Þ
The consequences of adopting this definition are illustrated in Fig. 7. The values of apparent contact modulus interpreted as described are compared with those obtained using the results of indentation using a perfectly sharp cone. It is clear that even a degree of blunting as small as b = 0.1 results in severe deviation from the perfect cone results at shallow loads, while the use of even blunter indenters results in massive overprediction of the apparent contact modulus. Fig. 8. The apparent contact modulus versus relative contact radius for a perfectly bonded and a freely sliding film respectively, indented with a series of flat punches.

Flat punch

The flat punch solution is a particularly simple case to analyse, since the extent of contact is known, and inversion of Eq. ( 24) does not require iteration. Also, since a is independent of indentation depth in this case, Eq. ( 27) shows that the indentation load depends linearly on displacement, and from Eq. ( 29) the apparent contact modulus is constant for a given punch. As a consequence, the analysis of load-displacement traces for any fixed punch radius does not reveal any depth-resolving information about the material properties of the system.

However, if the punch radius a (we imply that it is normalised by layer thickness h, as all other linear dimensions) is allowed to vary, then the full range of the coated system response to be examined. This is illustrated in Fig. 8, where the apparent contact modulus is plotted against the relative contact radius, which in this context replaces the relative indentation depth. Simple qualitative analysis confirms that when a / h ≪ 1, the zone of highest stresses is confined to a region in the vicinity of contact, and lies entirely inside the surface Consequently, the apparent contact modulus can be expected to have the value close to that of the layer. In the other extreme, when a flat cylindrical punch has the radius which greatly exceeds the layer thickness, then the apparent contact modulus is likely to approach that of the substrate.

For very large punch radii we remark that in Fig. 8 some variation is observed in the apparent contact moduli. The source of this fluctuation is numerical instability linked to the large values of a/h. It can be eliminated, to a certain extent, by choosing a larger number of discretisation points. However, the calculation becomes too inefficient and prone to accumulation of errors. Proper asymptotic analysis would be the correct approach to this class of problems.

The detailed nature of the transition between the two extremes was determined using the solution method described in the previous section. The apparent contact modulus in this case is defined as

E flat 4 a ð Þ = P 2da :
Two cases were considered: that of perfect bonding at the layer/ substrate interface, and that of frictionless adhesion, i.e. when only normal tractions are transmitted. The results show that the model response function based on the relative contact radius,

E 4 a ð Þ = e 2 4 + e 1 4 -e 2 4 1 + 1 β 0 a h η ð36Þ
once again affords a good representation of the transition. The identified parameters e 1 ⁎ , e 2 ⁎ lie close to the apparent contact modulus of the layer and substrate (see Table 4).

This last result suggests that the apparent contact modulus could be represented as a function of the relative contact radius for other indenter shapes as well. Indeed, the results for the conical indenter presented in Fig. 9 show a pattern similar to the previous plots (cf. also Tables 11 and12).

For the sake of completeness, we list the results for apparent contact modulus and model response function parameter determination in the tables given in Appendix B.

Conclusions

We have presented the application of a simple model response function to the analysis of indentation of elastic coated systems. We consider the utility of this approach to lie in the simplicity with which the elastic response can be described over several decades of relative indentation depths or loads. The use of the model response function to fit the numerical results reveals the asymptotic behavior of the apparent contact modulus of coated systems.

It may be expected that in the extreme of very shallow indentation, β ≪ 1, the corresponding parameter e 1 ⁎ ought to approach the plane strain modulus of the coating E 1 ⁎ . Our results appear to confirm this expectation. Similarly, one might also expect that in the extreme of very deep indentation, β ≫1, the corresponding parameter e 2 ⁎ ought to approach the plane strain modulus of the substrate E 2 ⁎ . This statement is much more difficult to justify. It also appears that our results do not always follow this pattern. A definitive argument on this subject may only be mounted on the basis of detailed asymptotic analysis.

The elastic models used in the present study were chosen due to their elegance and compactness: they allow the analysis of parameter and indenter shape sensitivity in a way that is far more efficient than, for example, the use of finite element modeling. As a means of encapsulating the results of these modeling efforts in a form that is readily accessible and easy to visualize, we used once again the Model Response Function (MRF) approach. Previously this method has been applied to the consideration of experimental data on the depth-dependence of composite hardness.

The fundamental approach used in previous studies was to consider the expenditure of energy during indentation. This was chosen because energy analysis allows one to obviate the need for detailed complete analysis of deformation fields that in the case of specific material response can be quite laborious. Nevertheless, Table 4 The apparent contact moduli for a freely sliding layer computed from the dependence of load on the relative indentation radius for a flat punch. elastic-plastic finite element modeling was also carried out, and very good agreement with energy-based predictions was obtained, enhancing the level off confidence in the choice of predictive model. In the case of purely elastic material response, the present results demonstrate that a combined analytical-numerical treatment of the problem becomes possible, so that detailed characterization of the deformation and stress fields can be obtained. For that reason energy analysis is not employed here.

v 1 E 1 ⁎ v 2 E 2 ⁎ e 1 ⁎
Another important conclusion that is substantiated and formulated in the present paper is that the same Model Response Function successfully applied to plasticity analysis (hardness-depth variation) can also be used to describe the apparent contact stiffness variation with depth. The physical meaning of the adopted function remains the same: it provides a compact, heuristic description of the relationship between important parameters of the indentation experiment.

Let us introduce the following notations for different expressions involving the elastic moduli of the layer and of the substrate:

μ 1 = E 1 2 1 + m 1 ð Þ μ 2 = E 2 2 1 + m 2 ð Þ η = a 2 μ 1 a 1 μ 2 β = μ 1 μ 2 ð37Þ
a 1 = 1 -m 1 b 1 = 1 -2m 1 k 1 = 3 -4m 1 ð38Þ
a 2 = 1 -m 2 b 2 = 1 -2m 2 k 2 = 3 -4m 2 : ð39Þ
The preceding expressions will be combined to form new symbols:

B 1 = 4a 1 a 2 r1 + b 2 β ð Þ β ð40Þ
B 2 = -1 -3b 2 β + k 2 -2b 2 2 β 2 + b 2 k 2 β 3 ð41Þ
B 3 = 2a 2 1 + 2 k 1 + 2b2k 1 β + k 2 β 2 β ð42Þ
B 4 = k 1 + -b 2 + 2b2k 1 + 4a 1 a 2 ð Þ β + k 1 + 4a 1 a 2 k 2 + 2b 1 b 2 2 β 2 + b 2 k 2 β 3 ð43Þ
B 5 = -1 -1 + 4b 2 ð Þ β + 1 -2b 2 ð Þ k2β 2 + k 2 2 β 3 ð44Þ
B 6 = k 1 + 1 + 2k 1 k2 ð Þ β + k2 + 4a2b2k1 + 16a 1 a2 2 + 2b 1 b 2 2 β 2 + k 2 2 β 3 ð45Þ
C 1 = 4a 2 1 1 + b2β ð Þ ð 46Þ
C 2 = 1 + 3b 2 β -k 2 -2b 2 2 β 2 -b 2 k 2 β 3 ð47Þ
C 3 = k 1 -1 -2k 1 ð Þ b 2 β + k2 + 16a 1 a2 2 + 2b 1 b 2 2 β 2 + b 2 k 2 β3 ð48Þ
C 4 = 4a 2 1 + 4a 1 2a2 + a 1 k 2 ð Þ β + 8a 1 a 2 b 1 β 2 ð49Þ
C 5 = 1 + k 2 + 2b 2 ð Þ β -1 -2b 2 ð Þ k 2 β 2 -k 2 2 β 3 : ð50Þ
Finally, for a perfectly bonded layer the function g is expressed as:

g b = 1 - B 1 + B 2 w + B 3 sinh 2 w e -w + B 4 + B 5 w + B 6 sinh 2 w sinhw C 1 + C 2 w 2 + C 3 sinh 2 w À Á e -w + C 4 + C 5 w 2 + B 6 sinh 2 w À Á sinh w ð51Þ
and for a freely sliding layer as:

g s = 1 - sinh 2 w + η w + sinh w cosh w ð Þ w + sinhw coshw + η sinh 2 w -w 2 À Á : ð52Þ
Appendix B. Additional tables

Fig. 1 .

 1 Fig. 1. A schematic diagram of the coated system during indentation.

Fig. 2 .

 2 Fig. 2. Different shapes of indenters considered in this paper.

Fig. 7 .

 7 Fig. 7. The apparent contact modulus versus relative indentation depths for a freely sliding film indented with a perfect cone, and with blunted cones with blunt spot radii b =0.1, 1, 10.

Table 5

 5 Computed apparent contact moduli for the indentation of a perfectly bonded layer with a sphere of radius R = 1.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	484.	109.	0.115	0.593
	0.3	500	0.3	100	484.	82.	0.117	0.554
	0.4	500	0.2	100	491.	69.	0.121	0.535
	0.2	100	0.4	500	103.	474.	2.234	0.797
	0.3	100	0.3	500	102.	426.	3.024	0.734
	0.4	100	0.2	500	101.	417.	4.074	0.702

Table 6

 6 Computed apparent contact moduli for the indentation of a freely sliding layer with a sphere of radius R = 1.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	482.	71.	0.076	0.565
	0.3	500	0.3	100	478.	71.	0.077	0.566
	0.4	500	0.2	100	482.	71.	0.076	0.565
	0.2	100	0.4	500	101.	426.	6.067	0.722
	0.3	100	0.3	500	101.	424.	6.042	0.722
	0.4	100	0.2	500	101.	426.	6.067	0.722

Table 7

 7 Computed apparent contact moduli for the indentation of a freely sliding layer with a sphere of radius R = 10.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	495.	62.	0.738	0.526
	0.3	500	0.3	100	491.	62.	0.745	0.527
	0.4	500	0.2	100	495.	62.	0.738	0.526
	0.2	100	0.4	500	101.	452.	73.	0.695
	0.3	100	0.3	500	101.	449.	73.	0.695
	0.4	100	0.2	500	101.	452.	73.	0.695

Table 8

 8 Computed apparent contact moduli for the indentation of a perfectly bonded layer with a sphere of radius R = 10.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	494.	99.	1.147	0.549
	0.3	500	0.3	100	492.	74.	1.173	0.526
	0.4	500	0.2	100	497.	63.	1.217	0.515
	0.2	100	0.4	500	102.	486.	24.181	0.772
	0.3	100	0.3	500	101.	434.	32.131	0.717
	0.4	100	0.2	500	101.	424.	43.309	0.688

Table 9

 9 Computed apparent contact moduli for the indentation of a perfectly bonded layer with a sphere of radius R = 10.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	439.	125.	0.0153	0.733
	0.3	500	0.3	100	432.	101.	0.0162	0.700
	0.4	500	0.2	100	434.	91.	0.0172	0.685
	0.2	100	0.4	500	23.	658.	0.617	0.368
	0.3	100	0.3	500	70.	561.	0.854	0.450

Table 10

 10 Computed apparent contact moduli for the indentation of a freely sliding layer with a sphere of radius R = 10.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	394.	96.	0.013	0.818
	0.3	500	0.3	100	391.	96.	0.014	0.820
	0.4	500	0.2	100	394.	96.	0.013	0.818
	0.2	100	0.4	500	90.	514.	1.269	0.550
	0.3	100	0.3	500	90.	510.	1.261	0.549
	0.4	100	0.2	500	90.	514.	1.269	0.550

Table 11

 11 Computed apparent contact moduli for the indentation of a freely sliding layer with a cone.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	489.	92.	0.189	1.376
	0.3	500	0.3	100	485.	92.	0.190	1.377
	0.4	500	0.2	100	489.	92.	0.189	1.376
	0.2	100	0.4	500	100.	483.11	3.690	1.283
	0.3	100	0.3	500	100.	478.94	3.669	1.284
	0.4	100	0.2	500	100.	483.11	3.690	1.283

Table 12

 12 Computed apparent contact moduli for the indentation of a perfectly bonded layer with a cone.

	v 1	E 1 ⁎	v 2	E 2 ⁎	e 1 ⁎	e 2 ⁎	β 0	η
	0.2	500	0.4	100	491.	125.	0.242	1.367
	0.3	500	0.3	100	488.	99.	0.244	1.318
	0.4	500	0.2	100	492.	88.	0.248	1.296
	0.2	100	0.4	500	102.	494.	1.907	1.526
	0.3	100	0.3	500	100.	479.	2.568	1.300
	0.4	100	0.2	500	100.	485.	3.213	1.220

Acknowledgments

We would like to express our gratitude to Dr Sean Yu of the Naval Research Lab, Washington, DC, for providing us with corrected formulas for the integral equation kernel, and for valuable discussions. AC wishes to thank the Royal Society of London for the provision of a visiting fellowship to enable him to visit Oxford, and Trinity College, Oxford, for their hospitality. AMK wishes to thank the Nuffield Foundation for the provision of the award NUF-NAL.

Appendix A. Additional formulas for the elastic problem


Andrei Constantinescu
email: andrei.constantinescu@lms.polytechnique.fr

Abdelbacet Oueslati
email: abdelbacet.oueslati@univ-lille1.fr

On the reconstruction of residual stresses after matter removal in rods

Keywords: Solids and structures, Residual stresses, Rod theory, Warping Solides et structures, Contraintes résiduelles, Poutres, Gauchissement

This article discusses the reconstruction methods for the inversion of the residual stresses in rods under the Navier-Bernoulli or the Vlassov (containing warping) kinematic conditions. The proposed framework permits us to express the residual stresses in terms of the Airy stress potential expressed as the solution of an ordinary differential equation. The particular case of the Navier-Bernoulli kinematics leads to a closed-form solution.

Résumé

Sur la reconstruction de champs des contraintes résiduelles dans les barres après enlèvement de la matière. Cet article présente les méthodes de reconstruction des contraintes résiduelles dans des poutres en considérant des hypothèses cinématiques de type Navier-Bernoulli ou Vlassov (avec gauchissment). La technique proposée permet d'exprimer les contraintes résiduelles à partir d'une fonction d'Airy, solution d'une équation différentielle ordinaire. La cinématique de Navier-Bernoulli conduit à une solution explicite de ce problème.

Version française abrégée

La détermination des contraintes résiduelles dans les structures est un sujet de grande importance puisque leur présence affecte considérablement la résistance à la fatigue et à la rupture des pièces mécaniques.

Deux classes de méthodes d'identification peuvent être utilisées : les méthodes non destructives (mesures par diffraction des rayons X, par ultrasons, etc.) et destructives (enlèvement de la matière, perçage de trou, etc.). Ces dernières provoquent un relâchement de contraintes (et de déformation) conduisant ainsi à un problème inverse de reconstruction de champ de contraintes initial à partir de mesures effectuées sur la frontière du solide.

Dans cet article, nous présentons les méthodes de reconstruction des contraintes résiduelles après enlèvement de la matière dans le cas des barres élancées de section rectangulaire. L'hypothèse principale adoptée est celle d'une relaxation purement élastique des contraintes après enlèvement successif de couches d'épaisseurs h. Cette hypothèse permet de ramener la résolution du problème de reconstruction à l'étude de trois problèmes élastiques distincts (voir Fig. 2) : le premier problème P1 est celui de l'équilibre de la poutre Ω libre de charges extérieures et contenant le champ de contraintes résiduelles inconnu σ . Pour le deuxième problème noté P2, le solide est constitué par la barre Ω h obtenue à partir de Ω après enlèvement d'une couche d'épaisseur h, c'est sur cette configuration que la mesures des contraintes par diffraction des rayons X est effectuée sur la face supérieure. Enfin, P3 est le problème élastique auxiliaire traduisant l'équlibre de la barre Ω h , soumise aux tractions surfaciques (-σ • e z ) sur sa face supérieure (σ étant le champ de contraintes résiduelles à identifier). L'application du principe de superposition à ces trois problèmes conduit à la relation [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF].

Il est impératif de noter que le problème de reconstruction possède une solution explicite si le système (2) et (3), ainsi que le problème auxiliaire ont des solutions explicites. Rappelons que pour le cas d'un demi-espace élastique, la solution du problème auxiliare P3 est fournie d'une manière explicite par la solution de Boussinesq-Cerutti et le système (2) est résolu après application de la transformée de Fourrier [START_REF] Ballard | On the inversion of subsurface residual stresses from surface stress measurements[END_REF].

Pour le cas des barres de dimensions finies, une solution explicite analogue à celle de Boussinesq-Cerutti est difficille à établir. L'idée est de résoudre le problème auxiliaire dans le cadre de la théorie des poutres de Navier-Bernoulli puis celle de Vlasov et d'introduire le potentiel de contraintes d'Airy (13) pour résoudre les Éqs. ( 2) et [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF].

Des formules intégrales explicites (16) pour la reconstruction de champ de contraintes résiduelles sont obtenues pour la cinématique de Navier-Bernoulli. Elles constituent une généralisation des formules de Morre et Avans [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF] établies pour les barres (plaques) de longueurs infinies.

Introduction

Identification of residual stresses is an important task in many engineering fields such as fatigue or fracture mechanics, where their presence can significantly increase or decrease the apparent strength of mechanical components.

The problem of stress redistribution after matter removal has instructively been illustrated by Prof. H.D. Bui who has sectioned a commercial aluminium rod. The two parts were highly deformed after the cut due to the presence of residual stresses as displayed in Fig. 1. The present article is intended to bring a modest contribution to this problem.

Accurate measurements of residual stresses are performed by non-destructive techniques X-ray diffraction measurements at the surface of structures or by destructive methods: hole-drilling, layer removal, etc., if one wants to explore the residual stress field in depth of the solid. This relies essentially on measuring the stress redistribution produced by matter removal on a body containing residual stresses produced by an incompatible initial thermal or plastic strain distribution. In such a case, one supposes that matter removal leads to an elastic redistribution of stresses and measures generally the strain or the stresses on the boundary after the matter removal by means of X-ray diffraction or by neutron diffraction. This technique conducts therefore to an inverse reconstruction problem of the initial residual stress before mater removal from measurements performed after matter removal.

The study of the matter removal problem can be related to an evolution problem and uniqueness of its solution has already been proven in [START_REF] Ballard | Reconstruction d'un champ des contraintes résiduelles à partir des contraintes mesurées sur des surfaces succesives : Existence et unicité[END_REF] for bodies of arbitrary shape.

For a practical purpose, Moore and Evans [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF] proposed for the matter removal technique, in their pioneering paper, a series of exact integral formulas for the reconstruction of residual stress field in the following cases: (a) solid rod and hollow cylinder (tube) with rotationally symmetric stresses independent of the axial coordinate z; (b) flat plate with stresses varying along the thickness of the plate. These formulas apply only in cases where the residual stress field depends on just one spatial variable (corresponding to depth). However, many configurations arise in applications where the stress fields are varying along the surface of the body. Explicit formulas for such cases have been obtained only in the following cases: (i) hollow cylinder (tube) with rotationally symmetric stresses and varying along the z-coordinate (exact infinite series expression given by Nishimura [START_REF] Nishimura | On axisymmetric residual stresses in tubes with longitudinally nonuniform stress distribution[END_REF]); (ii) two and three dimensional half-spaces with stresses depending of all spatial variables (closed-form expression obtained by Ballard and Constantinescu [START_REF] Ballard | On the inversion of subsurface residual stresses from surface stress measurements[END_REF]). The main difficulties in finding closed-form solutions of this problem for bodies of different shapes can be related to: (i) the existence of a fundamental solution of Boussinesq-Cerruti type of the elastic problem for the body under discussion; and (ii) the integration of the evolution problem defined previously.

This article proposes to fill the gap between axially homogenous stresses in infinite bars (plates) and varying stresses in a semi-infinite domains by discussing the case of slender finite dimensional rectangular solids within the framework of rod theory.

Rod theories, can be interpreted in terms of power series developments of the three dimensional elastic solution in term of a small parameter, i.e. the diameter of the cross-section of the rod. Expanding the three dimensional 'exact' elastic solution up to different orders one can obtain explicit successive approximations including bending, torsion, warping and Poisson effects of the beam under study. A detailed mathematical proof can be found in Trabucho et al. [START_REF] Trabucho | Mathematical modelling of rods[END_REF].

The asymptotic expansions of the three dimensional displacement and stress field are given in the following power series of odd power terms:

u(e) = u 0 + e 2 u 2 + e 4 u 4 + h.o.t., σ (e) = σ 0 + e 2 σ 2 + e 4 σ 4 + h.o.t. ( 1
)
where h.o.t. are higher order terms and e is the diameter of the cross-section.

The main conclusion of Trabusho et al. [START_REF] Trabucho | Mathematical modelling of rods[END_REF] is that without further assumptions on the loads or geometry, the zero order term of the asymptotic expansion i.e. (u 0 , σ 0 ) coincides with the Navier-Bernoulli model and the second order terms (u 0 , σ 0 ) + e 2 (u 2 , σ 2 ) permits the justification of the Saint-Venant model for simple torsion, Timoshenko's bending equations and Vlasov's rod theory [START_REF] Vlassov | Pièces longues en voiles minces[END_REF].

This paper discusses the reconstruction formulas for the inversion of the residual stresses in rods under the Navier-Bernoulli or the Vlassov (containing warping) kinematic conditions. The proposed framework permits one to express the residual in term of the Airy stress potential expressed as the solution of an ordinary differential equation. The particular case of the Navier-Bernoulli kinematics conducts to a closed-form solution.

The reconstruction problem

Let Ω be a free rod of rectangular section and dimensions a × b × c referred to a Cartesian coordinate system (x,y,z) (see Fig. 2). The neutral fibre coincides with (Ox) and (Oy, Oz) are principle axes of inertia of the cross-Fig. 2. The three problems for the rod.

section. This simplifying assumptions are made for the clearness of the presentation and do not represent a restriction to the generality of the setting of this problem.

From Ω, layers of constant thickness h are removed continuously, leaving for each h an e wr o dΩ h . During the matter removal procedure we assume that: H1 The relaxation of the body due to matter removal is elastic.

In order to define a method of solution, as in the three dimensional case, one can split the reconstruction problem into three elastic equilibrium problems (see Fig. 2): P1 The equilibrium of the initial solid Ω with the unknown distribution plastic strain field ε p creating residual stress field to be reconstructed. This field σ satisfies the balance equations:

div σ = 0, σ xx,x + σ xz,z = 0 σ xz,x + σ zz,z = 0 (2)
P2 The equilibrium of Ω h , the rod after matter removal containing the undisturbed plastic strain field ε p .Thisisthe configuration on which surface stress measurements are performed by means of X-ray diffraction method. P3 The auxiliary problem defined on Ω h is a purely elastic problem, with a stress distribution σ h created by surface tractions on the upper surface given by:

t =-σ (x, h) • e z =-σ xz (x, h)e x -σ zz (x, h)e z
where σ is the residual stress field of problem P1. In order to recall later this dependence we shall use the following notation:

σ h (x, z) = G[t](x, z) = G[σ xz ,σ zz ](x, z), z ∈[h, c]
Under the H1 assumption of an elastic redistribution of stresses, one can now apply the superposition principle and as a consequence the following equality between the components of the residual stress components σ xx of the three problems can be established:

σ xx (x, h) = σ m xx (x, h) + σ h xx (x, h) = σ m xx (x, h) + G[σ zx ,σ zz ](x, h) (3)
This implicit relation between the components of the residual stress and the measurement σ m xx represents together with the balance equation (2), a system of equations for the unknown components of the residual stress.

The problem has a closed-form solution provided both the expression G[σ zx ,σ zz ] and the system (2) and ( 3) have a closed form solution.

In the case of the half space problem G was computed using the Boussinesq-Cerutti stress potentials and the system has been transformed into an ordinary differential equation using an integral transform which could be integrated.

In the present case of a slender domain of finite extension, the main idea is to use the expression of G defined in the framework of an imposed rod kinematics for the problem. The balance equation will be eliminated by the use of the Airy stress potential Φ to express the residual stress σ .

Reconstruction for the Navier-Bernoulli rod

Let us consider a rod with the longitudinal axes along its mean fibre. According to the Navier-Bernoulli (Love-Kirchhoff) beam theory [START_REF] Vlassov | Pièces longues en voiles minces[END_REF], the rod kinematics is defined by the following form of displacement field:

u = u(x) -w ,x (x)z e x + w(x)e z (4)
where u(x) and w(x) are the axial translation and the transversal deflexion of the cross-section respectively. '• ,x ' denotes the partial derivative with respect to x. The angle of rotation of the cross-section is -w ,x (x) and signifies that the cross-section remains orthogonal to the mean fibre during the deformation of the rod. Moreover, the cross-section is rigid in its plane. The three dimensional strain of the rod is related with the displacement field under the small strain assumption by:

ε = ε xx e x ⊗ e x = (u ,x -w ,xx )e x ⊗ e x (5)
The stress field is expressed for a linear elastic isotropic material as:

σ = ε xx (λ + 2µ)e x ⊗ e x + λ(e y ⊗ e y + e z ⊗ e z ) (6)
where λ and µ are the Lamé moduli. One can remark that the kinematic assumption neglects the shear stress, which does not imply that the shear forces do vanish. The generalized stresses for the beam theory are related to the three dimensional stress components and to displacement field through:

N(x) = S σ xx (x, z) ds = ESu ,x (x), M y (x) = S zσ xx (x, z) ds =-EIw ,xx (x) (7)
where N(x) and M y (x) stand for the resultant force and the resultant bending moment of the section of coordinate x. S denotes the area of the cross-section and I = S z 2 dS is the moment of inertia.

The auxiliary problem Ω h (0 <h c, see Fig. 3) expresses the of equilibrium the following forces and stresses:

N h ,x + f h x = 0 M h ,xx + m h y,x + f h z = 0 (8)
The exterior loads are related to the components of residual stress σ by the following expressions:

f h (x) =-bσ xz (x, h)e x -bσ zz (x, h)e z (9)
m h y (x) = bσ xz (x, h) (c -h) 2 e y (10)
Supposing that the section x = 0 is traction free, the solution of the auxiliary problem can be defined as:

σ h xx (x, z) = N h (x) b(c -h) + 12M h y (x) b(c -h) 3 z - c + h 2 (11)
Using the previous formulae on the upper surface z = h (see Fig. 3) in the system of Eqs. ( 2) and (3) one obtains:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ xx (x, z) = σ m xx (x, z) + 4 c -z x 0 σ xz (u, z) du - 6 (c -z) 2 x 0 u 0 σ zz (u ′ ,z)du ′ du ∂σ xx ∂x + ∂σ xz ∂z = 0 ∂σ xz ∂x + ∂σ zz ∂z = 0 (12)
where σ m xx is the stress supplied by the X-ray measurements and σ are the residual stresses to be reconstructed. In order to solve this system, we shall introduce the Airy function Φ(x,z):

σ xx = ∂ 2 Φ ∂z 2 ,σ xz =- ∂ 2 Φ ∂x∂z ,σ zz = ∂ 2 Φ ∂x 2 (13)
Therefore equilibrium equations in (12) are automatically satisfied and the system reduces to following ordinary differential equation:

∂ 2 Φ ∂Z 2 - 4 Z ∂Φ ∂Z + 6 Z 2 Φ = σ m xx ( 14
)
where Z = cz. The conditions Φ(0,z) = 0, ∂Φ ∂z (0,z) = 0 and ∂Φ ∂x (0,z) = 0 can be obtained by integrating the boundary conditions for the stress field.

Eq. ( 14) is a second order Euler equation and the general solution has the form:

Φ(x,Z) = A(x)Z 2 + B(x)Z 3 + Z 2 h Z σ m xx (x, z ′ ) z ′ dz ′ -Z 3 h Z σ m xx (x, z ′ ) z ′2 dz ′ (15)
Using the boundary conditions σ xz = 0 and σ zz = 0 on the upper surface Z = c as well as σ xx = 0 on the end sections one obtains A(x) = B(x) = 0. A straightforward computations gives now the expression of the initial residual stresses as a function of the measured surface stress:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ xx (x, z) = σ m xx (x, z) + 2 z 0 σ m xx (x, z ′ ) c -z ′ dz ′ -6(c -z) z 0 σ m xx (x, z ′ ) (c -z ′ ) 2 dz ′ σ xz (x, z) = 2(c -z) z 0 σ m xx,x (x, z ′ ) c -z ′ dz ′ -3(c -z) 2 z 0 σ m xx,x (x, z ′ ) (c -z ′ ) 2 dz ′ σ zz (x, z) = (c -z) 2 z 0 σ m xx,xx (x, z ′ ) c -z ′ dz ′ -(c -z) 3 z 0 σ m xx,xx (x, z ′ ) (c -z ′ ) 2 dz ′ (16)
A simple inspection of the formulae show that they still hold if the width of the bar depends of the spatial variable x.

One can remark that imposing σ xx,x = 0 conducts formally to the original Moore and Evans formulae. The resemblance is only formal because the condition σ xx,x = 0 on a rod of finite length implies that σ xx = 0. Another remark is the apparent contradiction between the Navier-Bernoulli rod kinematics imposing σ xz = 0 and the final solution where σ xz = 0. However, a close inspection shows that this is directly related to the approximation given by the rod kinematics. In order to improve the accuracy in this respect we shall next propose to enrich the rod kinematics.

Reconstruction for the Vlasov rods with warping

General Vlasov rods with warping

The kinematical enrichment proposed by Vlasov's model in [START_REF] Vlassov | Pièces longues en voiles minces[END_REF] is based on the following form of displacement field:

u = u 1 (x) + u 2 (x)z + u 3 (x)ϕ(z) e x + w(x)e z (17)
where u 1 (x), w(x), u 2 (x) and u 3 (x) are the axial translation, the vertical deflexion, the rotation angle of the crosssection and the warping amplitude, respectively. The function ϕ represents the warping, i.e. out of plane axial displacement of the section. ϕ is supposed to be a known function and satisfies a serie of orthogonality conditions. As an illustration let us recall that warping of rods with an open profile section is necessarily related to the torsion, whereas for full or closed mean profile, warping can results also from bending and tension. In this theory, the small strain tensor is expressed as:

ε = ⎛ ⎝ ε xx 0 ε xz 000 ε xz 00 ⎞ ⎠ = ⎛ ⎝ u 1,x + u 2,x z + u 3,x ϕ 0 1 2 (u 2 + u 3 ϕ ,z + w ,x ) 000 1 2 (u 2 + u 3 ϕ ,z + w ,x ) 00 ⎞ ⎠ (18)
In linear elasticity, the three dimensional stress, is obtained from Hooke's law:

σ = ⎛ ⎝ (λ + 2µ)ε xx 02 µε xz 0 λε xx 0 2µε xz 0 λε xx ⎞ ⎠ (19)
The equilibrium equations are established by use of the Principle of Virtual Work and using the orthogonality conditions verified by the warping function ϕ:

|ϕ| S = 0 and |zϕ| S = 0 with |•| S = S • ds
They are expressed as:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ E 1 Su 1,xx + f x = 0 µSu 2,x + µSw ,xx + µ|ϕ ,z | S u 3,x + f z = 0 E 1 Iu 2,xx -µSu 2 -µ|ϕ ,z | S u 3 -µSw ,x + m y = 0 E 1 |ϕ 2 | S u 3,xx -µ|ϕ ,z | S u 2 -µ|ϕ 2 ,z | S u 3 -µ|ϕ ,z | S w ,x + r ϕ = 0 (20)
with

E 1 = 1-ν (1+ν)(1-2ν) E, E
the Young modulus and ν the Poisson ratio. The external loading is given by the distributed forces and respectively moments:

f = f x (x)e x + f z (x)e z ,m y (x) = f z (x)x, r ϕ (x, z) = f x (x)ϕ(z) (21)
The generalized stresses are related to the displacement and the three dimensional stresses by the following equations:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N = S σ xx (x, z) ds = ESu 1,x M y = S zσ xx (x, z) ds = EIu 2,x B = S ϕ(z)σ xx (x, z) ds = E 1 |ϕ 2 | S u 3,x T ϕ = S σ xz (x, z)ϕ ,z (z) ds = µ|ϕ ,z | S (u 2 + w ,x ) + µ|ϕ 2 ,z | S u 3 T = S σ xz (x, z) ds = µSu 2 + µ|ϕ ,z | S u 3 + µSw ,x (22)
where I =|z 2 | S and S =|1| S . N , T and M denote the axial force, the shear force and the bending moment respectively, defined as in the Navier-Bernoulli rod theory. B and T ϕ denote the bimoment and the second shear force, respectively. As for the previous section, one can write the solution of the auxiliary problem for each Ω h (0 <h c), within the Vlasov kinematics, given by:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ N h ,x + f h x (x) = 0 T h ,x + f h z = 0 M h y,x + T h + m ′ = 0 B h ,x -T h ϕ h + r h ϕ = 0 (23)
This provides the stress solution of the auxiliary problem:

σ h xx (x, z) = N h (x) b(c -h) + 12M h y (x) b(c -h) 3 z - c + h 2 + B h (x) |ϕ 2 h | S h ϕ h z - c + h 2
It is important to note that the warping functions are now parametrized by the height of the removed layer h and therefore denoted as ϕ h .

The system of Eqs. ( 2) is therefore reduced to the differential equation for the unknown residual stress field:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ xx (x, z) = σ m xx (x, z) + 4 c -z x 0 σ xz (u, z) du - 6 (c -z) 2 x 0 u 0 σ zz (u ′ ,z)du ′ du + ϕ z ( z-c 2 ) |ϕ 2 z | S z B z (x) ∂σ xx ∂x + ∂σ xz ∂z = 0 ∂σ xz ∂x + ∂σ zz ∂z = 0 (24)
As before, the expression of the stress field using an Airy potential function Φ(x,z) eliminates the equilibrium equations in (24) leading to the following expression of the bimoment:

B z (x) = b c z σ xx (x, z ′ )ϕ(z ′ ) dz = b c z ∂ 2 Φ ∂z 2 (x, z ′ )ϕ(z ′ ) dz ′
The system (24) becomes:

∂ 2 Φ ∂z 2 - 4 c -z ∂Φ ∂z + 6 (c -z) 2 Φ -b ϕ z ( z-c 2 ) |ϕ z | S z c z ∂ 2 Φ ∂z 2 (x, z ′ )ϕ z (z ′ ) dz ′ = σ m xx (x, c -z) (25)
Using the boundary conditions Φ(0,z)= 0 and ∂Φ ∂z (0,z)= 0, one finally obtains:

∂ 2 Φ ∂z 2 - 4 c -z ∂Φ ∂z + 6 (c -z) 2 Φ -b ϕ z ( z-c 2 ) |ϕ 2 z )| S h c z Φ(x,z ′ ) ∂ 2 ϕ(z ′ ) ∂z 2 dz ′ = σ m xx (x, z). (26)
Unfortunately, this integral-differential equation does not have a closed-form solution for a general form of the warping function ϕ. From the practical point of view, one has to choose first the warping function ϕ adapted to the rod profile and integrate afterwards Eq. (26) symbolically or numerically.

Vlasov rod with warping and Navier-Bernoulli assumption

An interesting result is obtained for the particular case of Vlasov kinematics with wrapping combined with the Navier-Bernoulli assumption of orthogonality of the cross-section and the mean fibre during deformation. This condition is expressed as: u 2 (x) =-w ,x (x), and the general form of the displacement is:

u = u 1 (x) -w ,x (x)z + u 3 (x)ϕ(z) e x + w(x)e z (27)
The second shear force T ϕ is in this case related to the shear force T by:

T ϕ = |ϕ 2 ,z | S |ϕ ,z | S T
The expression of generalized stresses is reduced to:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N(x) = b x 0 σ xz (u, h) du M(x) = b x 0 u 0 σ zz (u ′ ,h)du ′ du -b c -h 2 x 0 σ xz (u, h) du B(x) =-b |ϕ 2 ,z | |ϕ ,z | x 0 u 0 σ zz (u ′ ,h)du ′ du + bϕ h -c 2 x 0 σ xz (u, h) du (28)
Under the kinematic assumptions of this theory, the expression of the solution of the auxiliary problem σ h xx becomes:

σ h xx (x, h) = α h x 0 σ xz (u, h) du -β h x 0 u 0 σ zz (u ′ ,h)du ′ du where α h = 4 c -h + b |ϕ 2 z | S h ϕ 2 h h -c 2 ,β h = 6 (c -h) 2 + bϕ h h -c 2 |ϕ 2 h ,z | S h |ϕ h ,z | S h
The system (24) is now only reduced to:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ xx (x, z) = σ m xx (x, z) + α z x 0 σ xz (u, z) du -β z x 0 u 0 σ zz (u ′ ,z)du ′ du ∂σ xx ∂x + ∂σ xz ∂z = 0 ∂σ xz ∂x + ∂σ zz ∂z = 0 (29)
Again introducing Airy Φ(x,z) stress potential we obtain the following differential equation:

∂ 2 Φ ∂z 2 -α z ∂Φ ∂z + β z Φ = σ m xx (30)
with the boundary conditions Φ(0,z)= 0, ∂Φ ∂z (0,z)= 0 and ∂Φ ∂x (0,z)= 0. This second order partial differential equation can generally be solved in a symbolic way, provided the expression of the functions α z and β z are explicit. This can be achieved by a particular choice of the warping function ϕ.

Conclusion

The aim of this article was to investigate the problem of reconstructing residual stress from measurements after matter removal in slender bodies. The technique used was based on the superposition principle applied to the initial domain and the domain after matter removal. For all studies rod kinematics the residual stresses were directly related to the solution of an ordinary differential equation on vertical lines of the cross-sections of the rod.

An interesting point is to compare the information necessary to compute the initial residual stress field in an inner point of the body:

-for an infinite rod (plate), without stress gradient along the surface (Moore and Evans [START_REF] Moore | Mathematical correction for stress in removed layers in X-ray diffraction residual stress analysis[END_REF]): one needs to measure the stresses on the vertical line above the point; -for finite rod (plate) with stress gradient along the surface (this paper): one needs to measure the stresses and their gradient on the vertical line above the point; -half-space without stress gradient along the surface (Ballard and Constantinescu, 1994 [2]): one needs to measure the stresses in the all the removed volume.

This reflects on the one hand the growing complexity of the problem and on the other hand the strength of the rod kinematics hypothesis. One should, however, remember that the rod solution is only an approximate solution, when compared to the fully three dimensional solution of the problem. In order to construct this solution one should dispose of the fundamental stress solution for a rectangular body of finite extension, which is not available for the moment.

Fig. 1 .

 1 Fig. 1. Deformation of an aluminium rod after cutting, due to the presence of residual stresses, an experiment proposed by Prof. H.D. Bui.

Fig. 3 .

 3 Fig. 3. The coordinates for the rod after removal of a layer of thickness h.