Laurence Brassart

Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity from a fast explicit method

Keywords: laminate, anisotropic elasticity, incompatibility stresses, effective behavior, composite structure

Introduction

At a material surface of discontinuity like a grain boundary or an interface between two dissimilar materials, incompatibility stresses may arise due to heterogeneous thermo-elastic properties and/or difference of plastic strains. These internal stresses can be computed numerically thanks to finite element methods [START_REF] Roters | Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications[END_REF] or spectral approaches (A. [START_REF] Lebensohn | N-site modeling of a 3D viscoplastic polycrystal using Fast Fourier Transform[END_REF]. Analytical formulations also exist in the case of two-phase laminates [START_REF] Stupkiewicz | Modelling of laminated microstructures in stress-induced martensitic transformations[END_REF][START_REF] Milton | The Theory of Composites[END_REF][START_REF] Franciosi | Heterogeneous crystal and poly-crystal plasticity modeling from a transformation field analysis within a regularized Schmid law[END_REF][START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF] which can also be viewed as infinite bi-crystals or bi-materials subjected to a homogeneous far-field loading [START_REF] Gemperlová | Stresses in deformed cubic bicrystals with a tilt ⟨001⟩ grain boundary[END_REF][START_REF] Gemperlová | Compatibility stresses in deformed bicrystals[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF][START_REF] Richeton | Analytical expressions of incompatibility stresses at Σ3⟨111⟩ twin boundaries and consequences on single-slip promotion parallel to twin plane[END_REF][START_REF] Richeton | Incompatibility stresses and lattice rotations due to grain boundary sliding in heterogeneous anisotropic elasticity[END_REF]. Such analytical models provide instantaneous responses which are first-order estimations around interfaces subjected to more complex boundary conditions [START_REF] Peralta | Elastic stresses in anisotropic bicrystals[END_REF][START_REF] Tiba | Incompatibility stresses at grain boundaries in Ni bicrystalline micropillars analyzed by an anisotropic model and slip activity[END_REF]. Hence, it can be very useful to test, on a huge number of configurations, the influence of different parameters like material single crystal elastic constants, crystallographic misorientation, crystal volume fraction, interface inclination, amount of plastic deformation or loading type [START_REF] Gemperlová | Compatibility stresses in deformed bicrystals[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF][START_REF] Richeton | Analytical expressions of incompatibility stresses at Σ3⟨111⟩ twin boundaries and consequences on single-slip promotion parallel to twin plane[END_REF]. In the same way, explicit analytical models for N-phase laminates are very useful to investigate the influence of the previously mentioned parameters in problems involving several parallel interfaces, like in thin films (e.g., [START_REF] Welzel | The determination of stresses in thin films; modelling elastic grain interaction[END_REF][START_REF] Abadias | Properties of the Eshelby tensor and existence of the equivalent ellipsoidal inclusion solution[END_REF], in electronics packaging (e.g., [START_REF] Wong | A More Comprehensive Solution for Tri-Material Layers Subjected to Thermal Stress[END_REF], in nanolaminated materials (e.g., [START_REF] Mukhopadhyay | Effective mechanical properties of multilayer nano-heterostructures[END_REF][START_REF] Wang | Strength and plasticity of nanolaminated materials[END_REF] or in materials with welding joints which are classically modelled as a base metal, a heat affected zone and a weld metal (e.g., [START_REF] Du | Effect of rolling process on fatigue performance of stir zone of AA6061-T6 double-side friction stir welding joint[END_REF]. Actually, laminated materials are widely used in the industry because of the relative ease to control their properties from the choice of the materials combination,

Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity volume fraction and orientation [START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF]. For the design of laminates, fast estimates of their properties is thus required.

The objective of the present paper is to derive an explicit method, which can be easily implemented in a numerical code, in order to compute the effective behavior and the distribution of stresses in a general N-phase laminate made of parallel, planar and perfectly bonded interfaces and which is assumed to be subjected to a far-field homogeneous loading and/or homogeneous temperature variation, in addition to the possibility of having different uniform plastic strains in each phase. The method should be exact for a small number of phase (𝑁 = 2, 3, 4, ...) as well as for an infinitely large number of phases. Such an objective may appear relatively standard considering previous works on laminates in the literature (e.g., [START_REF] Milton | The Theory of Composites[END_REF][START_REF] Franciosi | Heterogeneous crystal and poly-crystal plasticity modeling from a transformation field analysis within a regularized Schmid law[END_REF][START_REF] El Omri | Elastic-plastic homogenization for layered composites[END_REF][START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF][START_REF] Glüge | Effective plastic properties of laminates made of isotropic elastic plastic materials[END_REF]). However, it is generally hard to find compact explicit formulas for direct use in engineering applications. Besides, at the difference of previous approaches where explicit formulas were restricted to a finite number of phases [START_REF] Milton | The Theory of Composites[END_REF][START_REF] Franciosi | Heterogeneous crystal and poly-crystal plasticity modeling from a transformation field analysis within a regularized Schmid law[END_REF] or to isotropic elasticity [START_REF] El Omri | Elastic-plastic homogenization for layered composites[END_REF][START_REF] Milton | The Theory of Composites[END_REF][START_REF] Glüge | Effective plastic properties of laminates made of isotropic elastic plastic materials[END_REF] or to purely elastic behavior [START_REF] Milton | The Theory of Composites[END_REF][START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF], this contribution provides explicit compact formulas for an arbitrary number of phases, a completely general anisotropic elasticity as well as different thermal and plastic strains in the phases. In particular, a simple direct analytical formula is derived in order to compute the stress in a given phase once the effective behavior of the laminate has been determined. Furthermore, the code to compute the effective behavior and the stress distribution in an arbitrary N-phase laminate is made freely available (see Supplementary Material at the end of the document).

The paper is organized as follows. Section 2 reminds the expressions of incompatibility stresses in a two-phase laminate from which the effective compliance, effective plastic strain and effective thermal expansion tensors can be deduced. Then, Section 3 presents an explicit method to find the effective behavior and the stress distribution for a general N-phase laminate based on the solutions of the two-phase laminate. Section 4 is dedicated to the validation of the proposed method by comparisons with finite element simulation results on a same boundary value problem. In Section 5, an application of the method is performed for a near-𝛽 titanium alloy with elongated grains, by comparing the level of internal stresses for different elastic loadings within a N-phase laminate made of 100 000 orientations and a 2-phase laminate with maximal elastic contrast. Then, in Section 6, comparisons with elastic self-consistent models considering oblate spheroidal grains of different aspect ratios are performed. In particular, the numerical performance of the laminate model compared to the self-consistent models is discussed. Finally, concluding remarks are given in Section 7.

In the following, the Einstein summation convention over repeated indices and the contracted Voigt notation (Voigt 1966) (11 → 1, 22 → 2, 33 → 3, 23 → 4, 31 → 5, 12 → 6) are used, i.e. indices range from 1 to 6. For consistency, an engineering convention is considered for strain components, i.e., 𝜀 4 = 2𝜀 23 , 𝜀 5 = 2𝜀 31 , 𝜀 6 = 2𝜀 12 , while the components of the elastic compliance tensor 𝑠 𝑖 𝑗 include the multiplying factors of 2 and 4 [START_REF] Voigt | Lehrbuch der Kristallphysik[END_REF]. Besides, the notation ⟨⟩ denotes a volume average and superscripts 1, 2 or 𝑚 denote fields in materials 1, 2 or 𝑚, respectively.

2 Two-phase laminate

Averaging rules, constitutive behavior and continuity conditions

In this section, an infinite 2-phase laminate composed of two alternating materials or crystals perfectly bonded along planar interfaces is considered. This rank-1 laminate is supposed to have been deformed under the action of a homogeneous temperature variation Δ𝑇 , as well as a macroscopic homogeneous and remotely applied stress Σ 𝑖 satisfying the averaging rule

Σ 𝑖 = ⟨𝜎 𝑖 ⟩, (1)
where 𝜎 𝑖 is the local Cauchy stress. Linear thermo-elasticity is assumed, as well as the possibility of having different uniform plastic strains 𝜀 𝑝 𝑖 in each phase. Under these assumptions and through consideration of strains compatibility and balance of linear momentum without body force in a static small strain setting, it can be shown that 𝜎 𝑖 and the local strains 𝜀 𝑖 are uniform in each

Thiebaud

Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity material [START_REF] Gemperlová | Compatibility stresses in deformed bicrystals[END_REF][START_REF] El Omri | Elastic-plastic homogenization for layered composites[END_REF][START_REF] Milton | The Theory of Composites[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF], so that

Σ 𝑖 = 𝑓 1 𝜎 1 𝑖 + 𝑓 2 𝜎 2 𝑖 , (2)
and

𝐸 𝑖 = ⟨𝜀 𝑖 ⟩ = 𝑓 1 𝜀 1 𝑖 + 𝑓 2 𝜀 2 𝑖 , (3)
where 𝑓 1 and 𝑓 2 = 1 -𝑓 1 are the material volume fractions. 𝐸 𝑖 is defined as the macroscopic strain. The strains in each material are expressed as

𝜀 𝑚 𝑖 = 𝑠 𝑚 𝑖 𝑗 𝜎 ℓ 𝑗 + 𝜀 𝑝 𝑚 𝑖 + 𝛼 𝑇 𝑚 𝑖 Δ𝑇 , (4)
where 𝑚 = 1, 2 while 𝑠 𝑖 𝑗 and 𝛼 𝑇 𝑖 are the components of the elastic compliance tensor and of the symmetrical thermal expansion tensor, respectively. Both 𝑠 𝑖 𝑗 and 𝛼 𝑇 𝑖 are supposed isotherm and uniform in each material. Besides, considering that the interface normals are parallel to the direction 𝒆 2 of a Cartesian frame (𝒆 1 , 𝒆 2 , 𝒆 3 ), Equations ( 2) and (3) along with tractions continuity and strains compatibility impose the following continuity conditions [START_REF] El Omri | Elastic-plastic homogenization for layered composites[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF][START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF])

         𝜎 1 2 = 𝜎 2 2 = Σ 2 , 𝜎 1 4 = 𝜎 2 4 = Σ 4 , 𝜎 1 6 = 𝜎 2 6 = Σ 6 , and
         𝜀 1 1 = 𝜀 2 1 = 𝐸 1 , 𝜀 1 3 = 𝜀 2 3 = 𝐸 3 , 𝜀 1 5 = 𝜀 2 5 = 𝐸 5 .
(5)

Stress partitioning

It was previously shown that the stresses in both materials can be written as [START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF][START_REF] Richeton | Analytical expressions of incompatibility stresses at Σ3⟨111⟩ twin boundaries and consequences on single-slip promotion parallel to twin plane[END_REF][START_REF] Tiba | Incompatibility stresses at grain boundaries in Ni bicrystalline micropillars analyzed by an anisotropic model and slip activity[END_REF][START_REF] Richeton | Incompatibility stresses and lattice rotations due to grain boundary sliding in heterogeneous anisotropic elasticity[END_REF])

𝜎 1 𝑖 = Σ 𝑖 + 𝑓 2 𝐺 𝑖𝑘 [(𝑠 2 𝑘 𝑗 -𝑠 1 𝑘 𝑗 )Σ 𝑗 + 𝜀 𝑝2 𝑘 -𝜀 𝑝1 𝑘 + (𝛼 𝑇 2 𝑘 -𝛼 𝑇 1 𝑘 )Δ𝑇 ], 𝜎 2 𝑖 = Σ 𝑖 -𝑓 1 𝐺 𝑖𝑘 [(𝑠 2 𝑘 𝑗 -𝑠 1 𝑘 𝑗 )Σ 𝑗 + 𝜀 𝑝2 𝑘 -𝜀 𝑝1 𝑘 + (𝛼 𝑇 2 𝑘 -𝛼 𝑇 1 𝑘 )Δ𝑇 ], (6)
where the non-zero components of the symmetric tensor 𝐺 𝑖 𝑗 are given by

𝐺 11 = ( s2 35 -s33 s55 )/𝐷 , 𝐺 33 = ( s2 15 -s11 s55 )/𝐷 , 𝐺 55 = ( s2 13 -s11 s33 )/𝐷 , 𝐺 13 = ( s13 s55 -s15 s35 )/𝐷 , 𝐺 15 = ( s15 s33 -s13 s35 )/𝐷 , 𝐺 35 = ( s11 s35 -s13 s15 )/𝐷 , (7)
with 𝐷 = s11 s2 35 + s33 s2 15 + s55 s2 13 -s11 s33 s55 -2s 13 s15 s35 and s𝑖 𝑗 = 𝑓 2 𝑠 1 𝑖 𝑗 + 𝑓 1 𝑠 2 𝑖 𝑗 . It is noteworthy that these expressions can also be viewed as the solutions in an infinite bi-material [START_REF] Gemperlová | Compatibility stresses in deformed bicrystals[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF][START_REF] Richeton | Analytical expressions of incompatibility stresses at Σ3⟨111⟩ twin boundaries and consequences on single-slip promotion parallel to twin plane[END_REF] and that the effect grain boundary sliding could be considered as well [START_REF] Richeton | Incompatibility stresses and lattice rotations due to grain boundary sliding in heterogeneous anisotropic elasticity[END_REF]). The differences 𝜎 1 𝑖 -Σ 𝑖 and 𝜎 2 𝑖 -Σ 𝑖 are called incompatibility stresses in materials 1 and 2, respectively. Besides, it must be noted that the tensor 𝐺 𝑖 𝑗 is a function of 𝑠 1 𝑖 𝑗 , 𝑠 2 𝑖 𝑗 and 𝑓 1 (or 𝑓 2 ) only so that we will use the notation 𝐺 𝑖 𝑗 (𝑠 1 , 𝑠 2 , 𝑓 1 ) in Section 3 where multiple tensors 𝐺 𝑖 𝑗 are considered.

Effective behavior

By definition of the effective behavior, the macroscopic strain satisfies the following relation

𝐸 𝑖 = 𝑆 𝑖 𝑗 Σ 𝑗 + 𝐸 𝑝 𝑖 + 𝐴 𝑇 𝑖 Δ𝑇 , (8)
where 𝑆 𝑖 𝑗 is the effective compliance tensor, 𝐸 𝑝 𝑖 the effective plastic strain tensor and 𝐴 𝑇 𝑖 the effective thermal expansion tensor. From Equations (3), ( 4) and ( 6), we have also

𝐸 𝑖 =⟨𝑠 𝑖 𝑗 𝜎 𝑗 ⟩+⟨𝜀 𝑝 𝑖 ⟩+⟨𝛼 𝑇 𝑖 ⟩Δ𝑇 =⟨𝑠 𝑖 𝑗 ⟩Σ 𝑗 -𝑓 1 𝑓 2 (𝑠 2 𝑖 𝑗 -𝑠 1 𝑖 𝑗 )𝐺 𝑗𝑘 ((𝑠 2 𝑘ℓ -𝑠 1 𝑘ℓ )Σ ℓ +𝜀 𝑝2 𝑘 -𝜀 𝑝1 𝑘 + (𝛼 𝑇 2 𝑘 -𝛼 𝑇 1 𝑘 )Δ𝑇 ) +⟨𝜀 𝑝 𝑖 ⟩+⟨𝛼 𝑇 𝑖 ⟩Δ𝑇 . ( 9
)
Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity

Then, by identification with Equation ( 8), the three effective tensors can be expressed as

𝑆 𝑖 𝑗 = ⟨𝑠 𝑖 𝑗 ⟩ -𝑓 1 𝑓 2 (𝑠 2 𝑖𝑘 -𝑠 1 𝑖𝑘 )𝐺 𝑘ℓ (𝑠 2 ℓ 𝑗 -𝑠 1 ℓ 𝑗 ) (10) 𝐸 𝑝 𝑖 = ⟨𝜀 𝑝 𝑖 ⟩ -𝑓 1 𝑓 2 (𝑠 2 𝑖𝑘 -𝑠 1 𝑖𝑘 )𝐺 𝑘ℓ (𝜀 𝑝2 ℓ -𝜀 𝑝1 ℓ ) (11) 𝐴 𝑇 𝑖 = ⟨𝛼 𝑇 𝑖 ⟩ -𝑓 1 𝑓 2 (𝑠 2 𝑖𝑘 -𝑠 1 𝑖𝑘 )𝐺 𝑘ℓ (𝛼 𝑇 2 ℓ -𝛼 𝑇 1 ℓ ). ( 12
)
3 N-phase laminate

In the following, an infinite N-phase laminate composed of 𝑁 alternating materials perfectly bonded along parallel planar interfaces is considered. All the other assumptions made for the 2-phase laminate in Section 2, including Equation (1) and Equation ( 4), are maintained so that stresses and strains are uniform within each phase (Glüge and Kalisch 2014).

Effective behavior

The solutions of the 2-phase laminate are first used to find the effective behavior of the N-phase laminate and then the stresses in each phase. The general methodology is described in Figure 1. It

prop 1,1 prop 1,2 prop 1,3 prop 1,4 . . . prop 1,𝑁 -3 prop 1,𝑁 -2 prop 1,𝑁 -1 prop 1,𝑁 ℓ = 1, 𝐾 1 = 𝑁 ℓ = 2, 𝐾 2 = 𝑁 /2 or (𝑁 + 1)/2 ℓ = 𝑛 -1, 𝐾 𝑛-1 = 2 ℓ = 𝑛, 𝐾 𝑛 = 1 prop 2,1 prop 2,2 prop 2,𝐾 2 -1 prop 2,𝐾 2 prop 𝑛-1,1 prop 𝑛-1,1 prop 𝑛,1
Figure 1 Description of the multiple steps procedure for determining the effective properties of a N-phase laminate and the stresses in each phase: ℓ denotes the step number and 𝐾 ℓ the number of phases at step ℓ. In every phase, the properties are denoted prop ℓ,𝑚 = (𝑓 ℓ,𝑚 , 𝑠 ℓ,𝑚 , 𝜀 𝑝 ℓ,𝑚 , 𝛼 𝑇 ℓ,𝑚 ).

consists in a multiple steps procedure, where at each step ℓ, the phases are grouped two by two and replaced by a new phase having the effective properties of a 2-phase laminate made of the corresponding two phases. The procedure is followed until only one phase is remaining, corresponding to

ℓ = 𝑛 = ⌈log 2 𝑁 + 1⌉. ( 13
)
The number of phases or materials of the laminate at a step ℓ is denoted 𝐾 ℓ . Hence, we have 𝐾 ℓ = 𝑁 , 𝐾 𝑛 = 1 and

𝐾 ℓ+1 = 𝐾 ℓ /2 if 𝐾 ℓ even (𝐾 ℓ + 1)/2 otherwise. ( 14
)
At step ℓ + 1, the properties of a new phase 𝑚, i.e. the volume fractions 𝑓 ℓ+1,𝑚 , elastic compliances 𝑠 ℓ+1,𝑚 𝑖 𝑗

, plastic strains 𝜀 𝑝 ℓ +1,𝑚 𝑖 and thermal expansions 𝛼 𝑇 ℓ +1,𝑚

𝑖

, are computed based on Equations ( 10)

Thiebaud Richeton
Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity to (12):

𝑓 ℓ+1,𝑚 =𝑓 ℓ,2𝑚-1 +𝑓 ℓ,2𝑚 𝑠 ℓ+1,𝑚 𝑖 𝑗 =𝑓 * 𝑠 ℓ,2𝑚-1 𝑖 𝑗 +(1-𝑓 * )𝑠 ℓ,2𝑚 𝑖 𝑗 -𝑓 * (1-𝑓 * ) (𝑠 ℓ,2𝑚 𝑖𝑘 -𝑠 ℓ,2𝑚-1 𝑖𝑘 )𝐺 𝑘ℓ (𝑠 2𝑚-1 ,𝑠 2𝑚 ,𝑓 * ) (𝑠 ℓ,2𝑚 ℓ 𝑗 -𝑠 ℓ,2𝑚-1 ℓ 𝑗 ) 𝜀 𝑝 ℓ +1,𝑚 𝑖 =𝑓 * 𝜀 𝑝 ℓ,2𝑚-1 𝑖 +(1-𝑓 * )𝜀 𝑝 ℓ,2𝑚 𝑖 -𝑓 * (1-𝑓 * ) (𝑠 ℓ,2𝑚 𝑖𝑘 -𝑠 ℓ,2𝑚-1 𝑖𝑘 )𝐺 𝑘ℓ (𝑠 2𝑚-1 ,𝑠 2𝑚 ,𝑓 * ) (𝜀 𝑝 ℓ,2𝑚 ℓ -𝜀 𝑝 ℓ,2𝑚-1 ℓ ) 𝛼 𝑇 ℓ +1,𝑚 𝑖 =𝑓 * 𝛼 𝑇 ℓ,2𝑚-1 𝑖 +(1-𝑓 * )𝛼 𝑇 ℓ,2𝑚 𝑖 -𝑓 * (1-𝑓 * ) (𝑠 ℓ,2𝑚 𝑖𝑘 -𝑠 ℓ,2𝑚-1 𝑖𝑘 )𝐺 𝑘ℓ (𝑠 2𝑚-1 ,𝑠 2𝑚 ,𝑓 * ) (𝛼 𝑇 ℓ,2𝑚 ℓ -𝛼 𝑇 ℓ,2𝑚-1 ℓ ), (15)
with 𝑓 * = 𝑓 ℓ,2𝑚-1 /(𝑓 ℓ,2𝑚-1 + 𝑓 ℓ,2𝑚 ). By imposing the consistency condition

𝑓 * 𝜎 ℓ,2𝑚-1 𝑖 + (1 -𝑓 * )𝜎 ℓ,2𝑚 𝑖 = 𝜎 ℓ+1,𝑚 𝑖 (16)
between stresses at step ℓ and at step ℓ + 1, the stresses at step ℓ can be related to the stresses at step ℓ + 1 by application of Equation ( 6):

𝜎 ℓ,2𝑚-1 𝑖 = 𝜎 ℓ+1,𝑚 𝑖 + (1 -𝑓 * )𝐺 𝑖𝑘 (𝑠 2𝑚-1 , 𝑠 2𝑚 , 𝑓 * )𝛽 𝜎 ℓ,2𝑚-1 𝑖 = 𝜎 ℓ+1,𝑚 𝑖 -𝑓 * 𝐺 𝑖𝑘 (𝑠 2𝑚-1 , 𝑠 2𝑚 , 𝑓 * )𝛽 (17) with 𝛽 = (𝑠 ℓ,2𝑚 𝑘 𝑗 -𝑠 ℓ,2𝑚-1 𝑘 𝑗 )𝜎 ℓ+1,𝑚 𝑗 + 𝜀 𝑝 ℓ,2𝑚 𝑘 -𝜀 𝑝 ℓ,2𝑚-1 𝑘 + (𝛼 𝑇 ℓ,2𝑚 𝑘 -𝛼 𝑇 ℓ,2𝑚-1

𝑘

)Δ𝑇 . In case where 𝐾 ℓ is odd, the last phase of step ℓ + 1 is simply set to be the same as the last phase of step ℓ, i.e.

𝑓 ℓ+1,𝐾 ℓ +1 = 𝑓 ℓ,𝐾 ℓ , 𝑠 ℓ+1,𝐾 ℓ +1 𝑖 𝑗 = 𝑠 ℓ,𝐾 ℓ 𝑖 𝑗 , 𝜀 𝑝 ℓ +1,𝐾 ℓ +1 𝑖 = 𝜀 𝑝 ℓ,𝐾 ℓ 𝑖 , 𝛼 𝑇 ℓ +1,𝐾 ℓ +1 𝑖 = 𝛼 𝑇 ℓ,𝐾 ℓ 𝑖 , 𝜎 ℓ+1,𝐾 ℓ +1 𝑖 = 𝜎 ℓ,𝐾 ℓ 𝑖 . (18)
From the procedure just described, it is noticeable that averaging rules and continuity conditions are automatically satisfied at each step, i.e.

∀ℓ,

𝑚 𝑓 ℓ,𝑚 = 1, ⟨𝜎 ℓ 𝑖 ⟩ = Σ 𝑖 , ⟨𝜀 ℓ 𝑖 ⟩ = 𝐸 𝑖 (19) together with ∀ℓ, 𝑚, 𝜎 ℓ,𝑚 2 = Σ 2 , 𝜎 ℓ,𝑚 4 = Σ 4 , 𝜎 ℓ,𝑚 6 = Σ 6 , (20)
and

∀ℓ, 𝑚, 𝜀 ℓ,𝑚 1 = 𝐸 1 , 𝜀 ℓ,𝑚 3 = 𝐸 3 , 𝜀 ℓ,𝑚 5 = 𝐸 5 . ( 21
)
As a consequence, the remaining phase at step 𝑛 satisfies

𝐸 𝑖 = 𝑠 𝑛,1 𝑖 𝑗 Σ 𝑗 + 𝜀 𝑝 𝑛,1 𝑖 + 𝛼 𝑇 𝑛,1 𝑖 Δ𝑇 (22)
which means, from Equation ( 8), that the effective tensors of the N-phase laminate are

𝑆 𝑖 𝑗 = 𝑠 𝑛,1 𝑖 𝑗 , 𝐸 𝑝 𝑖 = 𝜀 𝑝 𝑛,1
𝑖 , and

𝐴 𝑇 𝑖 = 𝛼 𝑇 𝑛,1 𝑖 . ( 23
)
In the above procedure which consists in grouping the phases in pair iteratively until one phase with the effective properties remains, the initial relative positions of the different phases have obviously no influence on the values obtained at the final step. This means that the effective properties of the laminate are not sensitive to the stacking order of the phases as it was already noted by [START_REF] Glüge | The effective stiffness and stress concentrations of a multi-layer laminate[END_REF].

Stress partitioning

The stresses in each phase of the initial N-phase laminate can be determined from a top-down algorithm, starting from 𝜎 𝑛,1 𝑖 = Σ 𝑖 at step 𝑛 and then going down to step 1 by application of Equation ( 17) step by step. However, a direct and much faster computation is possible by first noticing that for a 2-phase laminate

Σ 𝑖 + (1 -𝑓 1 )𝐺 𝑖𝑘 (𝑠 1 , 𝑠 2 , 𝑓 1 ) [(𝑠 2 𝑘 𝑗 -𝑠 1 𝑘 𝑗 )Σ 𝑗 + 𝜀 𝑝2 𝑘 -𝜀 𝑝1 𝑘 + (𝛼 𝑇 2 𝑘 -𝛼 𝑇 1 𝑘 )Δ𝑇 ] = Σ 𝑖 + 𝐺 𝑖𝑘 (𝑠 1 , 𝑆, 0) [(𝑆 𝑘 𝑗 -𝑠 1 𝑘 𝑗 )Σ 𝑗 + 𝐸 𝑝 𝑘 -𝜀 𝑝1 𝑘 + (𝐴 𝑇 𝑘 -𝛼 𝑇 1 𝑘 )Δ𝑇 ], (24)
Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity

where 𝑆 𝑖 𝑗 , 𝐸 𝑝 𝑖 , 𝐴 𝑇 𝑖 are, in this equation, the effective properties of the 2-phase laminate. Indeed, let us consider a fictitious cut of phase 1 parallel to the planar interface, such that one of two remaining parts, denoted hereafter 1 * , has a negligible volume fraction (𝑓 1 * → 0). Stresses in phase 1 * are equal to those in the first phase due to uniformity of stresses within a phase (see Section 2). This corresponds to the first member of Equation (24) (cf. Equation ( 6)). Then, as a consequence of the preceding procedure, it is known that these stresses can also be computed from a 2-phase laminate composed of phase 1 * and another phase of volume fraction equals to 1 which has thus necessarily the effective properties of the initial 2-phase laminate. This corresponds to the second member of Equation ( 24). The principle is illustrated in Figure 2. to establish Equation ( 24).

𝑓 1 , 𝑠 1 , 𝜀 𝑝 1 , 𝛼 𝑇 1 , 𝜎 1 𝑓 2 , 𝑠 2 , 𝜀 𝑝 2 , 𝛼 𝑇 2 , 𝜎 2 𝑓 → 0, 𝑠 1 , 𝜀 𝑝 1 , 𝛼 𝑇 1 , 𝜎 1 𝑓 = 1, 𝑆, 𝐸 𝑝 , 𝐴 𝑇 , Σ
Such a thought experiment can of course also be made for the N-phase laminate, i.e. performing a fictitious cut in a given phase 𝑚 such that one of the remaining part has a negligible volume fraction and then considering a 2-phase laminate made of the phase 𝑚 * with negligible volume fraction and another phase having the effective properties of the N-phase laminate. Accordingly, once the effective properties (𝑆 𝑖 𝑗 , 𝐸 𝑝 𝑖 , 𝐴 𝑇 𝑖 ) of the the N-phase laminate are known, the stresses in each phase can be directly computed as

𝜎 𝑚 𝑖 = Σ 𝑖 + 𝐺 𝑖𝑘 (𝑠 𝑚 , 𝑆, 0) [(𝑆 𝑘 𝑗 -𝑠 𝑚 𝑘 𝑗 )Σ 𝑗 + 𝐸 𝑝 𝑘 -𝜀 𝑝 𝑚 𝑘 + (𝐴 𝑇 𝑘 -𝛼 𝑇 𝑚 𝑘 )Δ𝑇 ]. (25)
Equation ( 25) shows that the stresses in one phase can be expressed only with respect to the effective properties of the laminate and the properties of the phase. Since the effective properties of the laminate are independent of the relative positions of the phases (cf. Section 3.1), this means that the stresses in the phases are also not sensitive to the stacking order of the phases. Besides, the numerical implementation of the proposed method is very straightforward and is provided as a Matlab code (see Supplementary materials). It contains two functions, one that computes the tensor 𝐺 𝑖 𝑗 (cf. Equation ( 7)) from two elastic compliance tensors and a volume fraction value and one that computes the effective properties of a 2-phase laminate according to Equation (15). The main code is just composed of a single iterative loop that calls the two aforementioned functions in order to compute the properties of the new phases step by step as described in Figure 1. At step 𝑛, the effective behavior of the N-phase laminate is known through Equation ( 23) and the stresses in any phase can then be directly computed from Equation ( 25).

Finally, it may be worth to indicate also the simple expressions of the stresses in the particular case of isotropic homogeneous thermo-elasticity. With 𝜇 the shear modulus and 𝜈 the Poisson's ratio, the plane components are indeed given by [START_REF] Rey | Slip heterogeneities in deformed aluminum bicrystals[END_REF][START_REF] Richeton | Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint[END_REF])

𝜎 𝑚 1 = Σ 1 + 2𝜇 1 -𝜈 [⟨𝜀 𝑝 1 ⟩ -𝜀 𝑝 𝑚 1 + 𝜈 (⟨𝜀 𝑝 3 ⟩ -𝜀 𝑝 𝑚
3 )],

𝜎 𝑚 3 = Σ 3 + 2𝜇 1 -𝜈 [⟨𝜀 𝑝 3 ⟩ -𝜀 𝑝 𝑚 3 + 𝜈 (⟨𝜀 𝑝 1 ⟩ -𝜀 𝑝 𝑚 1 )], 𝜎 𝑚 5 = Σ 5 + 𝜇 [⟨𝜀 𝑝 5 ⟩ -𝜀 𝑝 𝑚 5 ], (26)
while the anti-plane components are given by Equation ( 20).

Validation by finite element simulations

Several Finite Element (FE) simulations were performed in order to check the correctness of the methodology described in Section 3. The present Section 4 describes the particular case

Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity of an infinite 5-phase laminate with interface normals parallel to the direction 𝒆 2 and volume fractions 𝑓 1 = 0.25, 𝑓 2 = 0.1, 𝑓 3 = 0.3, 𝑓 4 = 0.2 and 𝑓 5 = 0.15. Cubic elastic constants of Cu at room temperature were considered with 𝑐 11 = 170 GPa, 𝑐 12 = 124 GPa and 𝑐 44 = 75 GPa for the five phases. Crystallographic orientations of the phases were randomly selected. For each phase, a slip of magnitude 0.02 was assigned to the most stressed system within a pure elastic loading, i.e. among the 12 FCC slip systems {111}⟨110⟩, the system with the highest absolute value of its resolved shear stress |𝜏 | was selected and affected 𝛾 = ±0.02 so that 𝜏𝛾 > 0. No thermal strain was considered. The FE simulation was performed thanks to the use of the open source software suite freefem++ [START_REF] Hecht | New development in freefem++[END_REF]. A structured mesh made of 270 000 triangular quadratic Lagrangian elements was used for the meshing of a bar of length 100 (see Figure 3), which corresponds to 4 nodes along 𝒆 1 and 𝒆 3 and 5001 nodes along 𝒆 2 . Periodic boundary conditions were imposed along the lateral sides. Arbitrary and fixed velocities were imposed on the bottom and upper faces of the bar (𝑢 1 = -0.05, 𝑢 2 = 0, 𝑢 3 = -0.1 on the bottom face, 𝑢 1 = 0.05, 𝑢 2 = 0.1, 𝑢 3 = 0.1 on the top face). Phase 1 was affected to the region 75 < 𝑥 2 < 100, phase 2 to 65 < 𝑥 2 < 75, phase 3 to 35 < 𝑥 2 < 65, phase 4 to 15 < 𝑥 2 < 35 and phase 5 to 0 < 𝑥 2 < 15. A first simulation was performed in pure elasticity in order to determine the most stressed systems and then a second simulation was performed considering static plastic strains in each phase.

The distribution of stresses obtained by the FE simulation is shown in Figure 3. For comparison, stresses were then computed from Equations ( 15), ( 23) and ( 25) considering the same elastic constants, the same crystallographic orientations of the phases, the same plastic strains and taking Σ 𝑖 as the average stress vector given by the FE simulation. In Figure 4, it can be seen that the matching between the stresses obtained numerically by a FE simulation and the ones computed analytically from Equations ( 15), ( 23) and ( 25) is excellent. The tiny differences are due 15), ( 23) and ( 25) (circles).

Thiebaud Richeton

Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity to the fact that the FE method is an approximate numerical method. In particular, it was checked that the stress differences decrease with the number of nodes along 𝒆 2 . These results, along with the several other FE simulations performed, prove the correctness of the analytical formulas established in Section 3.

Application to polycrystals with elongated grains

As an illustration of potential applications, the N-phase laminate model can be used as a quick way to estimate the level of internal stresses in metals and alloys having very elongated grains due to rolling or forging processes. For instance, near-𝛽 titanium alloys like Ti-1023 used in forged pieces of aircraft for landing gears or rotor systems can achieve high specific strength thanks to the building-up of complex 𝛼/𝛽 microstructures. However, these Ti-1023 alloys still contain very large domains of close 𝛽 orientations which corresponds to the prior-𝛽 grains that were deformed during the forging steps above the 𝛽 transus [START_REF] Lhadi | Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of 𝛽-Ti alloys[END_REF][START_REF] Lhadi | Elasto-viscoplastic tensile behavior of as-forged Ti-1023 alloy: Experiments and micromechanical modeling[END_REF]. These prior-𝛽 grains are millimeter size and are highly elongated along a same direction corresponding to the axial axis of the billet. Due to the significant elastic anisotropy of the cubic 𝛽 phase (Purushottam raj purohit et al. 2021), strong mechanical contrasts exist between these millimeter size regions which might be at the origin of the early crack initiation sometimes observed in those materials [START_REF] Lhadi | Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of 𝛽-Ti alloys[END_REF][START_REF] Lhadi | Elasto-viscoplastic tensile behavior of as-forged Ti-1023 alloy: Experiments and micromechanical modeling[END_REF]. In order to get a very first insight into such a mechanical issue, it can thus be interesting to simplify the microstructure as an infinite laminate composed of different 𝛽 orientations. By doing so, the focus is only on the mechanical interactions between prior-𝛽 grains and on a same kind of interfaces, i.e. the predominant long ones which are parallel to the axial axis of the billet. Hence, a laminate composed of 100 000 different 𝛽 orientations of equal volume fraction was considered as an application of the methodology described in Section 3. The orientations were generated from the open-source software package Neper [START_REF] Quey | Nearly uniform sampling of crystal orientations. Thiebaud Richeton Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity[END_REF]) so that their distribution was nearly uniform in presence of cubic symmetry. It is noteworthy that the effect of a crystallographic texture could have considered as well without difficulty. Elastic constants of Ti-1023, 𝑐 11 = 92.6 GPa, 𝑐 12 = 82.5 GPa and 𝑐 44 = 43.5 GPa, were taken from Purushottam raj purohit et al. ( 2021) who provided an estimation based on a method coupling Bayesian inference analysis, high energy X-ray diffraction and elastic self-consistent modeling. These constant values correspond to a strong elastic anisotropy ratio of 𝐴 = 8.6. The laminate interfaces normal was set parallel to the direction 𝒆 2 of the Cartesian frame (𝒆 1 , 𝒆 2 , 𝒆 3 ). A uniaxial macroscopic stress of magnitude 300 MPa was applied without consideration of any plastic or thermal strains. The direction of the uniaxial stress was varied by rotation around 𝒆 3 . This rotation was described by an angle 𝜃 , 𝜃 = 0°corresponding to an uniaxial stress along 𝒆 1 and 𝜃 = 90°along 𝒆 2 . The stresses in each grain were computed from Equations ( 15), ( 23) and ( 25) and then the maximum von Mises stress among the 100 000 grains was recorded and plotted with respect to 𝜃 , see Figure 5. For comparison, a 2-phase laminate made of grains of equal volume fraction was considered as well. The cristallographic orientations of the 2 grains were set so that the directions [111] and [ 110] in one grain and [100] and [010] in the other grain were parallel to 𝒆 1 and 𝒆 2 , respectively. ⟨111⟩ and ⟨100⟩ are actually the directions of maximum and minimum directional Young's modulus, respectively. Hence, the mechanical contrast should very strong for this 2-grain configuration.

