Charles Mareau
email: charles.mareau@ensam.eu

Camille Robert
email: camille.robert@ensam.eu

Different com posite voxel methods for the numerical homogenization of heterogeneous inelastic materials with FFT-based techniques

Keywords:

FFT-based homogenization methods aim at calculating the effective behavior of heterogeneous materials with periodic microstructures. These methods operate on a regular grid of voxels, and hence require an appropriate spatial discretization of periodic microstructures. However, when different microstructural length scales are involved, it is not always possible to have sufficient spatial resolutions to explicitly consider the influence of fine microstructural features (e.g. voids, second-phase particles). To circumvent this difficulty, one solution consists of using composite voxel methods to define the effective properties and the effective internal variables of heterogeneous voxels.

In this work, different composite voxel methods are proposed to deal with inelastic materials with multiple length scales. These methods use simple homogenization rules to calculate the effective behavior of heterogeneous voxels. The first part of this paper is dedicated to the description of the composite voxel methods, which are based either on the Voigt, laminate structure or Mori-Tanaka approximations. In the second part, these methods are used to model the elasto-plastic behavior of a pearlitic steel polycrystalline aggregate. According to the results, the Voigt approximation, which ignores morphological features, is not appropriate for treating heterogeneous voxels. When morphological information is accounted for, with either the laminate structure or Mori-Tanaka approximations, a better agreement with experimental observations is obtained. Though none of these methods is universal, they offer some possibilities to investigate the mechanical behavior of heterogeneous materials involving multiple length scales.

Introduction