In Figure 5, it is interesting to observe that the maximum von Mises stress of the 2-phase laminate is always below the one with 100 000 orientations despite the strong elastic contrast considered. However, it was shown in Section 3 that the stresses in a given phase of a N-phase laminate can also be deduced from a 2-phase laminate where the given phase is in interaction with the whole laminate (see Figure 2 and Equation ( 25)). Accordingly, the elastic contrast between the phase where the maximum von Mises stress is located and the effective laminate of 100 000 grains should normally not exceed the elastic contrast that exists between [111] ∥ 𝒆 1 and [100] ∥ 𝒆 1 -oriented grains. The difference of maximum von Mises stress is actually rather due to a volume fraction effect, incompatibility stresses in a phase of volume fraction 𝑓 scaling indeed as 1 -𝑓 , see Equation ( 6). Yet, the grain volume fraction is 𝑓 = 0.5 in the 2-phase laminate whereas it is negligible (𝑓 = 10 -5 ) in the laminate with 100 000 orientations. Similar results could also be obtained by considering strong mechanical contrasts resulting from difference of plastic or thermal strains. As a consequence, it can be inferred that considering a 2-phase laminate of equal volume fraction with extremal mechanical contrast may lead to underestimate the maximal level of internal stresses in polycrystals. Furthermore, the difference of positions of the local extrema between the 2-phase laminate and the one with 100 000 orientations can also be noticed in Figure 5.

Comparisons with elastic self-consistent models and discussion

It is noteworthy that the internal stress evaluation performed in Section 5 could also have been made from an elastic self-consistent model [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF][START_REF] Kröner | Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls[END_REF]Hill 1965a) with grains more realistic modeled as ellipsoidal inclusions. In a self-consistent model, the effective compliance tensor can be computed as

𝑆 𝑖 𝑗 = ⟨𝑠 𝑖𝑘 𝐵 𝑘ℓ ⟩⟨𝐵 ℓ 𝑗 ⟩ -1 , (27)
where the stress concentration tensor 𝐵 𝑚 𝑖 𝑗 for any phase 𝑚 is given by

𝐵 𝑚 𝑖 𝑗 = (𝑠 𝑚 𝑖𝑘 + M𝑖𝑘 ) -1 (𝑆 𝑘 𝑗 + M𝑘 𝑗 ), (28)
and where M𝑖 𝑗 is the constraint tensor defined by

M𝑖 𝑗 = (𝛿 𝑖𝑘 -𝑆 𝐸 𝑖𝑘 ) -1 𝑆 𝐸 𝑘ℓ 𝑆 ℓ 𝑗 . ( 29
)
𝛿 𝑖 𝑗 is the Kronecker delta and 𝑆 𝐸 𝑖 𝑗 is the interior Eshelby tensor [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF], here expressed as a 6 × 6 matrix consistent with the used Voigt notation (Barnett and Cai 2018). The stress tensor in a given phase is then directly deduced from the stress concentration tensor