FFT-based homogenization methods, originally proposed by [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] and [START_REF] Müller | Mathematical versus experimental stress analysis of inhomogeneities in solids[END_REF] , provide a convenient way of estimating the effective properties of heterogeneous materials with periodic microstructures. For instance, in a mechanical context, FFT-based methods have been used to investigate the behavior of different types of heterogeneous solids, including hyperelastic [START_REF] Lahellec | Analysis of Inhomogeneous Materials at Large Strains Using Fast Fourier Transforms[END_REF][START_REF] Kabel | Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations[END_REF], thermoelastic [START_REF] Vinogradov | An accelerated FFT algorithm for thermoelastic and non-linear composites[END_REF][START_REF] Anglin | Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions[END_REF], plastic [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], viscoplastic [START_REF] Lebensohn | N-site modeling of a 3d viscoplastic polycrystal using fast Fourier transform[END_REF], viscoelastic [START_REF] Figliuzzi | Modelling the microstructure and the viscoelastic behaviour of carbon black filled rubber materials from 3d simulations[END_REF] and elasto-viscoplastic [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF][START_REF] Eisenlohr | A spectral method solution to crystal elasto-viscoplasticity at finite strains[END_REF][START_REF] Mareau | Experimental and numerical study of the evolution of stored and dissipated energies in a medium carbon steel under cyclic loading[END_REF][START_REF] Robert | A comparison between different numerical methods for the modeling of polycrystalline materials with an elasticviscoplastic behavior[END_REF][START_REF] Mareau | Micromechanical modelling of twinning in polycrystalline materials: application to magnesium[END_REF] materials. Some recent applications have also used FFT-based methods to deal with nonlocal plasticity [START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms[END_REF] and dislocation mechanics [START_REF] Berbenni | A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics[END_REF][START_REF] Brenner | Numerical implementation of static field dislocation mechanics theory for periodic media[END_REF].

Independently on constitutive relations, FFT-based homogenization methods aim at finding the solution to the integral equations deriving from the application of compatibility and static equilibrium conditions. For the numerical resolution of the corresponding equations, different iterative procedures have been proposed [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF][START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. As shown in the comparative study of [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF] , the efficiency of these procedures depends on the property contrast. Also, instead of using a fixed-point algorithm, an alternative strategy, which is based on the conjugate gradient method, has been proposed by [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients[END_REF] to improve convergence properties. This idea has been extended to non-linear constitutive relations by [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF] .

Whatever the iterative procedure is, and except for some specific cases, FFT-based methods do not provide exact solutions to the integral equations, but only numerical approximations whose accuracy depend on different factors.

First, the numerical resolution of integral equations requires the introduction of a homogeneous reference medium to which a Green operator is associated. As shown in the work of [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Willot | Fast Fourier Transform Computations and Build-up of Plastic Deformation in 2D, Elastic-perfectly Plastic, Pixelwise Disordered Porous Media[END_REF] , alternative definitions of the Green operator may contribute to improving both the accuracy of the stress and strain fields and the convergence properties.

Second, FFT-based homogenization methods operate on a regular grid of voxels, and hence require an appropriate spatial discretization of periodic microstructures. Evidently, the accuracy of the numerical approximation is largely influenced by the spatial resolution used for discretizing the volume element.

Third, within the volume of a voxel, different constituents may exist. This situation is for instance encountered in interfacial voxels which connect multiple domains with different properties. For the specific case of interfacial voxels, with either a linear [START_REF] Gélébart | Filtering material properties to improve FFT-based methods for numerical homogenization[END_REF][START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF] or a non-linear behavior ( Kabel et al., 2016a;2016b ), the solution quality can be improved by using some simple homogenization rules to define the effective behavior of such voxels.

When multiple length scales are involved in a microstructure, different constituents may also exist within a voxel. Indeed, when fine microstructural details exist (e.g. nano-particles, micro-voids), high spatial resolutions are usually needed to explicitly consider such details. The computational cost is however often prohibitive and some alternative strategies have to be developed. In this work, different com posite voxel methods are proposed to circumvent this difficulty. These methods treat voxels as heterogeneous media and use homogenization techniques to compute the effective behavior.

This paper is organized as follows. In the first part, the equations governing the mechanical behavior of heterogeneous materials with periodic microstructures are given and the strategy for numerically solving these equations using FFT techniques is recalled. Then, different com posite methods for treating heterogeneous voxels are proposed for materials exhibiting an inelastic behavior. Finally, to demonstrate the interest of such methods, the elastoplastic behavior of a pearlitic steel polycrystalline aggregate is investigated. The different homogenization methods are employed to consider the characteristic two-phased fine microstructure of pearlite. The relevance of the proposed composite voxel methods is discussed from numerical results.

Field equations

To represent a heterogeneous material with periodic microstructure, a unit cell with volume V and boundary ∂V is considered. For the specific case of periodic boundary conditions, the displacement vector u ( x , t ), which is associated with a material point occupying a position x at time t , is decomposed as follows:

u ( x , t ) = ˜ u ( x , t ) + ( E ( t ) + ( t ) ) • x (1)
where ˜ u ( x , t ) is the periodic part of the displacement field, E ( t ) and ( t ) are respectively the macroscopic strain and rotation tensors. The effective response of the heterogeneous material is given by the evolution of the macroscopic stress tensor ( t ), which depends on microstructural features, local constitutive relations and prescribed loading conditions. The macroscopic quantities ( t ), E ( t ) and ( t ) are connected to their microscopic counterparts σ( x , t ), ( x , t ) and ω( x , t ) with the classical averaging relations of homogenization theory:

( t ) = 1 V V σ( x , t ) dV (2) E ( t ) = 1 V V ( x , t ) dV (3) ( t ) = 1 V V ω ( x , t ) dV (4)
The problem of finding the local fields σ( x , t ), ( x , t ) and ω( x , t ) consists of solving a particular set of differential equations resulting from compatibility and equilibrium conditions. More specifically, in the absence of volume forces, the respect of static equilibrium conditions imposes: div ( σ) | x ,t = 0

(5) Also, within the framework of infinitesimal strains and rotations, for compatibility conditions to be satisfied, the strain and rotation fields have to derive from the displacement field:

( x , t ) = sym ( grad ( u ) | x ,t ) (6) ω ( x , t ) = skw ( grad ( u ) | x ,t ) (7)
The resolution of the above differential equations requires specifying the constitutive relations associated with each material point. In the present work, we restrict ourselves to constitutive models for which the local strain and rotation tensors can be decomposed into elastic (superscript e ) and inelastic1 (superscript in ) contributions:

( x , t ) = e ( x , t ) + in ( x , t ) (8) ω ( x , t ) = ω e ( x , t ) + ω in ( x , t ) (9)
In the context of linear elasticity, the microscopic stress tensor σ( x , t ) is connected to the elastic strain tensor e ( x , t ) according to:

σ( x , t ) = C ( x , t ) : e ( x , t ) (10)
where C is the elastic stiffness tensor and S (with S = C -1 ) is the elastic compliance tensor.

The introduction of a reference elastic homogeneous medium, with stiffness C 0 , allows the reformulation of the above problem in the form of a set of integral equations:

( x , t ) = E ( t ) - V G sym x -x : τ x , t dV (11) ω ( x , t ) = ( t ) - V G skw x -x : τ x , t dV (12)
where G sym (respectively G skw ) is the symmetric (respectively skew-symmetric) modified Green operator associated with the reference medium and τ is the polarization field whose definition is:

τ( x , t ) = σ( x , t ) -C 0 : ( x , t ) (13)
In the present work, the reference medium is assumed to be isotropic with a bulk modulus k 0 and a shear modulus μ 0 . These properties are chosen according to the suggestion of [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] :

k 0 = k min k max (14) μ 0 = √ μ min μ max ( 15
)
where k min and k max (respectively μ min and μ max ) are the minimum and maximum values of the bulk modulus (respectively shear modulus) in V .

FFT-based resolution method

In most situations, no analytical solution is available for the integral Eqs. ( 11) and ( 12) . To obtain a numerical solution, the volume element is discretized into N voxels with volume v . For a given voxel, the corresponding barycentric position x is:

x = 1 v v x dv (16)
The discrete form of Eqs. ( 11) and ( 12) is then given by:

¯ ( x , t ) = E ( t ) - N G sym x -x : τ x , t (17) ω ( x , t ) = ( t ) - N G skw x -x : τ x , t (18)
where ¯ (

x , t ) , ω ( x , t ) , σ( x , t ) and τ( x , t ) = σ( x , t ) -C 0 : ¯ ( x , t )
represent respectively the strain tensor, the rotation tensor, the stress tensor and the polarization tensor associated with a voxel of position x at time t . A convenient way of evaluating the convolution operations involved in the above equations consists of using the discrete Fourier transform. Indeed, the application of the Fourier transform, which is denoted by a ' ˆ ' symbol, to Eqs. ( 17) and ( 18) leads to:

ˆ ¯ ξ, t = E ( t ) ξ = 0 -ˆ G sym ξ : ˆ τ ξ, t ξ = 0 (19) ˆ ω ξ, t = ( t ) ξ = 0 -ˆ G skw ξ : ˆ τ ξ, t ξ = 0 (20)
where ξ denotes the frequency vector. The integral Eqs. ( 17) and ( 18) are thus solved in the frequency domain, and the application of the inverse Fourier transform leads to the expression of the fields ¯ ( x , t ) and ω ( x , t ) in the original spatial domain. However, because the polarization field τ( x , t ) depends on the strain field ¯ ( x , t ) , an iterative procedure is generally required to solve (19) .

In this work, the procedure proposed by [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] , which is detailed in Appendix A , is used to solve (19) . Also, for the resolution of Eqs. ( 19) and ( 20) , the expressions of the modified Green operators in the frequency domain need to be known. In this work, a centered finite difference approximation is used to construct these operators from compatibility and equilibrium conditions [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

Composite voxel methods

When discretizing a volume element, the application of the above FFT-based resolution method requires to define the properties associated with each voxel. While the common method consists of assigning the properties of the predominant constituent, more sophisticated strategies can be used. Here, to deal with twoscale microstructures, voxels are decomposed into different complementary domains corresponding to different constituents. Each sub-domain, say α (with α = 1 to M ), has homogeneous properties.

For each voxel of volume v , the introduction of a set of indicator functions η α ( x , x , t ) allows the calculation of the volume fractions φ α of the different constituents:

φ α ( x , t ) = 1 v v η α ( x , x , t ) dv (21)
The indicator functions η α ( x , x , t ) take a unity value when the position x is included in the αth domain of the voxel with barycentric position x and a zero value elsewhere. For the αth domain of the voxel with barycentric position x , the average stress, strain and rotation fields at time t are denoted by σ α ( x , t ) , α ( x , t ) and ω α ( x , t ) .

Assuming the characteristic length scales associated with the different constituents are inferior to the size of a voxel, the effective stress, strain and rotation tensors are then defined from the averaging relations:

σ( x , t ) = α ( φ α ( x , t ) σ α ( x , t ) ) (22) ¯ ( x , t ) = α ( φ α ( x , t ) α ( x , t ) ) (23) ω ( x , t ) = α ( φ α ( x , t ) ω α ( x , t ) ) (24)
For each constituent, the stress tensor σ α ( x , t ) is connected to the strain tensor α ( x , t ) according to an inelastic type of constitutive relation:

σ α ( x , t ) = C α ( x , t ) : α ( x , t ) -in,α ( x , t ) (25)
where C α ( x , t ) and in,α ( x , t ) are the elastic stiffness tensor and the inelastic strain tensor for the αth constituent.

In the following, different homogenization methods for assigning properties to heterogeneous voxels are presented. For a given voxel, these methods aim at finding the effective behavior which is given by the effective elastic stiffness tensor C ( x , t ) and the ef- fective inelastic strain tensor ¯ in ( x , t ) . These tensors connect the effective stress tensor to the effective strain tensor according to:

σ( x , t ) = C ( x , t ) : ¯ ( x , t ) -¯ in ( x , t ) (26)
When inelastic strains are treated as eigenstrains, the problem of finding C ( x , t ) and ¯ in ( x , t ) is solved using the general equations for the homogenization of heterogeneous linear thermoelastic media. According to the classical result of homogenization theory of [START_REF] Rosen | Effective thermal expansion coefficients and specific heats of composite materials[END_REF] , and whatever the retained homogenization method is, the effective stiffness tensor C ( x , t ) and the effec- tive inelastic strain tensor ¯ in ( x , t ) are given by:

C ( x , t ) = α φ α ( x , t ) C α ( x , t ) : A α ( x , t ) (27) ¯ in ( x , t ) = α φ α ( x , t ) B α T ( x , t ) : in,α ( x , t ) (28)
where A α and B α are the elastic strain and stress concentration tensors for the αth constituent. These tensors, which need to be specified for each homogenization method, are connected to each other with the following relation:

B α ( x , t ) = C α ( x , t ) : A α ( x , t ) : S ( x , t ) (29)
where S = C -1 is the effective elastic compliance tensor.

In the following, for convenience, the dependence with position and time will be omitted unless needed.

Voigt approximation

For a given voxel, a simple approximation consists of considering the strain and rotation fields to be uniform within v . Thus, localization equations, which connect α and ω α to ¯ and ω , are:

α = ¯ (30) ω α = ω (31)
The elastic strain concentration tensor A α , which is needed to compute the effective elastic stiffness tensor and the effective inelastic strain tensor, is given for each domain by:

A α = I ( 32
)
where I is the fourth rank identity tensor.

While the Voigt approximation is straightforward to implement, it has two major drawbacks. First, the influence of morphological factors is completely ignored, only the volume fractions are considered. Second, it is an upper bound approximation which often leads to an important overestimation of internal stresses.

Laminate structure approximation

Another approximation is obtained by treating a heterogeneous voxel as a laminate structure. For such a voxel, the planar interfaces separating the different constituents have the same normal direction, which is defined by the unitary vector k . Considering the interfaces to be perfect, the continuity of the traction vector imposes:

σ α • k = σ • k (33)
Also, for the strain and rotation tensors to be compatible, the following equalities must hold:

α : ( l l ) = ¯ : ( l l ) (34) α : ( m m ) = ¯ : ( m m ) (35) ( α + ω α ) : ( l m ) = ( ¯ + ω ) : ( l m ) (36) ( α + ω α ) : ( k l ) = ( ¯ + ω ) : ( k l ) (37) ( α + ω α ) : ( k m ) = ( ¯ + ω ) : ( k m ) (38)
In the above equations, l and m are two mutually orthogonal vectors which are orthogonal to k . The complete set of strain and rotation localization equations is obtained by introducing constitutive relations into Eq. ( 33) :

C α : α -in,α • k = C : ¯ -¯ in • k (39)
For a given voxel with prescribed strain and rotation tensors ( ¯ and ω ), relations (34) to (39) allow the determination of the strain and rotation tensors associated with the αth constituent. The application of relation (39) however requires the knowledge of the effective elastic stiffness tensor C and the effective inelastic strain tensor ¯ in . To compute these tensors, relations (34) to (39) are grouped according to:

F α : α -H α : in,α = F : ¯ -H : ¯ in ( 40
)
The elastic strain concentration tensor A α for the αth domain is then deduced from the averaging condition (23) :

A α = F α -1 : α φ α F α -1 -1 (41)
The effective stiffness tensor C and the effective inelastic strain tensor ¯ in are finally obtained from relations ( 27) and ( 28) . With the above definition of the elastic strain concentration tensors, the expression of the effective stiffness tensor C is valid for any class of material symmetry; it is equivalent to the general expression proposed by [START_REF] Milton | The theory of composites[END_REF] for laminate composite materials.

In contrast with the Voigt approximation, the laminate structure approximation considers some kind of morphological information to determine the effective behavior of a given voxel. Though this approach does not require any sophisticated numerical procedure, it is restricted to lamellar structures.

Mori-Tanaka approximation

For some voxels, a matrix-inclusion type of structure can sometimes be assumed. In such situations, if the matrix constituent represents a large portion of the heterogeneous voxel, the Mori-Tanaka approximation [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] can be used to estimate the corresponding effective properties. For the matrix constituent, the superscript α = μ is used while inclusion constituents are denoted by the superscript α = β. According to the Mori-Tanaka approximation, the strain and rotation tensors for inclusion constituents are given by:

β = A β : μ + A β : E β sym : S μ : C β : in,β -C μ : in,μ
(42)

ω β = ω μ + E β skw : E β -1 sym : β -μ ( 43
)
where E β sym and E β skw are the symmetric and skew-symmetric Eshelby tensors associated with the βth inclusion. The tensor A β is calculated from:

A β = I + E β sym : S μ : C β -C μ -1 (44)
The strain and rotation tensors for the matrix constituent ( μ and ω μ ) are obtained from the averaging conditions ( 23) and ( 24) :

μ = A μ : ¯ -A μ : β (φ β A β : E β sym : S μ : (C β : in,β -C μ : in,μ )) (45) ω μ = ω - β φ β E β skw : E β -1 sym : β -μ (46)
with:

A μ = β φ β A β + φ μ I -1 (47)
The following expression of the elastic strain concentration tensor is used to determine the effective stiffness tensor C and the effective inelastic strain tensor ¯ in :

A α = A μ , α = μ A β : A μ , α = β (48)
In some specific situations, the Mori-Tanaka approximation provides an accurate estimation of the effective properties associated with a heterogeneous voxel. Nevertheless, the Eshelby tensor E α sym has to be known to compute the effective properties. While some analytical expressions are available for ellipsoidal inclusions when the elastic properties of the matrix material are isotropic [START_REF] Mura | Micromechanics of defects in solids[END_REF], a numerical procedure is needed to evaluate these tensors in the general case.

Application to pearlitic microstructures

In this section, the mechanical behavior of a fully pearlitic eutectoid steel is investigated. Pearlite is a lamellar aggregate with alternating layers of cementite and ferrite for which the characteristic length scale is typically about a few hundreds of nanometers [START_REF] Marder | The effect of morphology on the strength of pearlite[END_REF][START_REF] Dollar | Influence of deformation substructure on flow and fracture of fully pearlitic steel[END_REF][START_REF] Yahyaoui | Effect of interlamellar spacing on the elastoplastic behavior of c70 pearlitic steel: experimental results and self-consistent modeling[END_REF], much lower than the typical grain size of pearlitic grains. In equilibrium conditions, the volume fractions of ferrite and cementite are respectively about 88% and 12%. Since the explicit discretization of pearlitic microstructures is often not possible, composite voxel methods provide a solution for dealing with such fine microstructures.

Constitutive model

For a given constituent (i.e. ferrite or cementite), a ratedependent elasto-plastic behavior is assumed. Adopting the general framework of crystal plasticity for rate-independent processes [START_REF] Hill | Generalized constitutive relations for incremental deformation of metal crystals by multislip[END_REF], a set of S slip systems is introduced to express the inelastic contribution. Each slip system s (with s = 1 to S ) is defined with a slip direction b s, α and a slip plane normal n s, α . The inelastic strain and rotation tensors ( p, α and ω p, α ) are expressed as a function of the plastic shear strains γ s, α according to:

in,α = s sym b s,α n s,α γ s,α (49) ω in,α = s skw b s,α n s,α γ s,α (50)
The shear stress τ s, α acting on a given slip system is obtained from the projection of the stress tensor σ α :

τ s,α = b s,α • σ α • n s,α (51)
For each slip system, plastic flow is possible only when the critical shear stress g s, α is reached. The evolution of the critical shear stress is assumed to be linear:

g s,α = g α 0 + q ( H sq,α p q,α ) (52)
where g α 0 is the initial value of the critical shear stress, H α is the hardening matrix and p s, α is the accumulated plastic shear strain (i.e. ˙ p s,α = | ˙ γ s,α | ). Within a rate-independent plasticity framework, active slip systems are defined according to the following conditions:

˙

γ s,α = 0 , | τ s,α | = g s,α and | ˙ τ s,α | = ˙ g s,α (53) ˙ γ s,α = 0 , | τ s,α | < g s,α or | ˙ τ s,α | < ˙ g s,α (54)
At each time t , for active slip systems, the corresponding plastic shear strain rates ˙ γ s,α are determined to fulfill the consistency condition | ˙ τ s,α | = ˙ g s,α . According to the above constitutive relations, the set of material parameters for a given constituent consists of the independent thermoelastic constants associated with C α and α α , the initial critical shear stress g α 0 and the hardening matrix H α .

Representative volume element

To represent the pearlitic steel aggregate, the volume element presented in Fig. 1 is used. It has been constructed from a Voronoï tessellation and consists of 200 equiaxed grains with random crystallographic orientations. For the application of the FFT-based resolution method, the volume element is discretized into 128 3 voxels.

To describe pearlitic microstructures, three different assumptions are used: the Voigt approximation, the laminate structure approximation and the Mori-Tanaka approximation. For the application of the Mori-Tanaka approximation, ferrite is considered to be the matrix constituent and cementite lamellas are treated as penny-shaped inclusions with an aspect ratio of 20.

The numerical procedures used for computing the effective behavior of a volume element with prescribed boundary conditions are detailed in Appendix A .

Material parameters

The elastic constants associated with cementite, for which a purely elastic behavior was assumed, were extracted from the work of [START_REF] Henriksson | Carbides in stainless steels: results from ab initio investigations[END_REF] . For ferrite, both the {110} 111 and {112} 111 slip systems have been considered. The corresponding hardening parameters have been chosen to match the experimental results of [START_REF] Inal | Second-order stresses and strains in heterogeneous steels: self-consistent modeling and x-ray diffraction analysis[END_REF] who used diffraction techniques to determine the evolution of the average axial stress in each phase during an uniaxial tension test. The material parameters associated with both phases are listed in Table 1 .

Results

First, the effective response of the pearlitic polycrystalline aggregate has been computed for the case of uniaxial tension (up to an axial strain E 11 of 2%). The evolution of the macroscopic axial stress 11 as a function of the macroscopic axial strain E 11 is plotted in Fig. 2 . The evolution of the average stress is given for each metallurgical phase (ferrite and cementite) in Figs. 3 and4 . According to the results, the effective hardening rate in the plastic regime is much higher when the Voigt approximation is used. This difference between the Voigt approximation and the other approximations can be attributed to internal stresses, which are largely overestimated with the Voigt approximation. Indeed, while the behavior of the ferritic phase is correctly described with all methods, important discrepancies exist concerning cementite. More specifically, the bilinear type of behavior which is experimentally ob-Fig. 2. Evolution of the macroscopic axial stress 11 as a function of the macroscopic axial strain E 11 for a pearlitic steel polycrystalline aggregate subjected to uniaxial tension. Fig. 3. Evolution of the average stress in ferrite as a function of the macroscopic axial strain E 11 for a pearlitic steel polycrystalline aggregate subjected to uniaxial tension. The experimental results were obtained by [START_REF] Inal | Second-order stresses and strains in heterogeneous steels: self-consistent modeling and x-ray diffraction analysis[END_REF].

Fig.

of the average stress in cementite as a function of the macroscopic axial strain E 11 for a pearlitic steel polycrystalline aggregate subjected to uniaxial tension. The experimental results were obtained by [START_REF] Inal | Second-order stresses and strains in heterogeneous steels: self-consistent modeling and x-ray diffraction analysis[END_REF]. served for cementite is reasonably reproduced only when the laminate structure or the Mori-Tanaka approximations are used. For a given voxel, the Voigt approximation imposes the same strain state in both constituents. It thus fails in describing the fact that strains are preferably transferred to the soft constituent (ferrite) rather than the hard constituent (cementite). As a consequence, important internal stresses are generated in cementite because the imposed strain state is accommodated in a purely elastic manner. These internal stresses provide a supplementary contribution to hardening, resulting in an overestimation of the final axial stress in cementite. These results thus emphasize the importance of a proper description of morphological features when evaluating the response of heterogeneous materials with multiple length scales. The above computations have been performed on a DELL Pow-erEdge R710 server with 2 processors (Intel Xeon X5690) using 8 cores. The computation time required for performing uniaxial tension simulations is given for each method in Table 2 . The lowest computation time is obtained with the Voigt approximation, which does not involve any sophisticated numerical procedure. At the opposite, important computational efforts are needed for the application of the Mori-Tanaka approximation. Indeed, for this homogenization method, the determination of effective properties requires additional operations with a significant impact on the total computation time.

Finally, the response of the pearlitic polycrystalline aggregate has been calculated for plane strain loading conditions such that the imposed macroscopic strain state is:

E = +10% 0 0 0 0 0 0 0 -10% (55)
The loading conditions are chosen to represent a rolling operation with a reduction ratio of 10%. To show how rotations are impacted by the chosen homogenization method, the final {110} and {111} pole figures obtained for the ferritic phase after rolling are plotted in Fig. 5 . Though quite similar trends are obtained for all approximations, texture evolution is found to be slightly more pronounced when the Voigt approximation is used.

Conclusions

When using FFT-based techniques for modelling the behavior of heterogeneous materials with periodic microstructures, the common strategy consists of assuming the mechanical properties to be homogeneous within each voxel. This assumption is however not sufficient for microstructures with multiple length scales for which some specific strategies have to be developed. In the first part of this paper, different methods for considering multiple microstructural length scales have been proposed for materials with an inelastic behavior. More specifically, three different strategies, which use elementary homogenization rules to determine the effective behavior of heterogeneous voxels, have been described: the Voigt approximation, the laminate structure approximation and the Mori-Tanaka approximation.

Then, to demonstrate their interest, these methods have been used to calculate the effective behavior of a pearlitic steel polycrystalline aggregate. According to the numerical results, a specific care should be taken when using the Voigt approximation. Indeed, the Voigt approximation, which is an upper bound approximation, strongly overestimates macroscopic stresses and ignores the influence of morphological features. The impact of such features on the effective behavior can be considered by using either the laminate structure approximation or the Mori-Tanaka approximation. The application of the Mori-Tanaka approximation however requires more important computational efforts. [START_REF] Hielscher | A novel pole figure inversion method: specification of the MTEX algorithm[END_REF].

Appendix A. Implementation details

The global strategy used for computing the effective behavior of a volume element with prescribed boundary conditions is described in Table A.3 . The objective of the iterative procedure (from step 3a to step 3h) is twofold. First, it is needed to ensure that, once convergence is reached, the effective properties C and ¯ in used for steps 3c and d are the same as those determined from step 3f. Second, there are some situations where stress-controlled or mixed loading conditions are imposed to the volume element. However, whatever the boundary conditions are, step 3c uses the macroscopic strain and rotation tensors E and to solve the integral Eqs. ( 11) and ( 12) . The iterative procedure thus also aims at finding the macroscopic strain tensor E satisfying the imposed stress-controlled or mixed loading conditions.

Convergence is reached once the difference between the strain tensors obtained from two consecutive iterations is small enough. In the present work, a tolerance e c of 10 -8 has been used for determining whether convergence is reached or not.

Also, one should note that the iterative procedure, which is needed for the resolution of the integral Eqs. ( 11) and ( 12) , is contained in step 3c. In the present work, the resolution method pro- [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] to solve Eqs. ( 11) and ( 12) ¯ i and ω i d/ Use localization equations to determine the strain and rotation tensors within each constituent posed by [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] , which has been later reformulated by [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] , is used. A detailed description of the procedure is given in Table A.4 . For the calculations presented in this work, a value of 10 -9 is used for the convergence criterion e c .

Fig. 1 .

 1 Fig. 1. Periodic representative volume element used for the different calculations. The volume element is discretized into 128 3 voxels.

Fig. 5 .

 5 Fig.5. Final {110} and {111} pole figures for the ferritic phase of a pearlitic polycrystalline aggregate after rolling. The crystallographic orientations have been calculated with different approximations (Voigt, laminate structure or Mori-Tanaka). Pole figures have been plotted using the MTEX package[START_REF] Hielscher | A novel pole figure inversion method: specification of the MTEX algorithm[END_REF].

  the effective properties for each voxel C i and ¯ in i g/ Calculate the macroscopic stress state from relation (2) i h/ Evaluate the convergence indicator e : e = max V || α Description of the iterative procedure used for solving the integral Eqs. (11) and (12) with FFT-based methods. e ≥ e c : a/ Start a new iteration i = i + 1 b/ Calculate the polarization tensor field τi = C : ē i -1 -¯ in -C 0 : ē i -1 c/ Use the FFT algorithm to calculate the polarization tensor field in the frequency domain ˆ τi = FFT ( τi ) d/ Compute the convolution operation G sym * τi = FFT -1 ˆ G sym : ˆ τi e/ Propose a new estimation for the tensor field eē i = ē i -1 + 2 C 0 + C -1 : C 0 : ( E -ē i -1 -G sym * τi ) f/ Evaluate the convergence indicator e = max V ( || ē i -ē i -1 || ) 3Set the strain field ¯ equal to the tensor field ē ¯ = ē i 4Compute the convolution operation G skw * τi = FFT -1 ˆ G skw : ˆ τi 5Calculate the rotation field ω = -G skw * τi

Table 1

 1 Material parameters for the eutectoid steel grade. Only independent elastic constants are listed.

	Cementite
	C 11 (GPa)	C 12 (GPa)	C 13 (GPa)
	394	157	146
	C 22 (GPa)	C 23 (GPa)	C 33 (GPa)
	412	166	360
	C 44 (GPa)	C 55 (GPa)	C 66 (GPa)
	83	133	136
		Ferrite
	C 11 (GPa)	C 12 (GPa)	C 44 (GPa)
	231	135	116
	g 0 (MPa)	H ss (MPa)	H st (MPa)
	270	2 .2	2 .0

Table 2

 2 Total computation times required for performing uniaxial tension simulations for each method.

	Method	Computation time (s)
	Voigt	1 .27 × 10 5 s
	Laminate	2 .16 × 10 5 s
	Mori-Tanaka	2 .67 × 10 5 s

Table A .3

 A Description of the strategy used for computing the effective behavior of a volume element with FFT-based methods.

	Algorithm
	1	Set the state variables to their initial values
	2	For each time increment:
		i = 0
	3	Do while e ≥ e c :
		a/ Start a new iteration
		i = i + 1
		b/ Estimate the macroscopic strain and rotation tensors from
		boundary conditions
		E i and i
		c/ Use the iterative procedure of

No particular restriction applies to the origin of the inelastic contribution (e.g. thermal, plastic, viscoplastic).


Amine Ammar
email: amine.ammar@ensam.eu

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Rodolphe Heyd
email: rodolphe.heyd@univ-orleans.fr

Kevin H Knuth

Adom Giffin

Thermal Conductivity of Suspension of Aggregating Nanometric Rods

Keywords: nanofluids, thermal homogenization, aggregation kinetics, percolation PACS: 05.20.Dd, 05.10.Ln, 05.10.Gg

Enhancing thermal conductivity of simple fluids is of major interest in numerous applicative systems. One possibility of enhancing thermal properties consists of dispersing small conductive particles inside. However, in general, aggregation effects occur and then one must address systems composed of dispersed clusters composed of particles as well as the ones related to percolated networks. This papers analyzes the conductivity enhancement of different microstructures scaling from clusters dispersed into a simple matrix to the ones related to percolated networks exhibiting a fractal morphology.

Introduction

It is crucial today to explore technics of energetic and thermal performances enhancement in many fields of science to save on many energy resources. Many simple liquids (water, ethylene glycol, oils) are usually used as heat transfer fluids: in thermal engines (coolants, lubricating oils, etc.), in solar panels of individual water heaters, in the collectors of industrial solar power plants, in power transformers (industry, power distribution network), and in industrial installations (cutting oils, coolants). Unfortunately, these usual liquids all have low thermal conductivities at ambient or industrial temperature. We can find, for example, in [START_REF] Pop | Thermal properties of graphene: Fundamentals and applications[END_REF], some values of usual thermal conductivity of fluids: water, 0.58 W/(m•K); ethylene glycol, 0.18 W/(m•K), which are very low compared to usual metals such as: iron, 80 W/(m•K); copper, 400 W/(m•K); silver, 430 W/(m•K); graphite, 1000 to 1500 W/(m•K); carbon nanotubes, 2000 to 2500 W/(m•K); suspended graphene, 2500 W/(m•K) and diamond, 2500 W/(m•K). Any increase in the thermal conductivity promotes the exchanges of thermal energy and consequently the extraction ability of the heat from fluids.

The use of fluids containing metallic or semi-conductive particles in homogeneous and stable suspensions is one of the solutions studied for several years for the improvement of thermal transfer. Suspensions of microscopically sized conductive particles are contemplated and used for many years. However, due to the size of the suspended particles, they present problems that limit their industrial applications: abrasion and high risk of sedimentation and clogging. Moreover, the surface-to-volume ratio, which plays a predominant role in all conductive transport phenomena, can be significantly improved with particles of smaller sizes.

With the constant progress of nanotechnologies and the synthesis of specific nanoparticles, the idea of using thermal nanofluids has attracted a lot of interest in the scientific community [START_REF] Choi | Enhancing Thermal Conductivity of Fluids with Nanoparticles[END_REF], but also in industry. The scientific literature on the subject is vast and diverse. It covers both the fundamental (theoretical) aspects of the static thermal modeling, rheological properties of nanofluids and their synthesis, numerical modeling as well as the thermal and rheological characterization of suspensions in the laboratory.

In [START_REF] Xuan | Heat transfer enhancement of nanofluids[END_REF], a theoretical model is proposed to describe heat transfer performance of the nanofluids. Others works such as [START_REF] Milanese | An investigation of layering phenomenon at the liquid-solid interface in Cu and CuO based nanofluids[END_REF] are based on molecular dynamic simulations and aim to explain the differences in thermal conductivity between nanofluids based on metal (Cu) and metal oxide (CuO) nanoparticles. A large part of published works related to the thermal conductivities of nanoparticles suspensions have been based on experimental measurement. In [START_REF] Mintsa | New temperature dependent thermal conductivity data for water-based nanofluids[END_REF], effective thermal conductivity measurements of alumina/water and copper oxide/water nanofluids are performed and the effects of particle volume fraction, temperature and particle size were investigated. High aspect ratio fillers (carbon nanotubes (CNTs), silver nanowires, copper nanowires) are considered in [START_REF] Gu | Thermal conductivity of nanofluids containing high aspect ratio fillers[END_REF]. It is found there that silver nanofluid has the highest thermal conductivity compared with copper and CNT nanofluids. In [START_REF] Lacobazzi | An explanation of the Al 2 O 3 nanofluid thermal conductivity based on the phonon theory of liquid[END_REF], authors analyze several mechanisms affecting the thermal conductivity of Alumina based nanofluid, such as layering, Brownian motion, clustering, ballistic phonon motion, thermal boundary resistance and mass difference scattering. In [START_REF] Sundar | Experimental thermal conductivity of ethylene glycol and water mixture based low volume concentration of Al 2 O 3 and CuO nanofluids[END_REF], authors look for the specific case when the base fluid is a mixture of ethylene glycol and water. Thermal conductivity has been estimated experimentally at different volume concentrations and temperatures. Other works are related to the coupling with kinematics and rheological response. The effect of viscosity is studied in [START_REF] Colangelo | Thermal Conductivity, Viscosity and Stability of Al 2 O 3 -Diathermic Oil Nanofluids for Solar Energy Systems[END_REF]. The kinematic effects are also studied in [START_REF] Nnanna | Experimental model of temperature-driven nanofluids[END_REF], where an experimental method of studying the heat transfer behavior of buoyancy-driven nanofluids was addressed. These studies open opportunities for many application fields. For example, an interesting applications, related to wind turbines and solar collectors, is found, respectively, in [START_REF] De Risi | High efficiency nanofluid cooling system for wind turbines[END_REF][START_REF] Yousefi | An experimental investigation on the effect of Al 2 O 3 -H 2 O nanofluid on the efficiency of flat-plate solar collectors[END_REF]. The reader can find in [START_REF] Lomascolo | Review of heat transfer in nanofluids: Conductive, convective and radiative experimental results[END_REF] and in the references therein an overview of experimental results about the heat transfer capabilities of nanofluids, which is presented using widely scattered available information from diverse literature sources.

The main aim of the present work is the analysis of thermal properties, in particular the equivalent thermal conductivity, of suspensions involving manometric rods suspended in a Newtonian (purely viscous) fluid. When using particles highly conductive in a poorly conducting matrix, even at very low particle concentration (less than 1 vol.%), it has been observed that the resultant equivalent fluid (suspension) increases (sometimes more than 100%) its apparent thermal conductivity significantly.

In literature, many phenomena are found to be responsible for thermal conductivity enhancement: Brownian motion, ordered liquid layer formation at the nanoparticle surface, interfacial thermal resistance, fractal structuring, cluster formation, thermal conductivity of particles, viscosity of the host fluid, sonication, and surface treatment of the nanoparticles, among many others. In the vast literature about the thermal properties of nanofluids, we can also find a large dispersion of the experimental results. For example, in [START_REF] Philip | Thermal properties of nanofluids[END_REF], it is mentioned, in regard to the effect of the particle volume fraction for example, that enhancements within and beyond the effective medium theory (random dispersion of the nanoparticles through the host fluid) are reported in the literature. The effects of the nanoparticle sizes and nanofluid temperature have also been analyzed in this review; they also show a vast dispersion in the enhancements of the thermal conductivity of nanofluids.

The reasons for these contradictory reports are certainly multiple. We can mention, for example: the choices of the nanoparticle surface chemical treatments; the choice of the host fluid and of the method of the nanoparticle dispersion; and, even sometimes, the choice of the measurement technique.

Taking into account all of these aspects, we have chosen to adopt the following point of view, which consists of analyzing the sole influence of the organization of the nanoparticles within the nanofluid, on its thermal conductivity. The Brownian motion must be viewed in this study as a way to construct two realistic types of structured patterns: fractal skeletons (as observed experimentally and mentioned in [START_REF] Philip | Nanofluid with tunable thermal properties[END_REF]) and micrometric aggregates (as also observed experimentally), and not as a way to promote heat transfer at the microscopic scale.

Thus, we try in this work to answer the following fundamental question: to what extent is the organization of the nanoparticles inside nanofluids essential for improving the thermal conductivity of nanofluids?

This question is fundamental in many aspects. We will prove that the effect of organization is crucial in the enhancement. In addition, the percolation or not of the particles is determinant for final conductivity. This is demonstrated by calculating the heat transfer homogenization based on Fourier's law without taking into account the Brownian induced conductivity.

In most existing works, the values of the volume fraction considered was much lower than the values used in the present paper. As it has been previously indicated, the present work aims to show the effect of the structure on the thermal conductivity enhancement. In many experiments concerning low viscosity liquids, the Brownian induced conductivity is at the origin of the enhancement of thermal performances at low concentration. In this case, the simulations have to be done with another approach based on the Green-Kubo method, which is used to obtain the thermal conductivity of the nanofluid from the heat autocorrelation function (see for example [START_REF] Jain | Brownian dynamic simulation for the predictions of effective thermal conductivity of nanofluids[END_REF][START_REF] Jia | Investigation of the characteristics of heat current in a nanofluid based on molecular dynamics simulation[END_REF]). Our strategy in the present work is different. We address heat transfer from the homogenized thermal conductivity related to conduction (Fourier's law). Such an approach does not take into account Brownian motion for calculating the effective conductivity. Brownian motion has only been used here for developing microstructures (random fractals or clusters) as experimentally observed in many situations where nanoparticles agglomerate inside the host fluid. Moreover, the approach here proposed could be useful to explain thermal conductivity enhancements for some fluids where the viscosity effects are preponderant. In fact, in [START_REF] Xuan | Heat transfer enhancement of nanofluids[END_REF][START_REF] Estelle | Thermal conductivity of CNT water based nanofluids: Experimental trends and models overview[END_REF], we observe almost the same enhancement at a low concentration. Even if it is shown in [START_REF] Murshed | Superior thermal features of carbon nanotubes-based nanofluids-A review[END_REF] that high enhancements of the thermal conductivity can be found for low concentrations, compared to our simulations, it is not in contradiction with our study because we only focused here on structuration effects and not on the other different factors previously mentioned that could promote an extra-thermal conductivity enhancement.

It is well known that particles, suspended into a viscous medium, tend to create aggregates that could reach larger space scales, and, in some cases, they are resolvable at the eyes scale. The modeling of such aggregates requires the consideration of kinetic models at the meso-scale accounting for the flow induced aggregation/disaggregation mechanisms [START_REF] Ma | Rheological modelling of carbon nanotube aggregate suspensions[END_REF]. As soon as these particles are functionalized in order to avoid aggregation mechanisms, no aggregates are observed at the resolvable scales, and standard models developed for addressing dilute suspensions work well for describing their main rheological features [START_REF] Ma | The rheology and modelling of chemically treated carbon nanotube suspensions[END_REF].

In many cases when addressing materials or meta-materials (functionalized materials exhibiting gradients of properties originated by different concentrations and orientations of the dispersed charges) forming processes, after dispersing the particles, a microstructure develops inside and can be frozen by transforming the fluid phase into a solid (e.g., by polymerizing a thermoset resin).

The final properties of the resulting microstructured fluid or solid become radically different from a simple mixing rule as soon as percolation occurs. In these cases, an appropriate homogenization strategy is compulsory for evaluating the final properties. There are numerous works addressing homogenization strategies for systems consisting of perfectly dispersed particles into a matrix (i.e., in the absence of aggregates). In [START_REF] Perez | On the multi-scale description of electrical conducting suspensions involving perfectly dispersed rods[END_REF], a reinforced polymer consisting of a population of N nanometric conducting rods immersed into an insulating polymer matrix was assumed. At the microscopic scale, the system can be described from the specification of the position of the centre of gravity G i and orientation p i of each rod in the population, with i = 1, • • • , N .

The microscopic description, despite the fact of being the finest one, is too heavy from both the computational and the experimental points of view. For this reason, coarser descriptions are preferred. At the mesoscopic scale, the rods' individuality is substituted with an orientation distribution function Ψ(x, p) given the fraction of rods that at position x align along direction p. Because this probability distribution function-pdf-is defined in a space of high dimensionality (involving the physical space x and the conformation coordinates p), the solution of its associated balance equation (the so-called Fokker-Planck equation) remains cumbersome, even if advanced numerical discretization techniques have been proposed for alleviating or circumventing the curse of dimensionality that high dimensionality entails [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF].

Finally, macroscopic descriptions constitute the coarsest level of description but the one leading to faster simulations. This description is based on the use of the probability distribution function moments, also known as orientation tensors, being the expression of the second moment given by a(x) = p ⊗ pΨ(x, p) dp.

(

) 1
In order to quantify the rod network, we introduced in [START_REF] Perez | On the multi-scale description of electrical conducting suspensions involving perfectly dispersed rods[END_REF] the density of rod contacts C(x, p) for a rod with orientation p, depending on the two main microstructure descriptors: (i) the rods concentration φ(x); and (ii) the orientation distribution Ψ(x, p). Its evaluation follows the simple geometrical construction depicted in Figure 1 in which first approximation rods are assumed rigid and straight. Similar descriptors were proposed in [START_REF] Doi | Dynamics of rod-like macromolecules in concentrated solution. Part 1[END_REF][START_REF] Ranganathan | Fiber-fiber interactions in homogeneous flows of nondilute suspensions[END_REF][START_REF] Toll | On the tube model for fiber suspensions[END_REF][START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF]. It can be expected that the conductivity at position x along direction p scales with the number of contacts φ(x)C(x, p)Ψ(x, p). Thus, the conductivity becomes local and directional. One can expect that percolation along direction p occurs locally (at position x) when the number of contacts is higher than a threshold value T . Lower values imply an infinite local directional resistivity. Thus, percolation could be considered local and directional, allowing us to create a network in which to evaluate the effective conductivity.

When considering electrical properties, the calculation at the network level was carried out by using usual Kirchhoff's laws [START_REF] Perez | On the multi-scale description of electrical conducting suspensions involving perfectly dispersed rods[END_REF]. The question arising immediately concerns the existence of an appropriate fully-macroscopic description modeling thermal conductivity in a perfectly dispersed suspension composed of rods. This modeling framework, with some simplifying assumptions, is described in Section 2. Then, we will focus on more complex scenarios, the ones involving aggregation mechanics leading to the creation of macroscopic clusters (dispersed into the suspending matrix) and percolated networks, both addressed in Section 3. The computational thermal homogenization of such complex microstructures is addressed in Section 4 before then discussing some numerical results in Section 5.

Different assumptions are considered in the sequel: (i) advective heat fluxes are neglected when the matrix becomes solid in the final part in which the thermal properties are evaluated; (ii) the suspending fluid is assumed viscous and at rest, making possible the particle movement, but remaining insensible to thermal dependencies of its density when addressing thermal homogenization, the viscous forces being preponderant with respect to buoyancy, in all the cases considered here; and (iii) the thermal conductivity contrast between the matrix and the particles is assumed to be large.

On the Macroscopic Thermal Conductivity of Perfectly Dispersed Rod Suspensions

When the suspension involves perfectly dispersed rods, a homogenized thermal conductivity can be derived almost analytically. This model only involves the orientation distribution second moment, also called the "orientation tensor" introduced above. In the case of more complex particle shapes, computational homogenization (described in Section 4) represents a valuable alternative for evaluating the equivalent conductivity. When individual particles of any shape aggregate to create more complex microstructures, computational homogenization seems compulsory for calculating homogenized properties as discussed later.

Before addressing more complex scenarios, we start our analysis by describing the homogenization procedure in the case of a perfectly dispersed suspensions of rods, assuming no thermal interfacial resistance between the rods and the matrix as previously indicated, while ignoring advective heat transfer.

The matrix thermal conductivity is assumed, for the sake of simplicity and without loss of generality, isotropic, and expressed by the scalar k m , whereas the bulk isotropic conductivity of the rod particles is expressed by k p , with, as indicated above, k p k m . The thermal contrast c is defined from c = k p /k m , with, in our case, c

1. The heat flux Q m through the host matrix is calculated according to the standard Fourier law

Q m ≈ -(1 -φ)k m ∇T, (2)
where ∇T is the given macroscopic gradient of temperature, supposed to be constant within the representative volume element-RVE-ω. Now, in a given rod with orientation p, the heat flux is assumed directed in the rod direction (we are assuming high thermal contrast), and, consequently, the effective temperature gradient along the rod is written as

(∇T • p)p = (p ⊗ p) • ∇T, (3)
and, consequently, the heat flow inside the rod reads

q p (p) ≈ -k p (p ⊗ p) • ∇T. (4)
By adding the contribution of all the rods in ω, it results

Q p = S q p (p)φΨ(p)dp ≈ -k p φ a • ∇T. (5)
Now, with the total heat flux being

Q = Q m + Q p , we obtain Q = -k m (1 -φ)I + k p φ a • ∇T = -k m ((1 -φ)I + cφa) • ∇T, (6)
with I the identity tensor. Thus, we can define the homogenized thermal conductivity tensor K as

K = -k m ((1 -φ)I + cφa) . ( 7
)

Generating Microstructures and Networks

In this section, we consider different ways of creating more complex microstructures from individual particles when clustering mechanisms act. In particular, we are considering two different algorithms: (i) the DLAa (Diffusion Limited Aggregation algorithm) and (ii) the MDAa (Molecular Dynamics Aggregation algorithm)-both operating at different levels of concentration.

In all the cases studied in this paper, we performed different runs of simulation, in order to extract the average and variability of the computed results.

3.1.

From DLA to Enhanced-DLA First, we start with a fully empty space and put one particle somewhere. Now, a second particle is randomly introduced and it keeps moving randomly until it reaches the first one. Then, both are attached and assumed at rest while a third particle is introduced and randomly moves until reaching the cluster composed by the previously attached particles. The procedure repeats until all particles have been introduced and been attached to the particle network.

The main issues of such a strategy are that it does not address conveniently the fact that the systems proceed from a given constant particle concentration and the fact of addressing Brownian motion more consistently. To better address both questions, we propose considering a lattice model that contains all possible particle locations and defines consequently the discrete particle jumps during their motion. The process parameters are: (i) the particle concentration C, the diffusion D quantifying the Brownian motion and the characteristic length L of the representative volume element-RVE-with respect to the particle equivalent diameter. Thus, for example, L = 10 means that the RVE size is ten times the particle diameter. For non-spherical particles, an equivalent diameter can be defined from excluded volume considerations.

Thus, the enhanced-DLA proceeds by putting all the particles into the simulation box (RVE), a number that is determined from the concentration. One of the simulated particles is located at the center of the computational box and it is kept at rest all throughout the simulation. The remaining particles, all except the one attached at the centre of the RVE, move randomly. Periodic boundary conditions are assumed on the RVE boundaries, implying that one particle leaving the box from one of its faces, comes in from the opposite face. As soon as a moving particle collides with the one kept at rest at the box center, it is glued to the first and occupies the neighbor lattice position (associated with the arrival direction). Then, both (the initial one and the just glued) are kept at rest while waiting for the arrival of the third particle, and so on. Thus, the particle network continues growing until all particles join it.

The enhanced DLA algorithm is illustrated in Figure 2 for D = 5. In this figure, black particles are those belonging to the network assumed at rest. Around three moving particles (numbered white circles) are depicted, and, at the present time, all of them move D lattice units (D = 5 in this case), making directional random choices at each of the D -1 lattice jumps (we do not consider the possibility of coming back during a lattice jump). In what follows, unit jumps on the lattice will be denoted by N (north), S (south), E (east) or W (west). In the random particle movement, some constraints apply:

• The random path of a particle does not cross any other free (moving) particles.

• If a random path ends at a gluing point, the last being a neighbor of the already attached particles (black circles), symbolized by a cross in Figure 2, then the particle is attached to the growing network and kept at rest in the sequel (particle 2 in Figure 2). As soon as the particle joins the network, the set of gluing points is updated by including the neighbors of the just glued particle on the lattice. • If the random path crosses a gluing point, it is attached to the network, even if it does not complete its D -1 unit jumps (see particle 1 in Figure 2). • If the random path ends far enough from the growing network, it remains free and then ready to randomly move at the next time step (particle 3 in Figure 2).

The simulation finishes when all particles have been attached to the network. Figures 3 and4 present the resulting microstructures consisting of a single network composed of all the particles, when considering different values of the simulation parameters D and C. It can be noticed that the modification of the diffusion parameter D (from D = 1 to D = 16) does not affect the nature of the resulting fractal structure (the fractal dimension remains unchanged); however, its increase results in a faster development of the particle network. On the other side, by modifying the concentration, the fractal structure is significantly affected. Thus, for discussing the steady-state properties, in what follows, only the concentration will be considered.

Molecular Dynamic Aggregation Algorithm

When considering the so-called molecular dynamics aggregation algorithm-MDAa-particles are randomly placed into the RVE and a two-particle interaction potential (the Lennard-Jones potential in the simulations carried out later) is considered acting between each pair of particles i and j, according to

V(d) = 4ε σ d 12 - σ d 6 , ( 8
)
where d is the distance between the centers of two particles involved in the interaction, and σ and ε are two model parameters.

Based on this potential, the interaction force acting on the i-particle induced by the one located at position x j reads

F j→i = -∇V| d(i,j) n j→i , (9)
where d(i, j) is the mutual distance and n j→i the unit vector that points from x j to x i . Thus, the resulting force acting on particle x i is written as

F i = ∑ j =i F j→i . ( 10
)
The algorithm consists of updating particle velocity and position by using the dimensionless Newton's equation of dynamics, which allows for computing the particle acceleration

F i (t) = m i d 2 x i dt 2 t = m i a i (t), (11)
and then velocity and position from an adequate (simplectic) integration strategy, as the one of Verlet that reads

x i (t + ∆t) = x i (t) + v i (t)∆t + 1 2 a i (t)∆t 2 , ( 12
)
v i (t + ∆t) = v i (t) + 1 2 (a i (t) + a i (t + ∆t)) ∆t. ( 13
)
When using the MDAa, after introducing randomly all the particles into the RVE, interaction potential applies and particles start moving, creating different clusters within the RVE, with positions ands size evolving in time, instead of the single particle network created from DLA-based strategies.

Calculating the Equivalent Thermal Conductivity

Due to the microscopic heterogeneity (where matrix and particles coexist, the last cluster making part of a single network depending on the microstructure constructor chosen, MDAa or DLAa, respectively, the macroscopic thermal modeling needs a homogenized thermal conductivity that depends on the microscopic details.

We can define the macroscopic temperature gradient ∇T at the RVE level from

∇T = g(x) = 1 |ω| ω g(x) dx, (14)
where g(x) denotes the microscopic temperature gradient. We also assume the existence of a localization tensor relating the macroscopic temperature gradient (average of its microscopic counterpart) ∇T with the microscopic local temperature gradient g(x), L(x), such that

g(x) = L(x) • ∇T. (15)
Now, we consider the microscopic heat flux q(x) at x ∈ ω according to Fourier's law

q(x) = -k(x) • g(x), (16)
and its macroscopic counterpart Q that is written

Q = q(x) = -k(x) • g(x) = -k(x) • L(x) • ∇T, (17)
from which the homogenized thermal conductivity can be defined at the RVE level from

K = k(x) • L(x) . ( 18
)
As k(x) is perfectly known everywhere in the representative volume element ω, the definition of the homogenized thermal conductivity tensor only requires the computation of the localization tensor L(x). Several approaches are proposed in the literature to define this tensor, according to the choice of the boundary conditions. Our objective here is not to discuss such a choice. For the sake of simplicity, we use essential boundary conditions on ∂ω corresponding to the assumption of a uniform temperature gradient on the RVE ω. We consider the general 2D case that involves the solution of the two boundary value problems related to the steady state heat transfer model in the microscopic domain ω for two different boundary conditions on ∂ω:

∇ • k(x) • ∇Θ 1 (x) = 0 Θ 1 (x ∈ ∂ω) = x , (19)
and

∇ • k(x) • ∇Θ 2 (x) = 0 Θ 2 (x ∈ ∂ω) = y . (20)
It is easy to prove that both solutions verify

∇Θ 1 (x) = (1, 0) T ∇Θ 2 (x) = (0, 1) T , (21)
where (•) T denotes the transpose. Thus, the localization tensor finally results as

L(x) = ∇Θ 1 (x) ∇Θ 2 (x) . (22)
From a computational point of view, the main difficulty concerns the solution of the two boundary value problems using very fine meshes required to represent accurately all of the microscopic details. In the extreme case of meshes that are too fine, advanced solvers based on the use of separated representations (at the heart of the Proper Generalized Decomposition (PGD)) could be a valuable alternative [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Lamari | Routes for efficient computational homogenization of non-linear materials using the proper generalized decomposition[END_REF].

To check the proposed strategy just described, we consider the DLA-based microstructure depicted in Figure 5. Figure 6 shows both temperature fields Θ 1 (x) and Θ 2 (x).

For the solutions carried out in this study, the mesh size was chosen in order to ensure the numerical convergence, i.e., the mesh-size independence of the computed solutions. Figures 7 and8 analyze the solution dependence as a function of the number of nodes along the particle diameter, for one and many particles, respectively. It can be seen that when considering ten nodes along the particles' characteristic length, the conductivity becomes almost constant.

Numerical Results

In order to compare the homogenized thermal conductivity for different microstructures obtained by using the enhanced DLA and MDA algorithms, we proceed from the same initial particle distribution and different values of the associated concentration. The mesh considered ensured the numerical convergence of all of the thermal analyses carried out in the present study.

The value of concentrations used in the field of nanofluids is usually lower than the maximum range that we have considered in our study. It is generally on the order of 0.1% or to 2%. However, as mentioned previously, the dispersion of the experimental results is such that, in many cases reported in the literature, low volume concentrations of nanoparticles do not always lead to spectacular enhancements of the thermal conductivity. In some specific studies, the volume fraction can reach even values up to 9% [START_REF] Mintsa | New temperature dependent thermal conductivity data for water-based nanofluids[END_REF][START_REF] Philip | Nanofluid with tunable thermal properties[END_REF]. In the present study, our aim is to focus mainly on the structure-based enhancement effects. For lower concentrations, other numerical strategies are more convenient (the Green-Kubo method and heat autocorrelation function), but, as just indicated, this is not our purpose. In order to compare with some papers in literature, Figure 9 presents an extension of our study results to a volume fraction lower than 2%. The reported values show the increase of relative conductivity k suspension /k particle from 0.001 when there are no particles to 0.014 at 2% of volume fraction. These orders of magnitude are the same that are found in some existing works (e.g., [START_REF] Xuan | Heat transfer enhancement of nanofluids[END_REF][START_REF] Estelle | Thermal conductivity of CNT water based nanofluids: Experimental trends and models overview[END_REF]). Regarding this figure, modification of conductivity is almost the same for random distribution or MD simulation (they almost give the same values of conductivity because low concentration does not promote percolation). However, results of DLA are not significant here as they are related to a small aggregated structure in the middle of the simulation domain.

Figure 10 reports the results found for higher concentrations. It can be seen that, as expected for both kinds of microstructures (the one generated by the enhanced-DLAa and the MDAa), the equivalent (homogenized) conductivity increases with the particle concentration. Moreover, the thermal conductivity of the resulting microstructures is larger than the one related to the initial random distribution of particles, all of them being quite similar at very low concentrations. When the concentration increases, the percolated network development seems the most favorable for heat conduction. These numerical results are very similar to the experimental results reported in [START_REF] Philip | Nanofluid with tunable thermal properties[END_REF] on the basis of measurements of the thermal conductivity of magnetically polarizable nanofluids. As can be seen in Figure 10, at very large concentrations, our computations show that the MD aggregated microstructures are more favorable to heat conduction than the DL ones. This is a quite interesting result, in the sense that the agglomeration of nanoparticles in the form of individual clusters is not an obstacle to the heat conduction through the colloidal suspensions, as long as the sedimentation can be neglected. This facility can be easily understood by comparing the analyzed microstructure morphology depicted in Figures 111213, where the superiority of the percolated MD networks can be easily pointed out. For each concentration, we start with a random distribution for the first calculation of homogenized conductivity. Then, we let the random distribution evolve according either to the DLA or the DMA. The DLA tends to build as quickly as possible a network that looks like a fractal. During the simulation, and especially for high concentration, the fractal characterization seems to not be affected by the increase of the network. The topology of the network increases the conductivity by constructing preferential directions for heat conduction. Now, when the structure evolves according to the DMA, the observations are similar to what can be found in [START_REF] Lamas | Critical analysis of the thermal conductivity models for CNT based nanofluids[END_REF]. In fact, a homogeneous nanofluid is described as the dispersion of a discontinuous phase in a continuous one. However, the phase distribution can present three distinct geometric morphologies depending on the volume fraction and shape of the dispersed nanoparticles. These arrangements strongly modify the effective thermo-physical properties at the macro scale. The three geometrical morphologies can be defined as dispersed structure, aggregated structure, and percolated structure. The percolation threshold is defined as the minimum volume fraction that characterizes the formation of a continuous solid path. At low particle concentrations (less than 1%), it can be seen from our calculations (see Figure 9) that the organization of particles (in the form of clusters or fractals) has no noticeable effect on the thermal conductivity. The results are even comparable to those obtained with the effective medium theory-EMT-(random distribution of individual nanoparticles). Therefore, any unusual enhancement of the thermal conductivity at a low particle fraction must be explained by other arguments than structuration of the nanoparticles in the form of fractals or clusters. The results are quite different at higher particle fractions. It can be seen from our calculations (see Figure 10), for concentrations comprised of between 5% and 20%, that the structuration in the form of fractals (DLA) leads to thermal conductivity enhancements well above what is predicted by the EMT. Unfortunately, the formation of a fractal skeleton within the fluid may have a tragic influence on the rheological properties of the nanofluid, which could even make it unusable from a mechanical point of view. On the other hand, what seems very interesting to us is the positive influence of cluster formation on the thermal conductivity enhancement at fractions higher than 2, well above EMT and even fractal structuration. This is a very promising result in terms of nanofluid industrial applications because, in our opinion, the clusters' formation has a less damaging influence on the rheological properties of the nanofluid than the fractal skeletons (the validation of this constitutes a work in progress).

Conclusions

This paper analyzed thermal properties of microstructured fluids and solids, in which the initially dispersed conductive particles aggregate in order to create conductive percolated networks. Two microstructural constructors are analyzed, an enhanced DLAa and the MDAa, from the point of view of the morphology and the associated thermal properties. It has been proved that conductivity is significantly enhanced by the particles structured with respect to a perfectly random particle distribution.

Figure 1 .

 1 Figure 1. Evaluating rod contacts.

Figure 2 .

 2 Figure 2.Illustrating the DLA algorithm. Each particle moves D = 5 lattice units (horizontal and vertical coordinates refer to lattice units), in random directions, except if, before completing the D jumps, it reaches the aggregate (particle 1).

Figure 3 .Figure 4 .

 34 Figure 3. DLA-simulated microstructures when varying parameter D for C = 0.05.

Figure 5 .Figure 6 .

 56 Figure 5. DLA-based microstructure (again horizontal and vertical coordinates refer to lattice units).

Figure 7 .

 7 Figure 7. Convergence analysis in a volume containing a single particle. (a) Convergence curve; (b) Microstructure.

Figure 8 .

 8 Figure 8. Convergence analysis when considering many particles. (a) Convergence curve; (b) Microstructure.

Figure 9 .

 9 Figure 9. Evolution of equivalent conductivity with the particle concentration for low concentrations.

Figure 10 .

 10 Figure10. Evolution of equivalent conductivity with the particle concentration at higher concentrations.

Figures 11 -

 11  are generated using three concentrations equal to 0.1, 0.2 and 0.3, respectively. For each concentration, we start with a random distribution for the first calculation of homogenized conductivity. Then, we let the random distribution evolve according either to the DLA or the DMA. The DLA tends to build as quickly as possible a network that looks like a fractal. During the simulation, and especially for high concentration, the fractal characterization seems to not be affected by the increase of the network. The topology of the network increases the conductivity by constructing preferential directions for heat conduction. Now, when the structure evolves according to the DMA, the observations are similar to what can be found in[START_REF] Lamas | Critical analysis of the thermal conductivity models for CNT based nanofluids[END_REF]. In fact, a homogeneous nanofluid is described as the dispersion of a discontinuous phase in a continuous one. However, the phase distribution can present three distinct geometric morphologies depending on the volume fraction and shape of the dispersed nanoparticles. These arrangements strongly modify the effective thermo-physical properties at the macro scale. The three geometrical morphologies can be defined as dispersed structure,

Figure 11 .

 11 Figure 11. Initial random distribution and the resulting enhanced-DLA and MDA-based microstructures for C = 0.1. (a) Initial; (b) enhanced-DLAa; (c) MDAa.

Figure 12 .

 12 Figure 12. Initial random distribution and the resulting enhanced-DLA and MDA-based microstructures for C = 0.2. (a) Initial; (b) enhanced-DLAa; (c) MDAa.

Figure 13 .

 13 Figure 13. Initial random distribution and the resulting enhanced-DLA and MDA-based microstructures for C = 0.3. (a) Initial; (b) enhanced-DLAa; (c) MDAa.

Acknowledgments: Francisco Chinesta acknowledges the financial support of the ESI GROUP though its chair at the Ecole Centrale de Nantes.

Author Contributions: Amine Ammar was in charge of the proposal and implementation of the different microstructure constructors, Francisco Chinesta proposed the homogenization framework and Rodolphe Heyd defined the general physical framework. All authors have read and approved the final manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.


Amélie Buirette

Aurélie Halipré

Rachid Berbache

Joachim Schöpfel

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

National

The ANRT contribution to electronic theses and dissertations can be described on six levels.

  photos


Marie Chancel

Caroline Landelle

Caroline Blanchard

Olivier Felician

Michel Guerraz

Anne Kavounoudias

Mul- Tisensory Reweighting

Multisensory reweighting for kinesthesia in older adults

Submission guidelines for the International workshop on Aging in the Neuro-musculo-skeletal System, 15-17 March, 2016, in Marseille (France)

(1) Title of the abstract: Multisensory reweighting for kinesthesia in older adults

(2) Name of the authors: Marie CHANCEL 1 , Caroline LANDELLE 1 , Caroline BLANCHARD 1,[START_REF] Blanchard | Differential contributions of vision, touch and muscle proprioception to the coding of hand movements[END_REF] , Olivier FELICIAN [START_REF] Goble | Proprioceptive sensibility in the elderly: degeneration, functional consequences and plastic-adaptive processes[END_REF] , Michel GUERRAZ 4 , Anne KAVOUNOUDIAS 1 (3) Author affiliations: To perceive self-hand movements, the central nervous system (CNS) relies on multiple sensory information mainly derived from vision, touch, and muscle proprioception. Using psychophysical approaches, this study investigated how and to what extent the CNS relies on these sensory systems to build kinesthetic percepts when they all decline such as when aging. Illusory movement perceptions were induced by stimulating these three sensory modalities either separately or concomitantly. The perceptual responses reported by 19 healthy elderly adults (60-82 yrs) were compared to those of 12 young adults. Results suggest that reliance on sensory inputs for kinesthetic purposes is profoundly reshaped as early as 60 years old. Older people rely more on visual and tactile afferents for perceiving self-hand movements than younger adults. This could be due to a relative greater muscle proprioception impairment.

Introduction

To perceive self-hand movements, the CNS relies on multiple sensory information mainly derived from vision, touch, and muscle proprioception 1,[START_REF] Blanchard | Differential contributions of vision, touch and muscle proprioception to the coding of hand movements[END_REF] . Aging is associated with decline in all these sensory systems. Impairments of the somatosensory system, including touch and muscle proprioception, have been well described through numerous neurophysiological studies conducted in both animals and humans (see review [START_REF] Goble | Proprioceptive sensibility in the elderly: degeneration, functional consequences and plastic-adaptive processes[END_REF] ). Studies show alterations, in the elderly, of structural properties and density of mechanoreceptors, as well as peripheral and central nerve conduction. In the central nervous system, structural alterations such as the reduction of neuronal size, the number of synapses and the grey matter volume 4 with advancing age have been clearly demonstrated. These peripheral and central sensory damages may account for the functional deficits shown in older individuals, such as a decreased ability to perceive a movement [START_REF] Skinner | Age-related decline in proprioception[END_REF] , or to detect a tactile stimulus applied on the skin surface [START_REF] Desrosiers | Age-related changes in upper extremity performance of elderly people: a longitudinal study[END_REF] . Regarding visual motion perception, research shows that older observers are worse to discriminate speed and direction [START_REF] Bennett | The effects of aging on motion detection and direction identification[END_REF] . Most of the previous studies examined the alteration of each sensory source in isolation. The present work investigates to what extent such age-related plurimodal impairments are associated with a sensory reweighting to build kinesthetic percepts, and/or whether enhancement of multisensory integration could at least partly overcome the sensory decline.

Methods

Illusory sensations of right hand rotation were induced by stimulating separately or simultaneously the three modalities at two intensity levels. For this purpose, mechanical vibrations were applied to the pollicis longus muscle group in the subjects' wrists (proprioceptive stimulation), and a textured disk was rotated under the palmar skin of the subjects' right hands (tactile stimulation) while a background visual scene was projected onto the rotating disk (visual stimulation). The elicited kinesthetic illusions were copied by the subjects in real time. The perceptual responses of 19 healthy elderly adults (60-82 yrs) were compared to those of 12 young adults.

Results

1. At each level of intensity, the younger adults reported similar velocity and latency of hand rotation illusions in either modality stimulated. The older adults reported more salient and faster illusions for the visual and tactile conditions than for the muscle proprioceptive one. 2. Compared to those in the younger group, visual and tactile illusions were significantly more salient and faster in the older group. The onset latency of the illusions increased gradually with age. In contrast, the vibration-induced illusions were significantly less frequently evoked, less salient and more delayed in older adults. 3. For the three modalities confounded, increasing the intensity level of stimulation resulted in a smaller increase of illusion velocity in older adults than in younger adults. 4. Lastly, the velocity increase of illusory movement perception observed during multimodal stimulation in the elderly group was of the same order as that observed in the younger group.

Discussion

Very few studies have investigated the possibility of adaptive responses developed by the central nervous system in the elderly to compensate for multisensory system deterioration. A first hypothesis is that sensory system impairments in elderly people may be partly overcome by improving the central integrative processing of multisensory inputs. In a previous study, Liaurenti and colleagues [START_REF] Laurienti | Enhanced multisensory integration in older adults[END_REF] found that during visual and auditory discrimination tasks older individuals may take greater advantage of redundant audio-visual stimuli than younger adults by increasing the efficiency of the integrative processing. Although such increase was not evidenced in the present results, the multisensory benefit was found similar between older and younger participants. Alternatively, we found that the weighting of sensory modalities might change with aging: elderly people might rely more on visual and tactile afferents for perceiving self-hand movements than younger adults due to a relative greater impairment of muscle proprioception.

Conclusion

All in all, this study shows that reliance on sensory inputs for kinesthetic purposes is profoundly reshaped as early as 60 years old. Future studies should be conducted to confirm whether such multisensory reweighting is functionally adaptive.

1 . 2 .( 6 )

 126 Aix-Marseille Univ. CNRS, LNIA UMR 7260, Marseille, France University of Nottingham, School of Psychology, NG7 2RD Nottingham, UK 3. Aix-Marseille Univ, Inserm, INS, Marseille, France 4. Univ. Savoie Mont Blanc, CNRS, LPNC, F-73000 Chambéry, France (4) Address and contact information of the first author : Laboratoire de Neuroscience Intégrative & Adaptative (UMR7260) -3 Place Victor Hugo -13003 Marseille. Tel.: +33 4 13 55 03 65. Email: marie.chancel@univ-amu.fr (5) Type of communication: Oral presentation Poster Key words (5 maximum): Movement perception, Illusions, Muscle proprioception, Vision, Touch (7) Summary (10 lines maximum):

Figure 1 :

 1 Figure 1: Experimental set-up

Multisensory reweighting for kinesthesia in older adults

Marie CHANCEL 1 , Caroline LANDELLE 1 , Caroline BLANCHARD 1,2 , Olivier FELICIAN 3 , Michel GUERRAZ 4 , Anne KAVOUNOUDIAS 1