𝜎 𝑚 𝑖 = 𝐵 𝑚 𝑖 𝑗 Σ 𝑗 . (30)
For ease of direct comparisons with the results of the laminate model, this set of equations was also implemented in the Matlab software. The results of Figure 5 were then compared with the maximum von Mises stresses obtained from elastic self-consistent models considering oblate spheroidal grains (𝑎 = 𝑏 > 𝑐) of different aspect ratios (Figure 6), the small axis of the oblate spheroids being parallel to the interfaces normal in the laminate model. A convergence towards the laminate solution is observed by increasing the aspect ratio 𝑎/𝑐 of the oblate spheroids, both for the case with 100 000 orientations and for the specific case with two orientations only. For the case with 100 000 orientations, it is seen that, as long as elongated grains are considered, the laminate model provides fairly good qualitative trends for the variation of the maximal von Mises stress with a rotation of the uniaxial macroscopic stress. However, the maximal von Mises stresses are always overestimated compared to the values provided by self-consistent models with more realistic grain aspect ratios. These over estimations represent about 20 % of the uniaxial macroscopic stress for 𝑎/𝑐 = 10 and about 30 % for 𝑎/𝑐 = 5.

Besides, it must be underlined that Equation ( 27), from which is computed the effective compliance tensor 𝑆 𝑖 𝑗 in the self-consistent model, is an implicit equation since both the stress concentration tensor and the Eshelby tensor depend on 𝑆 𝑖 𝑗 , see Equations ( 28) and ( 29). In the present work, this implicit equation was solved through a fixed-point iteration method [START_REF] Walpole | On the overall elastic moduli of composite materials[END_REF]. This might however cause convergence issues for strong mechanical contrasts. By opposition, the laminate methodology of Section 3 is fully explicit and thus, is expected to be more robust and faster. Moreover, the Eshelby tensor 𝑆 𝐸 𝑖 𝑗 in an anisotropic medium cannot be expressed analytically and must be computed numerically. For this work, the Eshelby tensor was computed thanks to the Matlab code provided by Cai 2018 where 𝑆 𝐸 𝑖 𝑗 is obtained from the Hill P tensor [START_REF] Hill | Continuum micro-mechanics of elastoplastic polycrystals[END_REF]) which is itself computed from numerical integrals that are performed using the quadv function in Matlab. On the contrary, the expression of the tensor 𝐺 𝑖 𝑗 in Equation ( 7) is fully analytical and hence very simple to code. Actually, it can be shown that the Eshelby stress solutions can be directly obtained from the tensor 𝐺 𝑖 𝑗 by appropriate integrations [START_REF] Richeton | From bicrystals to spherical inclusions: A superposition method to derive analytical expressions of stress fields in presence of plastic strain gradients[END_REF].

In order to illustrate the numerical differences between the two types of model, the effect of the number of phases on the computation of the effective stiffness tensor, as well as the effect of the absolute tolerance value set in the quadv function in Matlab to evaluate the Eshelby tensor was analyzed. First, Figure 7 shows the relative difference of the norm of the effective stiffness tensor of 100 000 uniform orientations computed by the laminate model and an elastic self-consistent model considering oblate spheroids with 𝑎 = 𝑏 = 100𝑐 as a function of the absolute tolerance set in the quadv function. It is known that, the more the grains are elongated, the more accuracy is needed on the numerical evaluation of the Eshelby tensor [START_REF] Gavazzi | On the numerical evaluation of Eshelby's tensor and its application to elastoplastic fibrous composites[END_REF]. For the present extreme case with 𝑎 = 𝑏 = 100𝑐, Figure 7 shows that an absolute tolerance lower than 10 -12 should be set in order to obtain a convergence towards an accurate result. At the same time, Figure 7 exhibits also the average ratio between the CPU time needed by the elastic self-consistent to achieve the calculation over the one needed by the laminate model. This ratio is about 7 showing, as expected, that the laminate is much faster. Moreover, this ratio increases exponentially for absolute tolerance smaller than 10 -13 . This means that, when using a self-consistent model, a preliminary numerical study is needed to fix the value of the absolute tolerance (or the number of Gauss points if another numerical method is used as in [START_REF] Gavazzi | On the numerical evaluation of Eshelby's tensor and its application to elastoplastic fibrous composites[END_REF])) in order to both ensure a sufficient accuracy on the evaluation of the Eshelby tensor and also to maintain a reasonable CPU time. Then, Figure 8 shows that this relative difference of the norm of the effective stiffness tensor decreases with the number of phases considered. However, it must be underlined that the relative difference is already pretty small (i.e. smaller than 4 × 10 -5 ) when considering only two phases. This small difference agrees with the convergence towards the laminate solution displayed in Figure 6 when the aspect ratio was increased in the case considering two orientations only. This result can actually be understood from Equation ( 24). This equation shows indeed that the stresses within a phase of a two-phase laminate are the same as those in a phase of negligible volume fraction and same properties belonging to a two-phase laminate where the other phase has the effective properties of the first two-phase laminate. In the same way, the self-consistent model computes the stresses as a result of the mechanical interactions between one grain and the infinite homogeneous equivalent medium made of the two grains.

Conclusion

This work provides a fast explicit method to compute the effective behavior and the distribution of stresses in a general N-phase laminate made of parallel, planar and perfectly bonded interfaces.

The formulas are based on the solutions of the two-phase laminate. They are very compact which is convenient for numerical implementation and direct use in engineering applications. In contrast with previous approaches, this contribution considers, at the same time, an arbitrary number of phases, a completely general anisotropic elasticity as well as different thermal and plastic strains in the phases. In addition, a simple direct analytical formula is derived in order to compute the stress in a given phase once the effective behavior of the laminate has been determined. The developed method proves that both the effective properties of the laminate and the stresses in the phases are insensitive to the stacking order of the phases. The code to compute the effective behavior and the stress distribution in an arbitrary N-phase laminate is made freely available for download (link in the Supplementary materials section).

Besides, the correctness of the proposed method is checked by comparisons with finite element simulation results on a same boundary value problem, showing excellent agreements. An application of the method is also performed for a near-𝛽 titanium alloy with elongated grains, by comparing, for different elastic loadings, the level of internal stresses within a N-phase laminate made of 100 000 orientations and a 2-phase laminate with maximal elastic contrast. Interestingly, the maximum von Mises stress of the 2-phase laminate is always the lowest despite the strong elastic contrast considered, which is explained by a volume fraction effect. Similar results could be obtained by considering strong mechanical contrasts resulting from difference of plastic or thermal strains.

Finally, comparisons with elastic self-consistent models considering oblate spheroidal grains of different aspect ratios are performed. A convergence towards the laminate solution is observed by increasing the aspect ratio of the oblate spheroids. However, it is noteworthy that the laminate model is based only on fully analytical and explicit equations, which is more robust from a numerical point of view. It is also shown that the execution of the laminate model is much faster that that of self-consistent models.

As a perspective, it may be worth underlying that the developed method might be extended to obtain the effective strain energy potential of a N-phase viscoplastic laminate, i.e. grouping the N phases in pairs and applying the the existing compact expressions for the effective strain energy potential of a two-phase laminate [START_REF] Debotton | Transversely isotropic sequentially laminated composites in finite elasticity[END_REF][START_REF] Idiart | Modeling the macroscopic behavior of two-phase nonlinear composites by infinite-rank laminates[END_REF]) for each pair iteratively until one phase with the effective strain energy potential is obtained. The strains in each phase could then be deduced from a top-down algorithm.

Figure 2

 2 Figure 2 Description of the two configurations used to establish Equation (24).

Figure 3

 3 Figure 3 Distribution of 𝜎 1 (left), 𝜎 3 (middle) and 𝜎 5 (right) in MPa obtained from FE simulation.

Figure 4

 4 Figure4Comparisons of stresses: FE simulation (lines) and computed from Equations (15), (23) and (25) (circles).

Figure 5

 5 Figure 5 Comparisons of the maximum von Mises stresses obtained between a laminate made of 100 000 uniform orientations (dashed line) and a 2-phase laminate made of [111] ∥ 𝒆 1 and [100] ∥ 𝒆 1 -oriented grains (solid line). 𝜃 represents a rotation around 𝒆 3 of the uniaxial stress direction. Cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).

Figure 6

 6 Figure 6 Comparisons of the maximum von Mises stress obtained with the laminate model and elastic self-consistent models (SC) considering oblate spheroidal grains (𝑎 = 𝑏 > 𝑐) of different aspect ratios (𝑎/𝑐). Top: case of the two orientations considered in Figure 5. Bottom: case of 100 000 uniform orientations. 𝜃 represents a rotation around 𝒆 3 of the uniaxial stress direction, see Section 5. Cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).

Figure 7

 7 Figure 7 Black '+' (left axis): relative difference of the norm of the effective stiffness tensor of 100 000 uniform orientations computed by the laminate model and an elastic self-consistent model considering oblate spheroidal grains with 𝑎 = 𝑏 = 100𝑐 as a function of the absolute tolerance set in the quadv function in Matlab to compute the Eshelby tensor. Red '×' (right axis): corresponding ratios between the CPU time needed by the elastic self-consistent to achieve the calculation over the one needed by the laminate model. Simulations were repeated 10 times to provide average CPU time values. The Frobenius norm of a matrix was considered. The cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).

Figure 8

 8 Figure 8 Relative difference of the norm of the effective stiffness tensor computed by the laminate model and an elastic self-consistent model considering oblate spheroidal grains with 𝑎 = 𝑏 = 100𝑐 as a function of the number of phases considered. Phases correspond to crystallographic orientations randomly picked within a list of 100 000 uniform orientations. Simulations were repeated 1000 times to provide average values. The Frobenius norm of a matrix was considered. The absolute tolerance in the quadv function in Matlab was set to 10 -12 to evaluate the Eshelby tensor in the self-consistent model. The cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021)

Journal of Theoretical, Computational and Applied Mechanics March 2022 jtcam.episciences.org

Acknowledgements This work was supported by the French State through the program łInvestment in the futurež operated by the National Research Agency (ANR) and referenced by ANR-11-LABX-0008-01 (LabEx DAMAS).

Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.

Authors' contributions The author confirms sole responsibility for carrying the study, performing the numerical simulations and writing the manuscript.

Supplementary Material MATLAB codes to compute the effective behavior and the stress fields in an N-phase laminate, application to an alloy with elongated grains and comparisons with self-consistent models. Freefem++ scripts to check the finite element validity of the laminate solutions. Permalink: swh:1:dir:d7373f7f498776bffea0c0f0507a3a637cf39e1e; origin=https://github.com/AniPlas/Laminate;visit=swh:1:snp:07d3ff0677884a21549531f80eb9b316cefc3a97;anchor= swh:1:rev:6bfcf8de1452757d7b04ee17e24a5d3ef5667f97.

Competing interests

The author declares that they have no competing interests.

Journal's Note JTCAM remains neutral with regard to the content of the publication and institutional affiliations.


Wisama Khalil
email: wisama.khalil@lan.ec-nantes.fr

Philippe Lemoine

A Software package for the Calibration of Robots

This paper presents a software package for the simulation and the realization of the calibration of the geometric parameters of robots. This package which is called GECARO, GEometric CAlibration of RObots, contains new methods which can carry out the calibration of the geometric parameters and the gain transmission ratios of serial robots without the need of external sensors. It contains also classical methods which are based on measuring the position and the orientation of the terminal link using external sensors.

I. INTRODUCTION

The absolute accuracy of a robot depends to a large extent on the calibration of its geometric parameters. Geometric calibration of robots is the process by which the parameters defining the base parameters, link parameters and tool parameters are precisely identified. Classically, the geometric calibration is carried out by solving a system of linear or non-linear equations which is a function of the end-effector pose measurements and of the joint positions [1,...,4]. Recently, some papers have proposed autonomous calibration methods which do not need an external sensor ; -in [START_REF] Bennet | Self-calibration of singleloop, closed kinematic chains formed by dual or redundant manipulators[END_REF][START_REF] Edwards Ch | A Single-Point Calibration Technique for a Six Degree-Of-Freedom Articulated Arm[END_REF][START_REF] Zhong | A new method for autonomous robot calibration[END_REF][START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF] the observation matrix and the error vector have been constructed using the fact that different configurations of the robot can give the same pose of the terminal effector .

-in [START_REF] Nakamura | Robot autonomous error calibration method for off line programming system[END_REF] the observation matrix is based on calculating the vector product of a set of vectors in the same plane, which must be equal to zero. Second order terms appearing in the vector products are neglected, thus the obtained observation matrix contains more approximation than the classical one, and the problem of the identifiable parameters is not considered in this paper.

In this paper we present a software package which can simulate and realize the calibration of serial robots. The following methods are available:

-Position measure: which is based on using an external sensor giving the position of the terminal point. The corresponding motors positions are also needed.

-Frame measure : which is based on using an external sensor giving the position and the orientation of the terminal link. The corresponding motors positions are also needed.

-Position link : which is an autonomous method needing no external sensor [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF]. It can be applied for robots which can achieve the same position of the terminal point by multiple configurations. The data used in the calibration correspond to the motor variables of some sets of robot configurations, where the configurations of each set correspond to the same Cartesian position of the terminal point of the tool. No external sensor is needed.

-Frame link : which is an autonomous method needing no external sensor [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF]. This method can be used for robots which can achieve the same pose (position and orientation) of the terminal link by multiple configurations, thus having at least 6 d.o.f. The data used in the calibration correspond to the motor variables of some sets of robot configurations, where the configurations of each set correspond to the same Cartesian pose.

-Plane link: This is an autonomous method. It uses in the calibration a set of points in a plane. Two methods are available [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF] : in the first, the plane coefficients are supposed known, while in the second, the plane coefficients are unknown, they will be identified as well as the geometric parameters.

In all the methods, the geometric parameter are identified using a linearized model which is solved iteratively using least squares criterion and by updating the identified parameters and the observation matrix after each iteration.

The software package GECARO is running on PC computers and developed under MATLAB, using numerical methods. No symbolic model is needed.

The paper is organized as follows : The parameters defining the robot will be presented in section 2. Sections 3 , 4 and 5 describe the identification methods. Section 6 presents the organization of GECARO. Section 7 contains the conclusion.

2.DESCRIPTION OF THE ROBOT

Two types of parameters are used:

The geometric parameters defining the different frames and the geometric of the robot and general parameters concerning the joints limits, offset, gain rations, and the coupling between the motor and joints variables.

Description of the geometric parameters

We consider serial robots consisting of n joints and n+1 links. Link 0 is the base and link n is the terminal link, frame j is defined fixed on link j. We denote : frame -1 : the fixed reference frame, frame n+1: the tool frame.

The end-effector pose can be calculated with respect to the reference frame by the direct geometric model :

-1 T n+1 = -1 T 0 0 T n (q) n T n+1 (1)
where i T j is the 4x4 transformation matrix defining frame j with respect to frame i. The definition of the link frames will be carried out by Khalil and Kleinfinger notation [START_REF] Khalil | Autonomous calibration of robots using planar points[END_REF]. Frame j is defined such that z j is along the axis of joint j, and x j is perpendicular to z j and z j+1 . Frame j is defined with respect to frame j-1 by the matrix j-1 T j , function of the four parameters ( j , d j ,  j , r j ) .

It have been seen that, the calculation of -1 T n+1 can be obtained as a function of the geometric parameters of n+2 frames represented by four parameters for each frame ( j , d j ,  j , r j ) , except the first frame which is represented by only two parameters different than zero [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF]. The joint variable q j is equal to  j if joint j is rotational and r j if joint j is prismatic.

In the calibration process, we have to identify the deviation of the real parameters from the nominal values, thus to identify : ∆ j , ∆d j , ∆ j , ∆r j (for j = 0, n+1).

If the axis of joint j is parallel to the axis of joint j-1, an additional parameter  j must be considered, the nominal value of  j is equal to zero [START_REF] Hayati | Robot arm geometric link calibration[END_REF].

Using the same procedure, one can identify also the error in the transmission gain of the joints ∆K j (for j = 1, n j ) [START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF].

In GECARO the geometric parameters of the base frame, the robot links, and the tool frame will be represented by a unique table where each row, for j=1,…,nf, represents a frame as given in table 1. Beside the parameters  j ,  j , d j ,  j , r j , we see the following parameters :

.  j : which defines the type of the frame,  j = 0 for j rotational,  j = 1 for j prismatic, and  j = 2 if frame j is fixed.

. K j : the transmission gain ratio for joint j (not frame j), for j=1,...,nj, with nj is the number of joints.

The number of frames is denoted nf. The difference (nfnj) gives the number of rows representing fixed frames.

Definition of general parameters

The following parameters are also needed in the calibration of the geometric parameters : a-Offset: a vector of nf components, containing the nominal values of the Offset of joint variables, for a fixed frame the Offset is equal to zero. b-Qmax and Qmin: the vectors of joint limits, for a fixed frame the corresponding elements are equal to zero. c-Com: the (njxnj) coupling matrix of motor variables, d-Coa: the (nfxnf) coupling matrix of joint variables, e-K: the vector of gain transmission,

The joint variables are calculated as a function of motor variables using the following relation :

q = Coa*G(K)*Com*qm + Offset (2)
where :

q is the (nfx1) joint variables vector, for a fixed frame the joint variable is equal to zero, qm is the (njx1) vector of motor variables ,

G(K) is a (nfxnj) matrix, equal to :       diag (K) 0 (nf-nj)xnj
where 0 ixj is the zero matrix of dimension ixj .

f-The (nfx6) priority matrix. This matrix defines the order of priority of the geometric parameters. Columns 1,...,6 correspond to the parameters , d,, r, , K respectively. The elements of row j give the priority of the parameters of frame j, for instance the element [START_REF] Whitney | Industrial robot forward calibration method and result[END_REF][START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF] defines the priority of the parameter r 3 . During the calculation of the identifiable parameters or the execution of the identification algorithm, the parameters will be arranged as function of the corresponding priority number, the first parameters will be those having the highest priority number. If the elements of the priority matrix are all equal, for instance equal to 1, then the parameters will be arranged in the following order :

(1), ...,(nf), d(1),..., d(nf), (1),…, (nf), r(1),…, r(nf),  (1)
,…, (nf), K(1),…K(nj). The importance of this matrix comes from the fact that we construct the base of identifiable parameters as those corresponding to the first independent columns. j

 j  j d j  j r j  j Kj 1  1  1 d 1  1 r 1  1 K1 … … … … … … … … nj  nj  nj d nj  nj r nj  nj Knj nf  nf  nf d nf  nf r nf  nf 0 Table 1.
Definition of the robot parameters

CLASSICAL CALIBRATION METHODS

The end-effector pose can be calculated with respect to the reference frame by the direct geometric model (DGM) as given in equation ( 1).

Using a first order Taylor development of the parameters errors, a linear differential model defining the deviation of the end effector pose due to the differential error in the geometric parameters can be obtained as [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF] :

y = J (q) X (3) with X = X r -X
y represents the (6x1) vector of the position and orientation error. J is a (6xNp) matrix, Np = number of the parameters to be calibrated. X defines the (Npx1) vector of the errors of the geometric parameters. X r is the vector of the real unknown values of the geometric parameters, X is the vector of the nominal values of the geometric parameters,

The expressions of the columns of J are calculated as given in appendix A [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF].

To identify X, equation (3) will be applied for a sufficient number of configurations q 1 ,…,q m , the corresponding poses will be measured and the y i will be calculated. The resulting linear system of equations will be represented by :

Y = W X (4)
Equation ( 4) will be solved iteratively to get the least squares errors solution. After each iteration the geometric parameters will be updated.

If the orientation of the terminal link is not measurable, only the equations corresponding to the position error will be taken into account (the first three equations of Eq.( 3)).

It can be seen that if some columns of J are dependent in all configurations. Relation (3) can be reduced to [START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF]:

y = J b (q) X b (5)
where J b (q) contains the b independent columns of J, the corresponding parameters are known as identifiable parameters or base parameters they will be denoted by the vector X b.

The determination of the identifiable parameters X b must be done before the identification process. They can be obtained numerically using the QR decomposition of a matrix W similar to that defined in equation ( 4) but obtained using random configurations [START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF].

POSITION LINK AND FRAME LINK CALIBRATION METHODS

The main problem in the classical calibration methods, is the need to have an accurate, fast and not expensive equipment to measure the real end-effector pose. Many sensors have been proposed in the literature but neither fulfill these three conditions. The position (frame) link method can be applied for robots which can achieve the same position (frame) by more than one configuration [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF]. Thus, if q a and q b represent two different configurations giving the same pose of the tool then:

-1 T n+1 (q a , X r ) = -1 T n+1 (q b , X r ) [START_REF] Edwards Ch | A Single-Point Calibration Technique for a Six Degree-Of-Freedom Articulated Arm[END_REF] If the values of the geometric parameters in the model X are different from the real values X r , then using a first order development and by considering only the position equations, to facilitate the presentation, we obtain :

[J t (q b )-J t (q a )]∆X = [ -1 P n+1 (q a ,X)--1 P n+1 (q b ,X)] [START_REF] Zhong | A new method for autonomous robot calibration[END_REF] where J t is the position Jacobian, represented by the first three rows of J.

Eq.( 7) relates the deviation in the geometric parameters to the position error between q a and q b . To estimate the vector ∆X we have to repeat this process for sufficient number of couple of configurations.

If multiple configurations can be obtained for a given position and orientation of the tool, the system (7) will be replaced by :

[J(q b )-J(q a )]∆X =       [ -1 P n+1 (q a ,X) --1 P n+1 (q b , X)]  a,b (8)
where  a,b is the differential rotation between the orientation of the tool in the configurations q a and q b .

The calculation of the identifiable parameters can be carried out as given in Appendix B, by studying the QR decomposition of a matrix W calculated from ( 7) or ( 8) using sufficient random couple of configurations q a (j) and q b (j).

THE PLANAR CALIBRATION [10]

The methods, presented in this section, will be carried out using a set of configurations of the robot, where the terminal point of the robot is in the same plane.

The general equation of a plane is supposed as :

a x + b y + c z + 1 = 0 (9)
where a, b, c represent the plane coefficients,

As the terminal point of the robot is in the plane, then :

a Px(q) + b Py(q) + c Pz(q) + 1 = 0 [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF] where Px, Py, Pz represent the coordinates of the position of the terminal point in the world frame.

Two methods are developed.

The first method

Assuming the coefficients of the plane are known. Using a first order development for X in equation [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF], we obtain :

[a Jx(q) + b Jy(q) + c Jz(q) ] ∆X = -a Px(q) -b Py(q) -c Pz(q) -1 [START_REF] Khalil | Autonomous calibration of robots using planar points[END_REF] where :

Jx is the first row of the Jacobian matrix defined in (3), Jy is the second row of the Jacobian matrix , Jz is the third row of the Jacobian matrix, Equation [START_REF] Khalil | Autonomous calibration of robots using planar points[END_REF] gives a linear equation in ∆X, for each configuration.

Sufficient number of configurations must be used to obtain a system of equations similar to (4) which will be solved iteratively to identify ∆X, the geometric parameters will be updated after each iteration.

The calculation of the identifiable parameters can be carried out numerically, by studying the QR decomposition of a matrix W calculated from Eq.( 11) using sufficient random points in a given plane.

The second method

Assuming that the coefficients of the plane are unknown, from Eq. ( 10) the following relation can be obtained :

[Px(q) Py(q) Pz(q) aJx(q)+bJy(q)+cJz(q)]

        ∆a ∆b ∆c ∆X = -aPx(q)-bPy(q)-cPz(q)-1 (12)
using sufficient number of configurations we can get a linear system as in the standard form (4).

The coefficients of the plane are initialized by calculating the equation of the nearest plane to the terminal points of the given configurations. The Cartesian coordinates of the terminal point are calculated using the direct geometric model of the robot. The solution of the following equation gives the initial values of a, b, c :

      -1 -1 -1 =         Px 1 Py 1 Pz 1 … … … Px m Py m Pz m       a b c (13)
with Px j , Py j , Pz j are the coordinates of the terminal point calculated by the DGM for the configuration q(j).

The calculation of the identifiable parameters can be carried out numerically, by studying the QR decomposition of a matrix W calculated from (12) using sufficient random points in a given plane.

DESCRIPTION OF GECARO

The following menus are available :

Robot, Position measure, Frame measure, Point link, Frame link, Plane link, Help Robot menu is used to define or select the robot to be treated. The following functions are available as submenus :

New, Open, Save, Save as, Quit, Help

The other functions correspond to a given calibration method. Each of them has the following three principal sub-menus a-Identifiable parameters : Calculation of the identifiable geometric parameters (base parameters).

The identifiable parameters can be obtained from the classical parameters by eliminating those which have no effect on the observation matrix or those whose corresponding elements in the priority matrix are equal to zero, and by regrouping some others (whose corresponding columns in the observation matrix are not independent). These parameters are calculated using QR decomposition of an observation matrix which can be calculated using either random values, or by using a real experimental data file. This function can also be used to test if some experimental data are sufficiently exciting (giving a good condition number for the observation matrix, and giving the same identifiable parameters as the random values).

It is to be noted that the identifiable parameters are not unique. It is more practical to identify, if possible, the parameters which can be updated easily in the control system without changing the symbolic direct and inverse geometric models of the robot. Therefore we use the following rules to define default priorities as input to this function (the user can define arbitrarily his input priority matrix) :

-the highest is equal to 6 correspond to joint offsets and motor transmission gains.

-the distances R, and D which are not equal to zero take priority 5, -the angles which are not equal to ±k*pi/2, with k integer, will take the priority 4, -the parameters of the base location or the tool frames, will get the priority 2, -the rest of parameters will get the priority 1,

The user can put equal to zero the elements of the priority matrix corresponding to the parameters which he does not like to identify, either because their values are known accurately, or because he cannot make use of their identified values.

The output of this function is the regrouping relations of the identifiable parameters and a new priority matrix which is similar to the input one, but the non identifiable parameters priorities are put equal to zero. b-Points generation : this function generates a file which can be used to test and simulate the functions of the corresponding method, the generated file contains two matrices : Qm : of dimension (nptxnj). Each row contains the motor variables of a random robot configuration. X : contains the Cartesian pose of dimension nptX3 if the position only is needed or of dimension 3nptx4 if the position and orientation are needed.

The number of points denoted (npt ) is given by the user, this number must be chosen such that sufficient number of equations can be generated. -the file containing the priority matrix, -the file containing Qm and X in classical methods or motor variables matrix Qm for autonomous methods, sufficient number of configurations must be given. In the case of position or frame link methods the number of configurations of each set giving the same pose will be given in a vector (whose dimension is equal to the number of sets).

CONCLUSION

This paper presents the software package GECARO, GEometric CAlibration of RObots. It contains new methods which can carry out the calibration of the geometric parameters and the gain transmission ratios of serial robots without the need of external sensors. It contains also classical methods which are based on measuring the position and the orientation of the terminal link using external sensors. The parameters defining the robot and the main menus of this package are described in the paper. GECARO is running on PC computers and developed using MATLAB, any general open loop robot can be treated directly.

c- Identification :

 Identification This function gives the identified values and the precision of the obtained solution. Two input files are needed :

APPENDIX A: CALCULATION OF J

Assuming the matrix defining frame j with respect to the fixed frame as:

-1 T j =     s j n j a j P j 0 0 0 1 (A-1)

The calculation of the columns of the matrix J can be done as follows [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Autonomous calibration of robots using planar points[END_REF]:

Where :

x denotes the vector product, L i,n+1 is the (3x1) position vector between the origin of frame i and the origin of frame n+1 equal to P n+1 -P i , 0 (3x1) is the (3x1) zero vector.

jq j is equal to j j if j rotational or jr j if j prismatic.

All the vectors of equations (A-2…A-5) are referred to the measuring fixed frame.

APPENDIX B: NUMERICAL CALCULATION OF THE IDENTIFIABLE PARAMETERS

Let us rewrite the system (4) as:

where :

where c is the number of columns, and b is the number of independent columns of W .

We can write:

Using (B-2) , equation ( 4) can be written as:

with

In the identification process, equation (B-3) will be used instead of (4). The solution will give directly X b which is called the base parameters vector. The matrix B is not needed in the identification process.

Numerically, the study of the base parameters is equivalent to study the space spanned by the columns of (rxc) matrix W with r>> c, calculated from J at m random configurations q 1 ,…,q m satisfying the constraint of the problem (terminal point in a plane, ...)

Using QR decomposition the (rxc) matrix W can be decomposed as [START_REF] Dongarra | LINPACK User's Guide[END_REF][START_REF] Lawson | Solving Least Squares Problems[END_REF]:

Q is a (rxr) orthogonal matrix, R is a (cxc) upper triangular matrix. 0 ixj is the (ixj) matrix of zeros.

If the element | R ii | ≤ , the corresponding parameter in X i is not identifiable. is the numerical zero which can be taken as [START_REF] Dongarra | LINPACK User's Guide[END_REF]:

Where  is the machine precision.

Base Parameters as Function of the Standard Parameters

Let us permute the columns of W such that the first b columns are independent.

[

where:

P is a permutation matrix, W 1 represent the b independent columns of W. W 2 represent the (c-b) dependent columns of W.

A QR decomposition of [W P] gives :

Where R 1 is a bxb regular matrix. Then it comes:


Ph W Khalil
email: khalil@lan.ec-nantes.fr

M Lemoine

Gautier

AUTONOMOUS CALIBRATION OF ROBOTS USING PLANAR POINTS

Keywords: Robot calibration, geometric parameters, identification

This paper presents two methods for the calibration of the geometric parameters of serial robots. They are based on putting the terminal point of the robot in a plane, only the position sensors of the joints are required. The methods can be put in practice quickly and precisely. Some practical issues concerning their practical use are given.

I. INTRODUCTION

Geometric calibration of robots is the process by which the parameters defining the base parameters, link parameters and tool parameters are precisely identified. The absolute accuracy of a robot depends to a large extent on the calibration of its geometric parameters. Classically, the geometric calibration is carried out by solving a system of linear or nonlinear equations which is a function of the end-effector pose measurements and of the joint positions [1,...,4]. In practice all the industrial robots are equiped with joint sensors, thus the real problem is to find a precise external sensor to give the pose values precisely and quickly because a lot of configurations are generally needed. Recently, some papers have proposed autonomous calibration methods which do not need an external sensor :

-in [START_REF]Self-calibration of single-loop, closed kinematic chains formed by dual or redundant manipulators[END_REF][START_REF] Edwards Ch | A Single-Point Calibration Technique for a Six Degree-Of-Freedom Articulated Arm[END_REF][START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Nakamura | Robot autonomous error calibration method for off line programming system[END_REF] the observation matrix and the error vector have been constructed using the fact that different configurations of the robot could give the same pose of the terminal effector .

-in [START_REF] Zhong | A new method for autonomous robot calibration[END_REF] the observation matrix is based on equating to zero the vector product of a set of vectors in the same plane. Second order terms appearing in the vector products are neglected, thus the obtained observation matrix contains more approximation than the classical one, and the problem of the identifiable parameters in this method is not considered.

In this paper we use a set of points in a plane to carry out the calibration of the geometric parameters, the observation matrix is constructed using the equation of the plane, two methods are used. In the first, the plane coefficients are supposed known, while they will be identified in the second. The solution is obtained using a linearized model which is solved iteratively using least squares method.

It is to be noted that the idea of using points in a plane has been also used in [START_REF] Tang | A study of three robot calibration methods based on flat surfaces[END_REF], but we do not classify this method as autonomous one, because it is based on putting the end effector point in a plane which must be the fixed x-y plane, thus only the equation z= 0 is used. In our method the plane may be put any where and the constraint equation of the plane is used and its coefficients can be identified.

II. CLASSICAL CALIBRATION

We consider serial robots consisting of n joints and n+1 links. Link 0 is the base and link n is the terminal link, frame j is defined fixed on link j. We denote : frame -1 : the fixed reference frame, frame n+1: the frame fixed to the tool.

The end-effector pose can be calculated with respect to the reference frame by the direct geometric model (DGM) :

-1 T n+1 = -1 T 0 0 T n (q) n T n+1 (1)
where i T j is the 4x4 transformation matrix defining frame j with respect to frame i.

The definition of the link frames will be carried out by Khalil and Kleinfinger notation [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF].

Frame j is defined with respect to frame j-1 by the matrix j-1 T j , function of the four parameters ( j , d j ,  j , r j ) .

It have been seen that, the calculation of -1 T n+1 can be obtained as a function of the geometric parameters of n+2 frames represented by four parameters for each frame ( j , d j ,  j , r j ) [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF] (0 and d0 are identically equal to zero). The joint variable q j is equal to  j if joint j is rotational and r j if joint j is prismatic.

In the calibration process, we have to identify the deviation of the real parameters from the nominal values, thus to identify : ∆ j , ∆d j , ∆ j , ∆r j (for j = 0, n+1, except 0 and d0).

If the axis of joint j is parallel to the axis of joint j-1, an additional parameter  j must be considered, the nominal value of  j is equal to zero [START_REF] Hayati | Robot arm geometric link calibration[END_REF].

The error in the transmission gain of the joints ∆K j can be also taken into account [START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF].

The identification model

The tool frame pose can be calculated by the direct geometric model given as -1 T n+1 which is a non-linear function of the geometric parameters. Using a first order Taylor development of the parameters errors, a linear differential model will be used in the calibration. If the errors are big, an iterative procedure can be applied.

The differential vectors defining the deviation of the end effector due to the differential error in the geometric parameters can be given by [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF] :

y = J (q) X (2)
where: y represents the (6x1) vector of the position and orientation error. J is a (6xNp) matrix, Np = number of the parameters to be calibrated. -X defines the (Npx1) vector of the errors of the geometric parameters. It is equal to : X = X r -X where : X r is the vector of the real unknown values of the geometric parameters, X is the vector of the nominal values of the geometric parameters,

The expressions of the columns of J are calculated as given in [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF][START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF].

To identify X, equation (2) will be used for a sufficient number of configurations q 1 ,…,q m , the corresponding poses will be measured and the y i will be calculated. The resulting linear system of equations will be represented by :

Y = W X (3)
where:

W =         J(q 1 ) … J(q m ) , Y =         y 1 … y m (4)
Equation ( 3) will be solved to get the least squares errors solution.

Remark : If the orientation of the terminal link is not measurable, only the equations corresponding to the position error will be taken into account.

The identifiable parameters

It can be seen that if some columns of J are dependent in all configurations. Relation (2) can be reduced to [START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF][START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF]:

y = J b (q) X b (5)
where J b (q) contains the b independent columns of J, the corresponding parameters are known as identifiable parameters or base parameters they will be denoted by the vector X b.

The determination of the identifiable parameters X b must be done before the identification process. They can be obtained numerically using the QR decomposition of a matrix W like that given in (3) in which q 1 ,…,q m represent random values [START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF].

The choice of b independent columns to be kept is not unique. We choose the columns corresponding to the parameters which can be updated easily in the control system and do not change the existing inverse and direct geometric models

III. THE NEW CALIBRATION METHODS

The main problem in the classical method, is the need to have an accurate, fast and not expensive equipment to measure the end-effector pose.

The methods, presented in this section, will be carried out using a set of configurations of the robot, where the position of the terminal point of the robot is in the same plane. Practically this set of points can be easily done [START_REF] Zhong | A new method for autonomous robot calibration[END_REF].

The general equation of a plane is supposed as :

a x + b y + c z + 1 = 0 (6)
where a, b, c represent the plane coefficients, Since the terminal point of the robot is in the plane, thus :

a Px(q) + b Py(q) + c Pz(q) + 1 = 0 [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF] where Px, Py, Pz represent the coordinates of the terminal point in the world frame. Equation ( 7) will be exploited to carry out the calibration of the robot parameters. Two methods are developed :

The first supposes that the coefficients of the plane (a, b, c) are known, while the second method supposes that the coefficients of the plane (a, b, c) are unknown.

The first method

Assuming the coefficients of the plane are known, this can be practically useful in the case of interchanging two robots, such that a plane in the working area is already identified by the first robot. The equations of the plane can be also obtained easily using an external sensor where a very limited number of points is needed (theoretically only three points are needed) to determine the plane coefficients. If the values of the geometric parameters in the model X are different than the real values X r , then using a first order development for X in equation ( 7), we obtain : a (Px(q) + Jx(q) ∆X) + b (Py(q) + Jy(q) ∆X) + c(Pz(q) + Jz ∆X)+ 1 = 0 (8)

Relation ( 8) can be written as :

[a Jx(q) + b Jy(q) + c Jz(q) ] ∆X = -a Px(q) -b Py(q) -c Pz(q) -1 [START_REF] Zhong | A new method for autonomous robot calibration[END_REF] where :

Px(q) , Py(q) , Pz(q) represent the coordinates of the terminal point calculated as a function of q using the direct geometric model and the nominal values of the geometric parameters, Jx is the first row of the Jacobian matrix defined in (2), Jy is the second row of the Jacobian matrix defined in (2), Jz is the third row of the Jacobian matrix defined in (2), Equation ( 9) gives a linear equation in ∆X, for each configuration.

Sufficient number of configurations must be used to obtain a system of equations on the standard form ( Y = W X ). This system of equation will be solved iteratively to identify ∆X, the geometric parameters will be updated after each iteration.

The calculation of the identifiable parameters can be carried out numerically using a procedure similar to that given in [START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF], by studying the QR decomposition of a matrix W calculated from Eq.( 9) using m random points in a given plane.

The second method

Assuming that the coefficients of the plane are unknown, from equation ( 7) the following equation can be obtained :

(a+∆a) (Px(q)+Jx(q) ∆X)+(b+∆b) (Py(q)+Jy(q) ∆X)+(c+∆c)(Pz(q)+Jz(q)∆X)+1= 0 [START_REF] Tang | A study of three robot calibration methods based on flat surfaces[END_REF] neglecting the second order terms, we get :

[Px(q) Py(q) Pz(q) a Jx(q)+b Jy(q)+c Jz(q)]

      ∆a ∆b ∆c ∆X = -aPx(q)-b Py(q)-c Pz(q)-1 (11)
we have to repeat this process for sufficient number of configurations to get a linear system as in the standard form .

The coefficients of the plane are initialized by calculating the equation of the nearest plane to the terminal points of the given configurations. The Cartesian coordinates of the terminal point are calculated using the direct geometric model of the robot. The solution of the following equation gives the initial values of a, b, c :

           -1 -1 -1 =         Px 1 Py 1 Pz 1 … … … Px m Py m Pz m         a b c (12)
with Px j , Py j , Pz j are the coordinates of the terminal point calculated by the DGM for the configuration q(j).

The calculation of the identifiable parameters can be carried out numerically, by studying the QR decomposition of a matrix W calculated using m random points in a given plane, where each row represents an equation similar to [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF].

Practical issues

The following remarks are reported as a result of experiments and simulation of different robots using the given planar methods :

-The first method, which supposes that a, b, c are known, is the best. The convergence to the real parameters always takes place even if the initial errors in the geometric parameters are very big.

-The second method works well for moderate initial errors in the parameters such that the error pose is about 0,30 meters, which is far away of the performances of industrial robot.

IV. EXAMPLE

In this section we present the non-identifiable parameters, of a PUMA 560 robot with arbitrarily base location and terminal end effector, using four different methods of calibration. In the following bold letters gives the parameters having no effect on the model while the others are regrouped to some other parameters : i-classical method using position measure the non identifiable parameters are : 7, r3, r7, 7

ii-position link method without external sensor [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF], the non identifiable parameters are : 1, d1, r0, 1, 0, 7, r1, r6, r7, 7, r3, iii-using known plane, the non identifiable parameters are : 7, r0, r6, r7, 7, 0, r1, r3. iv-using unknown plane, the non identifiable parameters are : 7, 1, d1, r0, r6, r7, 1, 7, 0, r1, r3.

It is to be noted that in methods 2, 3 and 4 most of the non-identifiable parameters belong to the parameters whose axes are fixed, such as : r0, ..., 1. The use of an external sensor will permit to identify these parameters in a frame fixed with respect to this sensor, which is in general not the desired world frame. Thus actually, the problem of identifying these parameters exists also in the case of external sensor. To solve this problem, the coordinates of some known poses in the environment must be available in order to put the robot in these poses and use the classical method to identify them.

V. CONCLUSION

This paper presents new methods for the calibration of the geometric parameters of serial robots. They are based on putting the terminal point of the robot in a plane, only the position sensors of the joints are required. The methods can be put in practice quickly and precisely. The identifiable parameters given by the proposed method are slightly different than those of the classical method using an external position sensor. A software package "GECARO" (GEometric CAlibration of RObots) has been developed to carry out the calibration of real robots or to study the simulation of the calibration process using the classical methods and the the autonomous methods presented in this paper and those presented in [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF].