Étienne Malachanne

David Dureisseix

Patrick Cañadas

Franck Jourdan

Experimental and numerical identification of cortical bone permeability

Keywords: Permeability, Poroelasticity, Cortical bone, Structural simulation

Bone is a complex system, and could be modeled as a poroelastic media. The aim of this paper is to identify the macroscopic value of the cortical bone permeability coefficient. A simple experimental method was designed in order to determine the permeability coefficient. Two bone samples taken from different ox femurs were filled with water, to place them under internal pressure. The measurements gave both the fluid flow through the lateral surfaces and the internal pressure. The originality of this work is the coupling between an experimental process and a structural computation performed with a finite element method. The mean cortical bone permeability coefficient identified was about k = 1.1 × 10 -13 m 2 . This value tends to confirm other values found in the literature, obtained by different methods and often at macroscopic scale. It confirms also the domination of vascular permeability (Haversian and Volkmann's canals).

Introduction

In order to take into account the complex system of bone microstructure, many authors have relied on poroelasticity theory [START_REF] Biot | Problème de la consolidation des matières argileuses sous une charge[END_REF][START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF][START_REF] Coussy | Mechanics of porous continua[END_REF]. In this work, we focus on the identification of a poroelastic parameter : the cortical bone permeability. Some values of permeability are given in [START_REF] Cowin | Bone poroelasticity[END_REF][START_REF] Li | Permeability of cortical bone of canine tibiae[END_REF][START_REF] Zhang | Estimates of the peak pressures in bone pore water[END_REF][START_REF] Smit | Estimation of the poroelastic parameters of cortical bone[END_REF][START_REF] Johnson | Fluid flow in bone in vitro[END_REF], arising from experimental data, analytic developments or numerical simulations. They have been estimated at different scales, leading to a wide dispersion of values. In this communication we propose an identification of permeability at macroscopic level, and the objective is to provide a procedure to identify the macroscopic permeability of cortical bones, without using local samples at various locations in the bones, but with a bone slice considered as a struture in its own. Coupling experimental measurements and 3D structural simulations enables identification of the average material characteristics. This procedure is validated on two samples of different ox femurs for which the results are compared with those obtained in the literature.

Specimen preparation, experimental device and results

The principle of the test bed is to prescribe a water pressure gradient on the sample and to measure the resulting fluid flux throughout the cortical part. Samples were taken from several ox femurs. Ox had been dead for three days when the experimental procedure began and had been preserved at a temperature of 4 • C. Slices of approximately 20 mm height were cut out, and the marrow was removed. Samples were obtained in the middle of the long axis. This part of femur is less vascularized, which seems to reduce the risk of blood coagulation.

Figure 1 shows the first sample after preparation. Because of the structure of a femur, samples were approximately ellipticals with major axis about 70 mm long for the first sample and 65 mm long for the second sample and minor axis about 40 mm long for the first and 43 mm long for the second. On the lateral surface all the nerves and cartilages were also removed, to obtain a "clean" sample. Finally the cut surfaces are polished in order to ensure a sufficient flatness, and to reduce the roughness. The design of the experimental device, see Figure 2, came from permeameters used for permeability measurements of soils. This was simplified to take advantage of the rigidity of the sample. The bone was gripped between two steel plates. Both were bored at their center to connect the internal cavity of the sample to a cylindrical capillary tube of 2 mm calibrated inner diameter d, and 3 m height. The watertightness was ensured by two polymeric joints between the extremal surfaces and the plates. The bottom hole was used for the initial fill-in and to clear the device at the end of the test. Water was injected to fill the sample and the capillary tube, and to pressurize the internal cavity of the sample. After one or two minutes, drops of water appeared on the lateral surface of the bone as shown in Figure 3.

Q = rP ( 1
)
where r is the hydraulic resistance. Therefore, dh dt = -r s ρgh and:

h h 0 = e -t τ
where τ = s rρg is a characteristic time scale of the structure (material plus geometry). A test is reported as the evolution of the logarithm of h h0 versus time t. We performed five series of tests on five bone samples. The five bone Unfortunately such a watertightness was only ensured for two samples. When the water level on the capillary tube was low, the device was refilled, without changing any conditions. A pause of several minutes was observed between each test. During the first moments of measurement, the liquid passing through the bone was quite colored. Once the water had drained the bone, the outgoing liquid became transparent. We can first observe differences between the five measurements as shown in Figure 4. The linearity of these evolutions validates the hydraulic law [START_REF] Biot | Problème de la consolidation des matières argileuses sous une charge[END_REF]. Since the slope evolves during the tests differently for each sample, the interpretation of this phenomenon is still unexplained; we chose here to consider only the mean slope 1 τ . It is identified as 1 τ = 2.26 × 10 -4 s -1 for the first sample and 1 τ = 1.43 × 10 -4 s -1 for the second one. Since this is not a material characteristic, the two sets of results on 4 are not directly comparable. A structural computation is required to identify permeability as a material characteristic.

Front view

Left view

Geometrical and physical model

To take into account the complex geometry of the samples, CAD numerical models were built, using measured contour sampling points and Bezier-type interpolation within CATIA TM software, see Figure 5 for the first sample.

Since the tests were quasi-static, we consider that the steady state is always obtained. Darcy's law reads where W is the opposite of Darcy velocity and K c is the hydraulic permeability coefficient. The boundary conditions are:

p = P on ∂ 1 Ω (2) p = 0 on ∂ 2 Ω (3) W • n = 0 on ∂ 3 Ω (4)
n being the unitary outwards vector, ∂ 1 Ω is the internal surface where the inner pressure is applied, ∂ 2 Ω is the external surface where the outer relative pressure is nil. On the extremal surfaces ∂ 3 Ω, the impervious plates prescribed a nil fluid flow.

Since we used a finite element method, a mesh of the sample was generated with the CAD software (30,000 linear tetrahedral elements were generated). A weak formulation of the problem is to find p ∈ P such that:

∀p * ∈ P 0 , Ω grad p • K c grad p * dΩ = 0
where P is the set of pore pressure fields p with a finite energy, and p| ∂1Ω = P, p| ∂2Ω = 0. P 0 is the associated vector space. After spatial discretization we obtain:

Hp = q
H is the permeability matrix, p is the column vector of nodal pressure, q is the column vector of generalised nodal fluid flux. If the degrees of freedom are split into internal ones, with subscript i, and degrees of freedom on ∂ 1 Ω, ∂ 2 Ω and ∂ 3 Ω denoted with a subscript e one gets,

H ii H ie H ei H ee p i p e = q i q e ( 5
)
with q i = 0 and p e = P u e where u e is a boolean vector whose unitary entries correspond to nodes on ∂ 1 Ω. The overall fluid flow is then given by:

Q = u T e q e ( 6
)
Algebraic manipulations lead to the relationship between permeability and total fluid flux is:

Q = u T e H *
ee u e P with the condensed permeability matrix

H * ee = H ee -H ei H -1 ii H ie .
H and H * ee depends lineary on K c and:

H = K c H 0
as well as:

H * ee = K c H * ,0 ee
The superscript 0 denotes operators obtained with a unitary permeability coefficient. The structural effect is stored in u T e H * ,0 ee u e and the material characteristic is given by:

K c = r u T
e H * ,0 ee u e Figure 6 shows the different permeability values identified during the five tests on both samples. The identified permeability coefficient is finally obtained

K c = 1.1 × 10 -13 m 4 N -1 s -1
With the water dynamic viscosity µ = 1 Pa.s, one gets the intrinsic permeability of the cortical bone k = µK c = 1.1 × 10 -13 m 2 . With such a procedure, only radial permeability is found. Other tests with different fluid flow orientations would lead to the full identification of a permeability tensor.

Discussion

The published studies to determine bone permeability mainly differ in terms of the scale at which experiments were done, and the nature of the bone. For heterogeneous structures, the local permeability is related to a particular location on the bone sample, and several experiments are required to get a representative value. On the other hand, we expect a global permeability measurement to allow the identification of the average value. Zhang et al [START_REF] Zhang | Estimates of the peak pressures in bone pore water[END_REF] calculated permeability analytically and obtained 6.36 × 10 -13 m 2 (at osteon level) and 1.47 × 10 -20 m 2 (associated to lacunae-canalicular porosity). Smit [START_REF] Smit | Estimation of the poroelastic parameters of cortical bone[END_REF] used a finite element model of an osteon and found K c = 2.2 × 10 -19 m 4 N -1 s -1 . So we consider that these permeability values are related to microscopic scales, and are below the scale of our study. In [START_REF] Li | Permeability of cortical bone of canine tibiae[END_REF] an experimental process on cortical bone canine tibiae collected bone plugs at different levels (proximal, middle and distal), and locations of the diaphysis (anterior, medial, lateral and posterior).

The mean value of intrinsic permeability was 5 × 10 -15 m 2 . Our permeability value falls between the value of the vascular permeability given by [START_REF] Zhang | Estimates of the peak pressures in bone pore water[END_REF] and the mean permeability found experimentally on dog tibiae in [START_REF] Li | Permeability of cortical bone of canine tibiae[END_REF]. The discrepancy between these values is probably due to differences in the bone structure (canine tibiae, bone of ox femur and slice location). Finally, bone characteristics are subjected to variability due to the biological process of material creation. Moreover, our value tends to experimentally confirm the theoretical value given by [START_REF] Zhang | Estimates of the peak pressures in bone pore water[END_REF] and seems to be close to vascular permeability, higher than lacunae-canalicular permeability. This result agrees with the domination of vascular permeability [START_REF] Zhang | Estimates of the peak pressures in bone pore water[END_REF].

Cancellous bone was not removed from our samples but (i) its thickness is smaller than that of the cortical phase, and (ii) its permeability coefficient is considerably higher than the permeability of cortical bone: cancellous bone permeability values in [START_REF] Kohles | Direct perfusion measurements of cancellous bone anisotropic permeability[END_REF] are between 4.65 × 10 -10 m 2 and 2.33 × 10 -10 m 2 . Consequently, this should not entail the identified value we obtained.

Conclusion

The aim of this work was to identify a global macroscopic value of cortical bone permeability. In order to take into account the geometry of the sample used and also its mechanical parameters, we associated a structural simulation with an experimental process. Our result agrees with the permeability values obtained by several authors but with a different method and at a different scale, mainly microscopics. Moreover it validates the domination of vascular permeability with respect to flow through porous media. This experimental process gives a measure of global radial permeability of a bone sample through the cortical thickness.

Figure 1 :

 1 Figure 1: Photography of the first bone sample approximatively elliptical with major axis about 70 mm long and minor axis about 40 mm long

Figure 2 :

 2 Figure 2: Experimental device. Water is injected at the bottom of a 2 mm diameter capillary tube

Figure 3 :

 3 Figure 3: Bone sample under internal pressure. One may notice on the left picture that drops appear on the lateral surface

Figure 4 :Figure 5 : 6 W

 456 Figure 4: Measurements of water height evolution h normalised by the initial height h 0 vs. time; for the first sample (top); for the second sample (bottom)

Figure 6 :

 6 Figure 6: Permeability values identified on first and second sample

Conflict of interest

All authors disclose any financial and personal relationships with other people or organisations that could inappropriately influence (bias) their work.


David Dureisseix

Bertrand Marcon

A partitioning strategy for the coupled hygromechanical analysis with application to wood structures of Cultural Heritage

Keywords: fluid-structure interaction, multiphysics, hygroscopic, anisotropy, thermomechanics, Jacobi, paintings, Mona Lisa

One of the main causes for damage of panel paintings of cultural heritage is due to the environmental micro-climatic variations. Wood science and numerical modeling may help to analyze and calibrate restoration interventions, to predict the behavior of the artwork, taking into account the individual painted panels and the conservation environment. In this article, a partitioning numerical strategy is proposed to predict by simulation the behavior of such artworks. It is based on a multiphysics partition focusing on the interactions between the different physics that are described in a thermodynamical framework. It is applied to the case study of Mona Lisa, modeled as a strongly coupled hygromechanical structure. The strategy is designed to couple two different modelings: a plate model for the mechanical behavior of the panel, and one-dimensional transverse diffusion for the moisture evolution. This is the post-

INTRODUCTION

Wood mechanical behavior is strongly dependent on moisture, and this hygromechanical coupled phenomena is itself dependent on temperature. Therefore, painted wood panels of cultural heritage (furniture, musical instruments, paintings...) are often sensitive to several loadings such as mechanical restraint, as well as humidity and temperature cycles due to the hygrothermal variations of the surrounding environment [START_REF] Mecklenburg | Structural response of painted wood surfaces to changes in ambient relative humidity[END_REF][START_REF] Camuffo | The microclimate inside the Pollaiolo and Botticelli rooms in the Uffizi Gallery, Florence[END_REF][START_REF] Vici | Mechanical response of wooden boards subjected to humidity step variations: Climatic chamber measurements and fitted mathematical models[END_REF]. To improve preventive conservation and to guide restoration acts on wooden artworks, virtual testing via numerical simulation is necessary to assess a risk analysis, in order to predict if the object will remain safe under various scenarios.

In this article, we focus on a numerical strategy to solve the coupled problem of interactions between the panel support and the pictural layer and we consider humidity variations only as a first step, with isothermal evolutions. Of course, the long term goal of such studies also requires interdisciplinary experiments, at least to nurture the model with information on the initial state of the artwork. For the application we use in this article, the case of Mona Lisa, we rely on the study reported in [START_REF] Mohen | Mona Lisa, inside the painting[END_REF]. In all of the following, the small perturbation assumption is used for the structural model.

More sophisticated models are available in the literature, especially for wood drying [START_REF] Perré | A heterogeneous wood drying computational model that accounts for material property variation across growth rings[END_REF][START_REF] Perré | The use of macroscopic equations to simulate heat and mass transfer in porous media[END_REF][START_REF] Svensson | Simulation of drying stresses in wood. Part 1: Comparison between one-and two-dimensional models[END_REF][START_REF] Whitaker | Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying[END_REF]. We choose a simpler one herein in order to focus on the numerical strategy to solve the hygromechanical coupled problem via a partitioning strategy. Nevertheless, to design such a strategy, we need for a consistent thermodynamics approach, especially for selecting the couples of dual variables involved in the model. Apart from the coupled physics, the second difficulty lies in the anisotropy arising from the wood underlying micro-structure.

A three-dimensional hygromechanical model is therefore stated within the framework of thermomechanics of materials. This study is afterwards dedicated to one-sided painted panels for which main moisture gradient is along the thickness of the panel. The panel itself is modeled as a plate. The final hygromechanical coupled problem is solved with a partitioning approach relying on a structural solver on one hand, and a diffusion code on the other hand. For this strongly coupled problem, an inner iteration loop is used to recover the coupling effects.

This article is organized as follows: Section 2 recalls the basic behavior and notations for hygromechanical problems, and Section 3 details the model we are interested in. The computational strategy is designed in Section 4 in the case of identical 3D meshes for the elastic and for the diffusion problems, while Section 5 couples two dedicated models and discretizations for the application to one-sided painted wood panels, and displays the numerical results. The influence of the discretization parameters on the accuracy of the solution and on the convergence of the partitioning scheme is studied.

w

HYGROMECHANICAL BEHAVIOR 2.1 Water in wood

Water in wood material can occur at different levels:

• chemical water that cannot be extracted without transforming the material constituents;

• bound (or hygroscopic) water in cell walls;

• water vapor in lumens;

• free (or capillary) water in liquid form in the voids of the wood.

The moisture content w is the mass of water contained in the last three forms, in an elementary representative volume, expressed as a fraction of its oven-dry mass. Using apparent water specific density ρ w and specific density of the oven-dry wood ρ o , one gets

w = ρ w ρ o (1)
Wood panels of cultural heritage are usually loaded within hygroscopic regime, for which water mainly occurs as bound water and water vapor. In such cases, the moisture content is below the fiber saturation point (FSP) for which it is approximately w FSP = 30 % [START_REF] Siau | Transport processes in wood[END_REF][START_REF] Skaar | Wood-water relations[END_REF].

Variations of humidity are supposed to be sufficiently slow for bound water to be in local equilibrium with the relative humidity of the air, defined as H = p/p s , where p is the partial water vapor pressure in the air, and p s is the saturated vapor pressure (therefore, one gets 0 ≤ H ≤ 1).

In this study, the material behavior will be described at macroscopic scale, and only diffusion is considered as a water transfer mechanism (no pressure gradients are involved that can lead to a bulk viscous flow). Table 1: Elastic parameters at reference state (after [START_REF] Guitard | Mécanique du matériau bois et composites[END_REF][START_REF] Dureisseix | Mechanical modeling of the activity of the flexible frame, chapter II.6[END_REF]); sorption and shrinkage parameters; diffusivity parameters (after [START_REF] Olek | Inverse analysis of the transient bound water diffusion in wood[END_REF])

Elastic behavior

Concerning the elastic behavior, the stress field σ and the strain field ε are linked via the Hooke tensor. The elastic behavior is orthotropic; moreover in a wood piece, due to the radial growth of the tree, the directions of anisotropy are heterogeneous (roughly with a cylindrical symmetry). These directions are usually denoted with L, R and T for the longitudinal, radial and tangential directions, respectively. The elastic part of the material behavior is represented with Hooke tensor C. The compliant tensor S is classically expressed in the orthotropic base (L,R,T) with Voigt notations:

S = C -1 =         1/E L -ν LR /E L -ν LT /E L 0 0 0 -ν RL /E R 1/E R -ν RT /E R 0 0 0 -ν T L /E T -ν T R /E T 1/E T 0 0 0 0 0 0 1/G RT 0 0 0 0 0 0 1/G T L 0 0 0 0 0 0 1/G LR         (L,R,T) (2)
with symmetry, positiveness and definiteness conditions. E α are Young modulus in direction α = L,R,T, ν αβ are Poisson coefficients, G αβ are the shear modulus.

For wooden structures of cultural heritage, it is often impossible to extract a specimen of wood and perform a tensile test; only measurements without contact are allowed. Therefore, estimations of material characteristic parameters are useful. Concerning the case of poplar wood (Populus sp.), which is under concern herein, [START_REF] Guitard | Mécanique du matériau bois et composites[END_REF] provides estimation of elastic parameters for the reference state, as functions of density. Density can be obtained in a non destructive manner. We consider here a specific density ρ = 0.45 ± 0.025 g/cm 3 [START_REF] Uzielli | Physical and mechanical characterization of the support, chapter II.4[END_REF], and the reference state w = 12 %, that leads to the estimations of Table 1. They proved to reproduce experimental observations [START_REF] Dureisseix | Mechanical modeling of the activity of the flexible frame, chapter II.6[END_REF].

Sorption curve, swelling and shrinkage

The basic test for the hygromechanical behavior consists in placing an elementary representative volume (or a specimen with homogeneous properties) in a humidity-controlled and temperature-regulated air. Increasing slowly the humidity H to get an equilibrium with the moisture content w in all the body of the specimen allows to get the value of this moisture content by monitoring the weight of the specimen. This isothermal sorption curve is usually obtained under null stress. If not, an additional mechanism known as mechano-sorption is involved [START_REF] Ranta-Maunus | Impact of mechano-sorptive creep to the long-term strength of timber[END_REF][START_REF] Salin | Numerical prediction of checking during timber drying, a new mechanosorptive creep model[END_REF][START_REF] Navi | Some new aspects of wood mechanosorptive effects[END_REF]; this mechanism is not taken into account herein. Usually an hysteresis is obtained between adsorption and desorption, due to chemical interactions [START_REF] Ma | Prediction of vapor-moisture equilibriums for a wood-moisture system using a modified UNIQUAC model[END_REF][START_REF] Frandsen | A hysteresis model suitable for numerical simulation of moisture content in wood[END_REF]. We will not take this effect into account here, and we will consider the average curve of Figure 1. This Figure plots data obtained from different sources, and the Guggenheim, Anderson and Boer-Dent (GAB) model is identified to obtain an average curve, from data in [START_REF] Peralta | Sorption of moisture by wood within a limited range of relative humidities[END_REF]:

w = w m CKH (1 -KH)(1 -KH + CKH) (3)
The parameter values are referenced in Table 1.

To derive a constitutive model, we do need to define a reference state, on a reference configuration, from which the displacement of the structure U is measured (i.e. ε ref = 0). We choose here the standard state at a temperature of 25 • C, and a relative humidity H ref = 65 % (corresponding approximately to w ref = 12 %, in the isothermal sorption curve of Figure 1). Moreover, we assume that for the same Figure 1: Isothermal sorption-desorption curve (after [START_REF] Peralta | Sorption of moisture by wood within a limited range of relative humidities[END_REF][START_REF] Bratasz | Sorption of moisture and dimensional change of wood species used in historic objects[END_REF][START_REF] Garcìa Esteban | Comparison of the hygroscopic behaviour of 205-year-old and recently cut juvenile wood from Pinus sylvestris L[END_REF]23]) reference configuration, the stress is null: σ ref = 0. Concerning the sorption test, the relationship between moisture content and relative humidity (3) can be written as

∆w = f (H/H ref ) (4)
where ∆w = w -w ref .

During the same test, measurements of dimension variations of the specimen may characterize the swelling and shrinkage effect, that mainly appears during the hygroscopic regime. Since this effect is also anisotropic, all the principal components of the corresponding free strain ε have to be measured. A linear model is usually satisfactory, and the fitting leads to ε = A∆w (when σ = 0) [START_REF] Perré | A heterogeneous wood drying computational model that accounts for material property variation across growth rings[END_REF] where, in the orthotropic basis,

A =   a L 0 0 0 a R 0 0 0 a T   (6)
with constant adimensional coefficients a α . To estimate these coefficients, the approach in [START_REF] Préziosa | Contraintes internes de séchage dans le matériau bois : le tenseur des coefficients de contrainte de séchage comme caractéristique isotrope de la matière ligneuse[END_REF] is used: for a large set of wood species, A = Sβ1 where 1 is the identity tensor and β = K β ρ 2 o with K β ≈ 2.15 GPa/(g/cm 3 ) 2 . This leads to the values reported in Table 1 (for instance, with ∆w = 10 % the free tangential strain is a T ∆w = 4.8 %).

The shrinkage tensor A is indeed a coupling between mechanical behavior and hygrometry.

REFERENCE PROBLEM

Admissibility conditions

Concerning the structural problem on the domain Ω, a given displacement field U d is prescribed on a first part ∂ 1 Ω of the boundary ∂Ω. A surface force density F d on the complementary part ∂ 2 Ω, and a bulk force density f d are also prescribed. The conservation of momentum gives the static admissibility conditions on the stress field σ:

div σ + f d = 0 in Ω, σn = F d on ∂ 2 Ω (7)
n is the outward normal unitary vector at the boundary. The strain field ε should be compatible with a kinematically admissible displacement field U :

ε = (Grad U ) sym in Ω, U ∈ U, U = U d on ∂ 1 Ω ( 8
)
where U is the set of finite energy displacement field, typically:

U = H 1 (Ω).
The conservation of fluid mass leads to

∂ρ w ∂t = -div J (9)
where J is the water mass flux. With the expression (1) of the moisture content, this reads

ρ o ∂w ∂t = -div J (10)
The boundary conditions corresponding to the fluid quantities may be of several types depending on the problem. They are usually of one of the forms (Dirichlet, Neumann, or Robin-type):

w = w d , J • n = j d , ρ o k(w -w d ) = J • n (11)
where w d is a prescribed moisture content, j d is a given normal flux, and k is a surface emissivity (the inverse of an hydraulic resistance), where n is the normal to the boundary.

Strongly coupled constitutive model

The Young modulus are known to be dependent on the moisture content, i.e. C = C(∆w).

On the other hand, the diffusion of water in hygroscopic wood depends on diffusion of bound water through cell walls and of water vapor in the lumens, both being in equilibrium [START_REF] Baronas | Modelling of moisture movement in wood during outdoor storage[END_REF][START_REF] Siau | Transport processes in wood[END_REF][START_REF] Avramidis | Experiments in nonisothermal diffusion of moisture in wood[END_REF][START_REF] Cloutier | A model of moisture movement in wood based on water potential and the determination of the effective water conductivity[END_REF]. At the macroscopic scale, the overall effect is often modeled as a global diffusion with a unique diffusivity tensor D [START_REF] Barkas | The swelling of wood under stress[END_REF][START_REF] Merakeb | Couplage hygromcanique dans le processus de diffusion dans le bois[END_REF], and a Fick equation, for which

J = -ρ o D grad ∆w (12)
To keep the constitutive model consistent, a thermodynamics framework is compulsory. The state of the material in each elementary representative volume is described by a certain number of variables. For instance, for the structural part, the strain ε accounts for the deformation phenomena. The fluid part can also be described by a variable, the amount of water ρ o ∆w. We restrict ourselves here to these two variables, as sufficient ones to describe the internal material state.

To derive the state laws of the constitutive model, the Gibbs free energy is expressed as a function of the state of the material: Ψ(ε, ρ o ∆w). Following this description, the state laws can be derived as

σ r := ∂Ψ ∂ε (13)
A w := ∂Ψ ∂ρ o ∆w [START_REF] Olek | Inverse analysis of the transient bound water diffusion in wood[END_REF] where σ r is defined as a reversible stress, and A w is the thermodynamic dual variable of ρ o ∆w. The dissipation is derived from the Clausius-Duhem inequality, see [START_REF] Coussy | Poromechanics[END_REF] for instance, as

d = σ : ε + ∆µ ρ o ∆ ẇ -Ψ -J • Z 0 ( 15
)
where µ is the massic chemical potential of water [START_REF] Siau | Transport processes in wood[END_REF][START_REF] Siau | Nonisothermal diffusion model based on irreversible thermodynamics[END_REF], and Z is its spatial gradient: Z = grad ∆µ. At local equilibrium, µ is also the chemical potential of the water vapor, that can be linked to the relative humidity H, with the assumption of perfect gas,

µ = µ 0 + RT M ln H ( 16
)
µ 0 is the value of the massic chemical potential for a saturated vapor, R = 8.31 J K -1 mol -1 is the perfect gas constant, T is the absolute temperature, here T = T ref = 298 K, and M = 18 g/mol is the molar mass of water. The corresponding value at the reference state is µ ref , and

∆µ = µ -µ ref = RT M ln H H ref (17)
Using the previous state laws, the dissipation reads

d = (σ -σ r ) : ε + (∆µ -A w )ρ o ∆ ẇ -J • Z 0 ( 18
)
This dissipation is assumed to be expressed from the gradient of a potential, here chosen as

Φ(Z) = 1 2 Z • LZ ( 19
)
where L is a diffusion tensor for chemical potential (homogeneous to kg.m -3 .s 2 ). This choice is made to recover a Fick-like diffusion behavior. Indeed, with d = (∂Φ/∂Z) • Z, one gets the evolution laws:

σ = σ r (20) ∆µ = A w (21)
-J = LZ (22)
The fact that the moisture flux can be driven by the moisture gradient [START_REF] Uzielli | Physical and mechanical characterization of the support, chapter II.4[END_REF] or the chemical potential gradient [START_REF] Garcìa Esteban | Comparison of the hygroscopic behaviour of 205-year-old and recently cut juvenile wood from Pinus sylvestris L[END_REF] depends on the modeling, and on the fact that the corresponding diffusivity tensors may depend on the current state. For instance, with the state description with the couple (ε, ρ o ∆w), [START_REF] Garcìa Esteban | Comparison of the hygroscopic behaviour of 205-year-old and recently cut juvenile wood from Pinus sylvestris L[END_REF] can be written as

J = -L grad A w = -L ∂A w ∂∆w grad ∆w -L ∂A w ∂ε : Grad ε ( 23
)
where the last term (dependence on the stress strain, similar to the Gorsky effect with the stress gradient [START_REF] Larché | The interactions of composition and stress in cristalline solids[END_REF][START_REF] Glicksman | Diffusion in Solids: Field Theory, Solid-state Principles, and Applications[END_REF]) is usually neglected for the present application. Then one gets the relationship:

ρ o D = L ∂A w ∂∆w (24)
with emphasize the possible dependence of the diffusivity on the current state.

We now have to precise the constitutive model. First, we propose the modeling choice for the Gibbs potential:

Ψ(ε, ρ o ∆w) = 1 2 ε : C(∆w)ε -B 0 (∆w) : C(∆w)ε + Ψ 0 (∆w) (25)
for which the state laws are

σ r = ∂Ψ ∂ε = C(ε -B 0 ) (26)
A w = ∂Ψ ∂ρ o ∆w = 1 ρ o 1 2 ε : C ε -B 0 : Cε -B 0 : C ε + Ψ 0 ( 27
)
B 0 is homogeneous to a strain tensor, coupling elasticity and moisture, corresponding to the swelling / shrinkage phenomenon; Ψ 0 is a potential to be identified. Indeed, when considering an evolution at null stress, ε | σ=0 = B 0 (∆w) and the fitting [START_REF] Perré | A heterogeneous wood drying computational model that accounts for material property variation across growth rings[END_REF] gives

B 0 (∆w) = A∆w (28)
With the isothermal sorption curve at null stress (4) and the expression [START_REF] Navi | Some new aspects of wood mechanosorptive effects[END_REF], one gets

Ψ 0 (∆w) = ρ o RT M ln[f -1 (∆w)] + 1 2 ∆w(A : C A)∆w + (A : CA)∆w (29)
Note that the inverse function f -1 is required. With the GAB model used herein, this inversion can be done analytically (at least on the range -1 (C-1)K < H < 1 K ; to be able to use the logarithm, one may restrict himself to 0 < H < 1 K , which corresponds to w > 0). If not, it is still possible to represent the sorption curve as a tabulated function, given by experimental values, that can be inverted by switching entries of the table, with an interpolation between experimental points and a suitable management of out-of-range values. The first term is indeed the variation of chemical potential during the sorption test, see Figure 2.

Choice of the coupling parameters

To express the dependency of Hooke operator with respect to the moisture content, a practical approach consists in considering the components of the 6 × 6 compliant matrix with Voigt notations S ij , i = 1, ...6, j = 1, ..., 6. Each of them is usually linked to the moisture content w as follows: [START_REF] Gerhards | Effect of moisture content and temperature on the mechanical properties of wood[END_REF][START_REF] Kollman | Principles of wood sciences and technology[END_REF][START_REF] Jodin | Le bois, matériau d'ingénierie[END_REF] are used here: with the same notation as in Section 2.2, these parameters are stored in the matrix b:

1 S ij = 1 S ref ij (1 -b ij ∆w) (30)
b =         1.0 1.2 1.2 0 0 0 1.2 3.0 3.0 0 0 0 1.2 3.0 3.0 0 0 0 0 0 0 2 × 2.0 0 0 0 0 0 0 2 × 2.0 0 0 0 0 0 0 2 × 2.0         (L,R,T) (31)
This model is consistent for a limited range of moisture content values, especially to maintain the positive definiteness of the Hooke tensor For the considered values, one can check that this is the case up to ∆w ≈ 25 %, i.e. far within the range of the variation of moisture content for our application (see Figure 1). It is nevertheless checked in the code that intermediate solutions during the iterations do not exceed this value. As a consequence, C = S -1 is a non linear function of the moisture content: It decreases when the moisture content increases; therefore the following notation is useful for subsequent developments:

C(∆w) = C ref -C(∆w) (32)
where

C ref = C(∆w = 0).
The diffusivity dependence on the current state is usually determine with a phenomenological identification. Within the orthotropic basis,

D ref =   D L 0 0 0 D R 0 0 0 D T   (33)
For the dependence on moisture, under null stress, we choose to reuse the values obtained for European beech in [START_REF] Olek | Inverse analysis of the transient bound water diffusion in wood[END_REF][START_REF] Olek | Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood[END_REF], for which D α = D refα exp(c α ∆w), α = L,R,T, and the diffusivity values are recalled in Table 1. For the dependence on stress, we rely again on a phenomenological identification, though it is difficult to perform here due to the neglected mechano-sorptive effects, that nevertheless may occur when applying a prescribed stress on the specimen. According to the relationship between different expressions of the diffusion driving gradient that can be obtained (83), we use herein the following model:

D α = exp(c α ∆w) 1 + 1 σ0 A : σ D refα (34)
The value of σ 0 can be identified from experiments as those in [START_REF] Merakeb | Couplage hygromcanique dans le processus de diffusion dans le bois[END_REF] (for Douglas Fir), and we use the corresponding value of σ 0 = 1.9 MPa. Note that, with an iterative scheme, the diffusivity may have to be evaluated for non physical values of moisture content (i.e. values for non converged quantities during iterations). In order to increase the robustness of the approach, the exponential form is cutted out over the physical window of moisture content values, i.e. ∆w is replaced in the previous expression [START_REF] Gerhards | Effect of moisture content and temperature on the mechanical properties of wood[END_REF] with min[max(∆w, -w ref ), 1 -w ref ], (-w ref ) and (1 -w ref ) being the extrema for the allowable values of moisture content ∆w. This lead to use a dependence on the moisture content of the form described in Figure 3.

PARTITIONING STRATEGY FOR 3D PROBLEM

The main interest in using a partitioning scheme rather than a monolithic solver when dealing with multiphysics problems lies in (see [START_REF] Park | Stabilization of partitioned solution procedure for pore fluid-soil interaction analysis[END_REF][START_REF] Turska | On convergence conditions of partitioned solution procedures for consolidation problems[END_REF][START_REF] Felippa | Partitioned analysis of coupled mechanical systems[END_REF][START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems. Part II: Energy transfer analysis and three-dimensional applications[END_REF][START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Dureisseix | A LATIN computational strategy for multiphysics problems: Application to poroelasticity[END_REF] for instance):

• the possibility to couple different numerical codes (for the different physics);

• the increase in flexibility (useful when adding a new physics, or for code maintenance and evolution).

The main difficulties concern the case of strong coupling for which

• iterations are mandatory to recover the coupling effects;

• the robustness is not ensured, which may lead to a lack of convergence.

Herein, we focus on a class of partitioning schemes that promote the interactions between subsystems, i.e. the local coupled equations between the different physics, by emphasizing on material fields that participate to the interactions. Within this class of solvers, one may find the Large Time Increment (LATIN) approach [START_REF] Ladevèze | Nonlinear computational structural mechanics -New approaches and non-incremental methods of calculation[END_REF][START_REF] Néron | A computational strategy for thermo-poroelastic structures with a timespace interface coupling[END_REF][START_REF] Dureisseix | Two examples of partitioning approaches for multiscale and multiphysics coupled problems[END_REF], and for coupling different bodies with contact and friction, the Non Smooth Contact Dynamics (NSCD) approach [START_REF] Moreau | Numerical aspects of sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF].

These approaches lead to isolate the treatment of the interactions from the behavior of the subsystems themselves for which adequate techniques can be used in a modular fashion. Indeed, they deal with iterations involving a local coupled resolution when expressing the coupled constitutive equations on one hand, and admissibility conditions on the other hand, as does the proposed approach.

Splitting the constitutive relations

To focus on coupled and non-linear equations, constitutive relations are partitioned according to the following splitting:

σ = C ref ε + σ ( 35
)
where the pair ( σ, ε) is linked with ∆w:

σ = -C(∆w)ε -C(∆w)A∆w (36)
On the other hand, ∆µ = c∆w + μ [START_REF] Olek | Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood[END_REF] with c = 1 ρ o A : C ref A and the pair (μ, ∆w) is linked with ε:

μ = µ 1 (ε, ∆w) + µ 2 (∆w) (38)
where µ 2 (∆w) collects nonlinear terms, and µ 1 (ε, ∆w) the coupling terms:

µ 1 (ε, ∆w) = 1 ρ o [- 1 2 ε : C ε -A : Cε + A : C ε∆w] (39)
µ 2 (∆w) = RT M ln[f -1 (∆w)] - 1 ρ o [(A : CA)∆w + 1 2 (A : C A)(∆w) 2 ] ( 40
)
Finally for the diffusion part, using a constant diffusivity D 0 as a parameter for the approach,

J = -ρ o D 0 F + J ( 41
)
where F = grad ∆w. The pair ( J, F ) is linked with σ and ∆w:

J = -ρ o [D(σ, ∆w) -D 0 ]F (42)

Admissibility closure equations

The admissibility conditions arise from the conservation principles. They are linear when a suited choice of variables is used, and they will be used here to close the problem by adding equations linking the unknown couples ( σ, ε) on one side, independently of (μ, ∆w) and ( J, F ) on the other side, while being kept linear (but global on the studied spatial domain).

Solid part.

This structural problem is close to an elastic analysis, modeled with elliptic partial differential equations. Therefore, finite elements are a suited discretization tool, and they will be used in the following [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. With a displacement-oriented approach, the displacement field is described with a finite number of values, defined at the nodes of the finite element mesh, while strain and stress fields are evaluated at integration points within each element.

The conservation of momentum is

B T σ σ = f d ( 43
)
where σ is here the column vector of values of the stress tensor field σ, defined at integration points of the solid mesh, and f d are the generalized external forces. The operator B σ is classically defined as

∀u, Ω ε(u) : σ dΩ = u T B T σ σ
where u is the column vector of nodal displacement u. On the other hand, the strain column vector is

ε = B ε u (44)
With the previous constitutive-like relation [START_REF] Kollman | Principles of wood sciences and technology[END_REF], the previous admissibility equations give

[B T σ C ref B ε ] K u + B T σ σ = f d ( 45
)
where K is the stiffness matrix. With appropriate boundary conditions on the displacement field, this problem is regular, and the strain is

ε = B ε K -1 f d -X S σ (46)
with

X S = B ε K -1 B T σ .
Note that, from a computational point of view, X S should not be assembled explicitly, due to the large cost of producing the inverse of the stiffness matrix; on the other hand, its product with a stress requires a global elastic solve.

Fluid part.

This fluid problem is a transient diffusion problem, modeled with parabolic partial differential equations. Finite elements (unless stabilized [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origin of stabilized methods[END_REF][START_REF] Masud | A stabilized mixed finite element method for Darcy flow[END_REF]) are less suited to such a problem, but they will nevertheless be used in the following. Using more dedicated discretization tools, such as finite differences, finite volumes or even discrete elements and discontinuous Galerkin approach, would lead to different discrete matrices, but the coupling terms would be similar as those used herein. Coupling different discretization schemes is not under the scope of this article, but is currently under development. Here, the discretization of the moisture content field lead to select its values at the nodes of the mesh as primary unknowns. The integration points are the support of several dual fields: the moisture gradient associated to the water mass flux, and the chemical potential associated to the moisture content itself. Concerning time evolutions, a classical incremental scheme is used.

The conservation of mass for the fluid phase ( 10) is written at the discretized level as

-B T q ρ o B w ∆ ẇ + B T J J = g d (47)
The column vector of nodal moisture content is ∆w, and the vector of moisture gradient values is F = B F ∆w. These last values are expressed at the integration points of the fluid mesh, as for the vector of water mass flux J. By a simple extrapolation with finite element shape functions, the moisture content itself can be expressed at the same integration points: W = B w ∆w. Finally, the operators B q and B J are defined as ∀∆w,

Ω ∆w q dΩ = (∆w) T B T q q and Ω grad(∆w) • J dΩ = (∆w) T B T J J
With [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF] and after algebraic manipulations, one gets

[B T q ρ o B w ] M ∆ ẇ + [B T J ρ o D 0 B F ] H ∆w -B T J J = -g d (48)
M is a 'capacity'-like matrix, and H is the diffusivity matrix. A time integration scheme is needed (here the θ-method):

∆ ẇ = 1 θδt (∆w -∆w i ) -1-θ θ ∆ ẇi ,
where δt is the time step and θ ∈ [0, 1] is the parameter of the scheme, ∆w i and ∆ ẇi being the values at the previous time step. Therefore,

1 θδt M + H M ∆w = -g d - 1 θδt M [∆w i + (1 -θ)δt∆ ẇi ] gd +B T J J (49)
M is sometimes referred to as a 'dynamical' matrix. It is regular, and one gets

F = -B F M -1 gd + X F J ( 50
)
with

X F = B F M -1 B T J .
As previously, X F is not explicitly assembled, and its product with a flux requires the solve of a global diffusion problem.

A non linear Jacobi-like iterative partitioned scheme

At each time step of the incremental scheme, the reference problem consists in finding the fields ( σ, μ, J; ε, ∆w, F ) satisfying:

• the linear global solid structure behavior [START_REF] Dureisseix | Two examples of partitioning approaches for multiscale and multiphysics coupled problems[END_REF],

• the linear global fluid behavior ( 50) and (37),

• the nonlinear coupled local constitutive behaviors [START_REF] Jodin | Le bois, matériau d'ingénierie[END_REF], [START_REF] Park | Stabilization of partitioned solution procedure for pore fluid-soil interaction analysis[END_REF], [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF].

To derive a partitioning scheme, we propose to split the global operators into X S = X D S +(X S -X D S ), and

X F = X D F + (X F -X D F ). The operators X D
S and X D F are supposed to play the role of search directions, and to keep the coupled non-linearities local in space variables: they should be decoupled from each integration point to the other. They can be viewed as 'diagonals' of initial operators, or more simply, one can choose

X D S = C -1 ref and X D F = 1 ρo D -1 0 .
Concerning the parameter D 0 , one can choose the value of diffusivity for null stress and ∆w = 0, i.e. D ref .

An iterative procedure consists in updating an iterate ( σ, J) (k-1) following ( 50), [START_REF] Dureisseix | Two examples of partitioning approaches for multiscale and multiphysics coupled problems[END_REF], such that:

ε (k) + X D S σ(k) = B ε K -1 f d -(X S -X D S )σ (k-1) (51)
F (k) -X D F J(k) = -B F M -1 gd + (X F -X D F ) J(k-1) (52)
together with the constitutive relations [START_REF] Jodin | Le bois, matériau d'ingénierie[END_REF], [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF]. No admissibility condition involve the pair (μ, ∆w): Once ε and ∆w are known, ∆µ can be evaluated with [START_REF] Olek | Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood[END_REF] and [START_REF] Park | Stabilization of partitioned solution procedure for pore fluid-soil interaction analysis[END_REF]. This leads to the following (coupled but local at each integration point) system:

[(X D S ) -1 -C(W (k) )]ε (k) = σ d + C(W (k) )AW (k) (53) [(X D F ) -1 + ρ o (D(σ (k) , W (k) ) -D 0 )]F (k) = J d (54)
σ (k) = C(W (k) )(ε (k) + AW (k) ) (55) μ(k) = µ 1 (ε (k) , W (k) ) + µ 2 (W (k) ) (56)
where

σ d = (X D S ) -1 [B ε K -1 f d -(X S -X D S )σ (k-1) ] (57)
J d = (X D F ) -1 [-B F M -1 gd + (X F -X D F ) J(k-1) ] (58)
W (k) = B w ∆w (k-1) (59)
are constant during an iteration. Note that their evaluations require to solve two decoupled global and linear equations: one for an elastic structure problem, one for a linear diffusion problem. This problem is easily solved with a given (σ, ∆w, J) (k-1) : first, (σ d , J d , W (k) ) are computed with (57), ( 58), [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]; second, ( 53) is solved to get ε (k) ; third, ( 54), ( 55) are used to get F (k) , and finally, if needed, μ(k) is obtained from [START_REF] Verpeaux | CASTEM 2000 : une approche moderne du calcul des structures[END_REF]. Algorithm 1 gives the overall procedure. Note that, apart from the incremental scheme, this algorithm requires independent solves for each physics (a classical elastic problem and a classical diffusion problem), and local computations to recover the coupled effects. If the previous parameters X D S and X D F are selected as before, the global problems to be solved have constant left hand sides, and are factorized only once during the initialization phase.

Note that due to a strong coupling between fluid and structure, there is a need for an iterative scheme between the two physics for each time step, contrary to staggered schemes [START_REF] Felippa | Staggered transient analysis procedures for coupled mechanical systems: Formulation[END_REF][START_REF] Zienkiewicz | Unconditionally stable staggered solution procedure for soil-pore fluid interaction problems[END_REF][START_REF] Turska | On consistency, stability and convergence of staggered solution procedures[END_REF][START_REF] Schrefler | Standard staggered and staggered newton schemes in thermo-hydro-mechanical problems[END_REF][START_REF] Dureisseix | A LATIN computational strategy for multiphysics problems: Application to poroelasticity[END_REF]. At convergence for each time step, the produced solution is the same as for a full coupled solution scheme. The stability of the proposed partitioning approach is therefore entirely inherited from the time integration scheme that is used, and is not degraded.

For the termination criteria, apart from a classical stationary criteria, we proposed to use the following error estimators: the error indicators η S for the elastic (solid) part, and η F for the diffusion (fluid) part:

η (k) F = F (k) -F (k) F F (k) + F (k) F and η (k) S = ε (k) -ε(k) S ε (k) + ε(k) S (60)
with the energy norms:

• S = 1 2 Ω • C ref • dΩ 1/2 and • F = 1 2 Ω • ρ o D 0 • dΩ 1/2 (61)
Indeed, if fluxes F (k) and F (k) are equal ( F (k) = F (k) = F ), as well as strain fields ε(k) and

ε (k) (ε (k) = ε (k) = ε)
, with the local coupled problem solving, the fields satisfying the constitutive behaviors, J(k) and σ(k) are computed as functions of F and ε. Then, with the global solid problem, the new flux satisfying the admissibility conditions with σ(k) is again F . The global fluid problem ensures as well that ε is admissible with σ(k) . Therefore the set (F , J(k) ; ε, σ(k) ) satisfies both the constitutive relations and the admissibility conditions. If the reference problem is well posed, it is therefore the solution.

The convergence is declared when η S < S and η F < F with the user-given thresholds S and F .

Test case and performance comparisons

In this Section, the proposed strategy is compared to a monolithic scheme with a direct coupled solidfluid solver. In the following simplified case, we implement a classical modified Newton resolution in a standard 3D finite element code, namely Cast3M [START_REF] Verpeaux | CASTEM 2000 : une approche moderne du calcul des structures[END_REF]. For testing purposes, we simplify the geometry with an homogeneous 3D block of size 531 mm × 787 mm × 13 mm, and an orthotropic material with a local basis defined with a uniform angle α, Figure 4. We also consider simplified couplings: (i) a constant Set ∆w i ← 0, ∆ ẇi ← 0 3:

Set σi ← 0 4:
Set Ji ← 0 5: {Loop on time steps} 6: for i = 2, ..., n do 7:

Compute right-hand sides f di and

8: ∆ wi-1 ← ∆w i-1 + (1 -θ)δt∆ ẇi-1 9: gdi ← g di -1 θδt M ∆ wi-1 10:
Initial guess (k ← 1)

11:

Initiate σ(k) ← σi-1

12:

Initiate J(k) ← Ji-1

13:

{Loop on inner iterations} 14:

repeat 15: k ← k + 1 16:
Global solid problem 17:

Solve

Ku (k) = f di -B T σ σ(k-1) to get u (k) 18: Compute ε(k) ← B ε u (k) 19: Compute σ d ← (X D S ) -1 ε(k) + σ(k-1) 20:
Global fluid problem

21:

Solve M ∆w (k) = -g di + B T J J(k-1) to get ∆w (k) 22:

Compute W (k) ← B w ∆w (k)

23:

Compute F (k) ← B F ∆w (k)

24:

Compute

J d ← (X D F ) -1 F (k) -J(k-1)
25:

Local coupled problem (at each integration point)

26:

Set σd ← σ d + C(W (k) )AW (k) , and

27: Solve [(X D S ) -1 -C(W (k) )]ε (k) = σd to get ε (k) 28: Update σ(k) ← -(X D S ) -1 (ε (k) -ε(k) ) + σ(k-1) 29: Set σ ← C ref ε (k) + σ(k) and D ← D( σ, W (k) ) and 30: Solve [ρ o ( D -D 0 ) + (X D F ) -1 ]F (k) = J d to get F (k) 31: Update J(k) ← (X D F ) -1 (F (k) -F (k) ) + J(k-1) 32:
Compute η

S and η

(k) F 33:
{Convergence criteria}

34: if (η (k)
S < S and η

(k) F < F ) then 35: Set ∆w i ← ∆w (k) , ∆ ẇi ← 1 θδt (∆w i -∆ wi-1 )

36:

Set σi ← σ(k)

37:

Set Ji ← J(k)

with N = B T σ C ref A, and (ii) for not a too strong coupling, we linearized the dependence of the diffusion coefficients on the stress, and we neglect the dependence on the current moisture content, as:

ρ o D(σ) = 1 - 1 σ 0 A : σ ρ o D ref (63)
With a modified Newton scheme as a monolithic solver, the tangent operator has to be computed at the beginning of each time step, but is kept constant for the inner loops. Due to the norm in [START_REF] Engelund | The influence of moisture content on the water vapour resistance of surface coated spruce[END_REF], this expression is not differentiable at σ = 0, and so is the left-hand-side contribution H(σ) = B T J ρ o DB F . We therefore use a generalized Newton approach [START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF][START_REF] Alart | Méthode de newton généralisée en mécanique du contact[END_REF] for systems that are piecewise differentiable, for which

∂ σ H : σ = -B T J 1 σ 0 A : σ ρ o D ref B F (64)
The upper surface is prescribed to a null moisture flux (impervious boundary condition), while the remaining part of the boundary is subjected to a prescribed moisture content ∆w The comparisons are performed for different values of σ 0 (expected to change the strength of the coupling). The spatial discretization uses 1152 8-node hexahedral elements. The termination criteria for the partitioned approach uses S = F = = 10 -6 , and for the Newton solver, for both solid and fluid, the relative residual should not exceeds = 10 -6 as well. Both the monolithic and the partitioned approach lead to the same solution (up to the convergence precision ).

d (t) = f (H(t)/H ref ) with H(t) = H ref -0.
In each case, the cost is expected to be driven by the linear global solves. The linear problem to solve at each Newton iteration is half decoupled (the fluid problem can be solved first, then the solid problem can be solved in a second hand). This is therefore optimistic for the Newton approach, since then, the cost is equivalent for both the two approaches, at each global solve. Nevertheless, the tangent linear operator is updated at each time step for the modified Newton scheme (that requires as many factorizations as time steps, and as many forward-backward solves as inner Newton iterations), but the left-hand-side is constant for the partitioned approach (that requires only one factorization for each physics, and as many forward-backward solves as inner fixed-point iterations).

Table 2 compares the number of inner iterations for both approaches with respect to the value of the coupling strength σ 0 (decreasing its value leads to a stronger coupling). In the convergence domain, both approaches require roughly the same number of inner iterations. Since left-hand-sides for the partitioned approaches are factorized only once, this approach is therefore the most efficient. There is indeed a domain where the fixed-point like approach failed to converge when the coupled non-linearities become larger, but this convergence failure appears at the same time as for the Newton solver.

MODEL COUPLING

Due to the nature of the application, a full 3D model is not required. The aim of this Section is to adapt the previous approach to a particular case of plate model for elasticity and 1D model for diffusion. Therefore, this Section emphasizes the modularity of the proposed approach, since the only modification concerns the solver black-box. The same modularity feature could be used when different finite element meshes are used for each physics.

Structural model

Painted panels of cultural heritage often exhibit a small thickness to in-plane length ratio. for instance, in the case of Mona Lisa, the thickness e of the poplar panel ranges between 12.4 mm and 13.8 mm, its width between 53.3 cm and 53.4 cm, and its height h between 79.1 cm and 79.4 cm [START_REF] Dureisseix | Mechanical modeling of the activity of the flexible frame, chapter II.6[END_REF]. Therefore, neglecting also the curvature effects, the structure is usually modeled as a plate satisfying Kirchhoff-Love assumptions [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]. The 3D displacement field is therefore described with the displacement of the average surface displacement denoted with u for the displacement in the plane directions, and v for the transverse displacement, as

U = u + ve 3 -z grad v ( 65
)
e 3 is the normal vector along the thickness, and z is the associated coordinate. Here, the derivatives are performed with respect to the planar coordinates. The corresponding generalized strains are (e, χ), with

e = (Grad u) sym χ = -Grad(grad v) (66)
which allow to recover the plane part of the strain as

[ε] = e + zχ (67)
Moreover, the plane stress state is also assumed. Using Hellinger-Reissner variational formulation allows to define (i) the consistent generalized stresses (N , M ):

N = < [σ] > M = < z[σ] > (68)
where [•] and < • > denote the planar part and the average with respect to thickness, respectively; and (ii) the generalized constitutive relation, considering [START_REF] Kollman | Principles of wood sciences and technology[END_REF]:

N = < C cp ref > e+ < zC cp ref > χ + Ñ M = < zC cp ref > e+ < z 2 C cp ref > χ + M ( 69
)
where C cp ref is the elastic tensor that links the plane parts of stress and strain tensors with the plane stress assumption:

(C cp ref ) -1 = [C -1 ref ], and Ñ =< [ σ] >, M =< z[ σ] >.
To adapt the previous algorithm to this model, the most modular approach consists in using the discretized operators K, Bσ , Bε corresponding to the discretized plate model, here with DKT plate elements [START_REF] Batoz | An explicit formulation for an efficient triangular plate bending element[END_REF]. u now denotes the column vector of plate degrees of freedom.

The previous global solid box is replaced with the plate solver of algorithm 2, where the generalized forces for the plate problem are computed with [START_REF] Brémand | Measuring the relief of the panel support without contact, chapter II.3[END_REF], and by reciprocity, the plane strain ε is computed from ( 67), [START_REF] Allegretti | Barrier effect to water vapour of early European painting materials on wood panels[END_REF]. The corresponding structural operator is X S = Bε K -1 BT σ . Since the eventual outof-plane components of stresses and strains are not taken into account in the plate model resolution, X S is not invertible. Nevertheless the partition X D S = C -1 ref is still regular. Finally, for the termination criteria, C cp ref is used in (61).

Moisture evolution model

Concerning the moisture content evolution, on the present studied case, only one side of the panel is painted. This side will exhibit a much lower emissivity than the other side. Therefore, apart from some edge effects on the lateral sides, the main moisture gradient occurs along the thickness and neglecting the other directions lead to consider only 1D diffusion along the thickness, and the 1D mesh is a vertical segment. The scalar diffusivity for this problem is then D = e T 3 De 3 , e 3 being the normal to the plate. All the previous developments still hold for the fluid part in this case; the in-plane components of F are null. The previous global fluid box is therefore replaced with the 1D diffusion solver of algorithm 3.

The overall model coupling uses therefore the diagram of Table 3.

Algorithm 2 Global solid problem for the plate model

1:
Compute plate loading Ñ (k-1) , M (k-1) from 3D stress σ(k-1)

2:

Solve plate problem Ku

(k) = f di - BT σ [ Ñ (k-1) , M (k-1) ] to get u (k)
3:

Compute plate generalized strains [ẽ (k) , χ(k) ] ← Bε u (k)

4:

Compute 3D strain ε (k) from plate generalized strains ẽ(k) , χ(k)

5: Compute σ d ← (X D S ) -1 ε(k) + σ(k-1)
Algorithm 3 Global fluid problem for the 1D diffusion model

1:
Compute normal component of flux J(k-1)

n from 3D flux J(k-1) 2: Solve 1D diffusion problem M ∆w (k) = -g di + BT J J(k-1) n to get ∆w (k) 3: Compute W (k) ← B w ∆w (k) 4:
Compute normal component of moisture gradient

F (k) n ← BF ∆w (k) 5:
Compute 3D moisture gradient F (k) from its normal component

F (k) n 6: Compute J d ← (X D F ) -1 F (k) -J(k-1) Solid (plate model) Interface (3D) Fluid (1D) generalized strains (e, χ) ---------→ ε projection generalized stresses ( Ñ , M ) ← --------- σ restriction W ← --------→ W identity F ← ---------F n projection J ---------→ J n restriction
Table 3: Overall sketch of the model coupling strategy

Pictural layer

With aging, the pictural layer often became fragile, and a network of fine cracks (craquelures) usually developed and hopefully became stable during the past history of the painting [START_REF] Ravaud | The complex system of fine cracks, chapter II.2[END_REF]. This coating possesses its own permittivity to moisture. Here, we assume the pictural layer to be homogeneous at macro scale (homogenized from the craquelure network), and the thickness of the layer negligible when compared to the thickness of the panel, therefore, without adding additional mechanical stiffness.

As for the bulk diffusion, there are several modeling choices for the diffusion through the thin pictural layer. The normal massic flux j = J • n is linked to the jump in the state variable from one side (the wood panel side) to the other (the environmental air) [START_REF] Goossens | The measurement of the moisture transfer properties of paint films using the cup method[END_REF]. Denoting with • the jump of a quantity, i.e. its value at point M 2 minus its value at point M 1 , Figure 6, and since there is a water vapor on both sides, one can model the flux as j = -l ∆µ or j = -m H/H ref or, when using the sorption curve for the wood side, and a same fictitious wood on the other side, j = -ρ o k ∆w . The difference in all these models is the dependency of the material parameter l, m or k on the current state of the layer.

Note that within the framework of the proposed algorithm, each model can be implemented quite easily. For the problem we consider here, the state value at point M 2 is prescribed, usually in term of relative humidity: at z = e/2, H = H d (t). Therefore the flux can be split into a prescribed value and a value depending on the state at point

M 1 : j = -j d + mH 1 /H ref , where j d (t) = mH d (t)/H ref , or j = -j d + l∆µ 1 or j = -j d + ρ o k∆w 1 .
The last case, with a constant parameter k is the simplest one, since it directly depends on the primal unknown ∆w 1 = ∆w(z = e/2); therefore the flux

j d = -ρ o k∆w d , with ∆w d (t) = f (H d (t)/H ref )
, is added to the right hand side, and the diffusion matrix on the left hand side is summed up with ρ o k for the suited degree of freedom at z = e/2. The other cases can be implemented as well. The easiest way to implement this pictural layer effect is to add a special element to the diffusion mesh: an element with a null length (similar to a joint element for 2D or 3D problems), that will add a degree of freedom for moisture content (the one to be prescribed to ∆w d ), and for which the gradient F is the jump in moisture content [∆w d -∆w 1 ]. There is an additional unique integration point for this element, and provided that operators B J , B F and D take into account this additional element, the framework is unchanged.

Here, we consider the model using the moisture content jump, with a parameter k dependent on the surface strain s = Tr[ε(z = e

2 )] = Tr(e + e 2 χ). In this case, we proceed to the splitting (similar to ( 41), ( 42)):

j = -ρ o k( s )[∆w d (t) -∆w 1 ] = -ρ o k 0 [∆w d -∆w 1 ] +  (70)
where

 = -ρ o [k( s ) -k 0 ][∆w d -∆w 1 ] (71)
which is computed during the local coupled problem as (k) , and added to the right hand side during the global fluid problem, while the left hand side is added to ρ o k 0 . The evaluation of surface emissivity ρ o k is the most difficult part. Since each workpiece of cultural heritage is a unique object, and due to the difficulty to proceed to tests on these masterpieces, experiments can be conducted on mockups, designed as close as possible to the original studied object. In the scope of this article, we restrict ourselves to and estimation of its value considering the available data in the literature, and some crude micromechanical scenarios.

To get the form of the phenomenological model linking the emissivity k and the strain on the painted surface, we therefore consider a thin layer of either paint and varnish or only varnish (at the fine crack positions), with a given crack density. With a one-dimensional normal diffusion assumption (which is not the case due to the 3D diffusion at the layer thickness scale, but which is sufficient to derive the searched expression), a mixture law gives k = (1 -η s )k 1 + η s k 2 where k 1 (resp. k 2 ) is the emissivity of the painted surface, of area S 1 (resp. the cracked area S 2 ), and the surface ratio η s = S 2 /(S 1 + S 2 ). The surface area modification due to the strain is s . If there is a perfect mechanical link between the pictural layer and the panel, S i = S 0 i (1 + s ), where S 0 i is the area at reference configuration (i = 1, 2); therefore η s and the emissivity are constant: η s = η 0 s and k = k 0 . The other extreme alternative considers that the deformation is localized in the varnish only:

S 1 = S 0 1 while S 1 + S 2 = (S 0 1 + S 0 2 )(1 + s ).
In this case, one gets η s = η 0 s + s and k = k 0 (1 + κ s ), with κ = (k 2 -k 1 )/k 0 . Therefore, we propose the following dependency expression

k = k 0 (1 + ω s κ s ) (72)
with 0 < ω s < 1. In the following, we choose ω s = 0.5, κ = 0.2. With this dependency, a second kind of coupling between moisture and elasticity occurs: a structural coupling via boundary conditions, in addition to the previous state and evolution coupled constitutive equations.

In [START_REF] Engelund | The influence of moisture content on the water vapour resistance of surface coated spruce[END_REF], for spruce (Picea abies) and a water based white acrylic paint (60 g/cm 2 ), the coefficient m is identified: m = 2.84 10 -6 kg m -2 s -1 . Since it is obtained from a single state point (H 1 = 82 % and H 2 = 93 %), the corresponding value at that point of the other parameter is ρ o k = 8.8 10 -6 kg m -2 s -1 . Note that, in [START_REF] Olek | Inverse analysis of the transient bound water diffusion in wood[END_REF][START_REF] Olek | Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood[END_REF], on a straight wood surface without treatment, for Scots pine wood (Pinus sylvestris L), a surface emissivity is nevertheless measured on a whole range of relative humidity (from H 1 = 40 % to H 1 = 75 %, while H 2 = 75 %) and the corresponding parameter is ρ o k = 2. 10 -4 kg m -2 s -1 , i.e. an hydraulic resistance of 4 % of the previous value. Therefore, we neglect the hydraulic resistance of back panel side in all of the following. Finally, [START_REF] Gettens | Painting materials, a short encyclopedia[END_REF] reports for a cellulose acetate coating with a thickness of 0.5 mm, an equivalent to the parameter ρ o k with a value approx. 1.7 10 -7 kg m -2 s -1 . Several typical painted layer replica have been tested in [START_REF] Allegretti | External resistance to water vapour transfer of varnishes on wood[END_REF][START_REF] Allegretti | Barrier effect to water vapour of early European painting materials on wood panels[END_REF] for which the values of ρ o k 0 range between 1.6 10 -6 kg m -2 s -1 and 8. 10 -5 kg m -2 s -1 . In the following computations, we select the average value ρ o k 0 = 4. 10 -6 kg m -2 s -1 .

Problem settings

The case of Mona Lisa is used here to test the proposed solving strategy. The material coefficients are those given in the previous sections and, to avoid numerical ill-conditioning, a particular unit system is used. Concerning the fluid part, the goal is to have the three following values close to unity (these last two quantities lead to have the same magnitude for the diagonals of the mass matrix and of the diffusion matrix):

• the diffusion characteristic time τ = 1 D R e 2 π 2 ,
• the characteristic surface mass density m = 1 3 ρ o e, • the characteristic diffusivity d = ρo e D R . The length unity has been chosen to be L F = 1 meter, the mass unity, M F = 1 kilogram, and the duration unity T F = 12 hours. In such a way, τ ≈ 1.59

T F , m ≈ 1.95 M F L -2 F , d ≈ 0.37 M F L -2 F T -1
F . Unfortunately, this system is not suited to the solid part. Since solvers are decoupled, a different system could be used for each physics. To keep the same geometry and loading time evolution, we choose L S = L F and T S = T F . Moreover we expect the characteristic stiffness r = E L /h to be close to unity. This leads to choose the unity for pressure equals to 1 GPa, to get r ≈ 1.27 M S L -2 S T -2 S . Finally, some care must be taken when coupling the physics, with different unit systems. For the concerned case, this is not an issue, since only the moisture content (adimensional) is needed for the solid part coupling. For the fluid part coupling, σ/σ 0 (adimensional) is required for the body diffusion and ε (adimensional) for the painting diffusion.

Concerning the geometrical parameters, we choose here to use a constant thickness e = 13 mm. A stabilized crack in the upper left part of the panel is taken into account (but not the double-dovetail walnut braces performed during an early restoration act [START_REF] Mohen | Mona Lisa, inside the painting[END_REF]). Note that the plate model coupled with the 1D diffusion model are unable to take into account 3D effects that may appear at the crack tip. Since the stability of this crack is not the goal of the present analysis, this is not an issue for the following results. Figure 6 depicts the diffusion problem, with associated boundary conditions. Figure 5 illustrates the position the single wide poplar panel has been initially sawn in the trunk [START_REF] Ravaud | The Mona Lisa's wooden support, chapter II.1[END_REF], and the current shape at reference state (outside the frame) [START_REF] Brémand | Measuring the relief of the panel support without contact, chapter II.3[END_REF]. Figure 7 depicts the elastic problem, for which the loading is prescribed displacements to model the interaction of the painting panel with its housing frame (considered as rigid). The contact areas have been located with examining the real workpiece and are either due to contact with the frame itself, or with the backside crosspieces that squeeze (up to a certain point) the distorted panel. The values of the prescribed displacement have been identified from force measurements [START_REF] Uzielli | Physical and mechanical characterization of the support, chapter II.4[END_REF] and an elastic analysis.

The studied time interval is [0, T ] with T = 24 h. Two loading parameters are involved: the first one is related to the prescribed displacement at contacts with the housing frame and with the horizontal crosspieces; the second parameter is the time evolution of the relative external humidity H(t), see Figure 8. The first stage consists in a drying evolution, as happens when a climatic regulation is set up. This stage ends with a rapid variation in relative humidity, similar to a failure in the regulation system. When the system is reconnected, the drying begins again, and two short-time increase in relative humidity mimic two inspections of the panel, for which it is getting out of the regulation system. No surface tractions are considered, F d = 0, and the body forces are specific gravity.

The initial form of the panel is taken into account from measurements performed with shadow Moiré [START_REF] Brémand | Measuring the relief of the panel support without contact, chapter II.3[END_REF], Figure 9.

NUMERICAL RESULTS

All the simulations are performed with an explicit Euler method (θ = 1) as a time integration scheme, and the time discretization is the same for both the solid and the fluid physics.

The mesh for the solid problem uses DKT plate elements with Hammer integration points (for bending, and the associated 3-node triangle for in-plane behavior).

For the fluid problem, the 1D discretization in the thickness of the plate uses linear elements. There are as many 1D 'segments' as there are integration points for the solid problem (i.e. 3 times the number of plate elements) and the same transverse mesh is used to perform averages in the thickness (for instance, for the plate behavior deduced from the 3D Hooke tensor field). Though the diffusion problem is a parabolic system, the problem is solved with finite elements. In such a case, it requires an 'accuracy condition' [START_REF] Hibbitt | Abaqus/Standard -User's Manual[END_REF][START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF] ∆t (∆l

) 2 ≥ 1 6θD ( 73
)
where ∆l is a characteristic element size, and ∆t the time step. With the previous expression of τ , the number of time steps n τ per characteristic time τ is limited by the number of elements in the thickness n F as n τ ≤ 6θn 2 F /π 2 . Concerning the splitting in the proposed strategy, the simplest choice is selected:

X D S = C -1 ref and X D F = 1 ρo D -1
ref (and the corresponding k 0 for the pictural layer). The termination thresholds for the inner iterations are chosen as S = F = 10 -6 .

To assess the performances of the proposed strategy, the influence of the discretization parameters are tested, both on the solution and on the behavior of the algorithm.

Pre-processing is performed with the finite element code Cast3M [START_REF] Verpeaux | CASTEM 2000 : une approche moderne du calcul des structures[END_REF], processing uses two codes (elastic, and diffusion) implemented within matlab TM platform, post-processing is performed with gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF].

Influence of the fluid spatial discretization

Considering first the plate discretization with n S = 509 finite elements (Figure 10, on the left); 4 different discretizations for the moisture diffusion in the thickness are tested: n F = 4, 6, 8 or 10 elements of regular size, with a linear interpolation, for a fixed time discretization with n T = 120 time steps.

The numerical results are reported first as the evolution of the deflection at the central point on the panel, measured from the initial position of the reference state, Figure 11, for the different discretizations. The different stages in the humidity loading are clearly visible. The scheme is indeed convergent with respect to the fluid discretization. With few fluid elements, the moisture profile in the thickness tends to overestimate the water penetration, and the swelling/shrinkage is itself overestimated: This is why the central deflection is larger with a coarse fluid discretization. With these results, a fluid discretization with n F = 8 elements is selected as sufficiently accurate and will be used in the following.

Table 4 reports the average number of iterations per time step, and the cumulative number of iterations n it , to declare the convergence of the partitioning scheme at each time step. The number of iterations hardly depends on the fluid discretization. It is almost constant for all the tested cases.

Influence of the time discretization

With a fixed number of fluid elements (n F = 8) for moisture diffusion, and still n S = 509 elements for the plate discretization, 3 different time discretizations were tested: n T = 63, 120 and 240 time steps, with as regular time-step lengths as possible while consistent with the occurrences of the different loading stages.

The obtained central deflection is reported in Figure 12. Few dependence is observed; therefore the time step length is already sufficiently small for the precision criteria (n T = 120 will be used in the following).

Nevertheless, for the previously mentioned criteria for the time step length selection with a parabolic equation [START_REF] Amini | Multi-scale domain decomposition method for large scale structural analysis with a zooming technique: Application to plate assembly[END_REF], this time step length should not be too small with respect to the spatial discretization. With θ = 1 and the initial value D ≈ D R = 30. 10 -10 m 2 .s -1 , this criteria leads to ∆t/(∆l) 2 > 5.55 10 7 s.m -2 . which is the case here (∆t = 12 min for n T = 120 and ∆l = 1.625 mm for n F = 8, for which ∆t/(∆l) 2 ≈ 27.3 10 7 s.m -2 ). Note that due to the coupled model, if the diffusion coefficients decrease, the characteristic time increases; therefore, the accuracy condition may be violated. This is also an indication of the interest for using an adaptive time step (this feature is not implemented here).

The convergence in terms of iteration count is also reported in Table 4. Obviously, this iteration number per time step decreases when decreasing the time step length, but not in large proportions.

Influence of the solid spatial discretization

Finally, several discretizations for the plate problem are now considered, with n S = 509, 1697 and 6788 elements, using n F = 4 finite elements for the fluid, and n T = 120 time steps, Figure 10. Again, Figure 13 reports the central deflection, and convergence of the iterations are given in the same Table 4.

Only a weak influence is exhibited. Indeed, in the case n S = 6788, the minimum size of the elements is already smaller than the panel thickness. Therefore, as a plate model, the previous discretization is As a general remark, the average convergence rate is only weakly dependent on the element size of the solid discretization. The number of iterations required to reach the convergence with the criteria S = F = 10 -6 are not constant with respect to the time step, but their distribution is always similar to the case depicted in Figure 14. The main convergence difficulties (the largest number of iterations per time step) happen when the moisture content does vary rapidly, and for its highest amplitude, i.e. when the coupling is maximum.

The convergence of the error indicators η S and η F , for same case, are reported in Figure 15, for the time steps corresponding to the smallest and the largest number of iterations to converge. Asymptotic convergence rate is close to 0.25 (i.e. log(η (k) /η (k+1) ) ≈ 0.25) both for the fluid and the solid quantities; The previous definitions of the error nevertheless lead to a higher value for the fluid error than for the solid one. When the number of iterations at convergence is larger, the same convergence rate was obtained, but the initial levels of both errors were larger. Therefore, a possible improvement would consist in a better initialization of the iterations. Nevertheless, this is not under the scope of this article.

Interpretation of the numerical results

Considering the results obtained with the discretization parameters n S = 1697, n F = 8 and n T = 120, the distributions of moisture ∆w in the panel thickness are depicted on Figure 16, at different time steps. Each integration point of the plate model has its own evolution, and each plot reports the ranges between the minimum and maximum values. The jump in the values at the right end of the curves are due to the pictural layer with its own permeability, modeled as a zero thickness layer. Finally, the minimum principal surface strain on the painted side (eventually compressive, which is the most dangerous for the pictural layer) is depicted on Figure 17, at several time steps. These values are obtained for low moisture content which are therefore the most penalizing.

Other versions of the proposed strategy

Slightly different choices can be made for the design of the previous scheme. Apart from different splittings in [START_REF] Kollman | Principles of wood sciences and technology[END_REF], [START_REF] Olek | Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood[END_REF] and [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems. Part II: Energy transfer analysis and three-dimensional applications[END_REF] that can be chosen, the choice of D 0 could be made as a uniform value.

Additionally, a relaxation can be added to the previous algorithm 1: With an additional relaxation parameter γ in the same way the classical Jacobi algorithm can be relaxed in the linear case (if 0 < γ < 1, this is an under-relaxed scheme; if 1 < γ < 2, this is an over-relaxed scheme; γ = 1 recovers the previous algorithm). To do so, the lines 18 and 19 have to be modified as follows:

ε(k) ← γB ε u (k) + (1 -γ)ε (k-1)
(74)

σ d ← (X D S ) -1 ε(k) + γ σ(k-1) + (1 -γ)σ (k-2) (75)
and the lines 22, 23 and 24 as

W (k) ← γB w ∆w (k) + (1 -γ)W (k-1) (76) F (k) ← γB F ∆w (k) + (1 -γ) F (k-1)
(77)

J d ← (X D F ) -1 F (k) -γ J(k-1) -(1 -γ) J(k-2) (78)
A relaxation parameter slightly less than 1 is usually suited to improve the performance. The Table 5 reports the total number of iterations for the reference case when changing the relaxation parameter.

Best improvement has been found to be close to 30 % with γ = 0.95.

CONCLUSIONS AND PROSPECTS

The partitioning framework that is proposed in this article is well suited to multiphysics problems, for which it focuses on the interactions between the different physics. The modularity of such an approach allows to couple different models for the different physics, and is expected to cope with different discretizations as well [START_REF] Dureisseix | Information transfer between incompatible finite element meshes: Application to coupled thermo-viscoelasticity[END_REF]. Its convergence properties were little influenced by the discretization parameters, for the non linear coupled problem of hygroelasticity in wood structures under study. The extension of this approach to other coupled multiphysics problems is on development.

Apart from the multimodel management ability of the proposed scheme, exemplified on the coupling of a plate problem with a monodimensional diffusion problem, work is on progress on its usage to couple two different discretization schemes, for instance finite volumes and finite elements. In each case, an adaptive discretization scheme (eventually independently for each physics) would be an interesting feature.

Concerning virtual testing of wooden artworks of cultural heritage, the material behavior used here mainly couples elasticity and hygrometry with the simplest model. For such applications, a meaningful model would take into account moisture-dependent visco-elasticity, mechano-sorption, and eventually their dependency on temperature. Additionally, more details of interest can be taken into account. For instance the crack tip may require a finer 3D discretization, and the previous approach could be coupled with a structural zooming technique, such as in [START_REF] Amini | Multi-scale domain decomposition method for large scale structural analysis with a zooming technique: Application to plate assembly[END_REF]. Concerning the contact points with the housing frame, due to the movements of the panel, their location may evolve on the panel edges. These unilateral interactions are expected to be integrated as well in the framework proposed in this article, similarly as in the NSCD approach [START_REF] Moreau | Numerical aspects of sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF].

A Constitutive behavior

Other choices can be made to express the constitutive behavior; for instance, the material state can be described with other state variables, such as the stress σ and the chemical potential ∆µ. The same reversible part of the model can be obtained with a Fenchel transform of the previous potential Ψ(ε, ρ o ∆w) [START_REF] Baronas | Modelling of moisture movement in wood during outdoor storage[END_REF], or with a direct modeling such as:

Ψ (σ, ∆µ) = 1 2 σ : C -1 σ + B 1 (∆µ) : σ + Ψ 1 (∆µ) (79)
Once the state laws are derived, an identification can be also performed as: and neglecting the effect of the stress gradient, one gets:

D = L Ψ 1 + 1 ρ o Ψ 1 A : σ (83)
This expression can be used to select a form of coupling between the different diffusion tensors.

E

  L = 10.06 GPa E R = 1.19 GPa E T = 0.58 GPa G RT = 0.20 GPa G T L = 0.64 GPa G LR = 0.86 GPa ν RT = 0.70 ν LT = 0.47 ν LR = 0.36 w m = 7.0 % K = 0.76 C = 6.0 a L = 0.8 % a R = 9.5 % a T = 48 % D L = 30. 10 -10 m 2 .s -1 D R = 2.5 10 -10 m 2 .s -1 D T = 0.54 10 -10 m 2 .s -1 c L = -1.74 c R = -0.66 c T = -0.16

Figure 2 :

 2 Figure 2: Evolution of massic chemical potential vs. moisture content

Figure 3 :

 3 Figure 3: Evolution of diffusivity vs. moisture content

Algorithm 1 algorithm 1 :

 11 Main Initialization (first time step, i ← 1, initial conditions) 2:

  2 sin(2πt/T ) and H ref = 65 %. The studied time interval [0, T ] with T = 20 days is discretized in 10 time steps. No external forces are considered, and only 6 displacement dofs are clamped to avoid rigid body motions.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 5678910 Figure 5: Mona Lisa panel sawing within the trunk, local material basis (R,T,L) at a current point M , and initial deformed shape (amplified 5 times)

Figure 11 :

 11 Figure 11: Influence of the moisture diffusion discretization on the central panel deflection evolution (n S = 509, n T = 120)

Figure 12 :Figure 13 :

 1213 Figure 12: Influence of the time discretization on the central panel deflection evolution (n F = 8, n S = 509)

Figure 14 :Figure 15 :

 1415 Figure 14: Number of iterations to reach convergence ( S = F = 10 -6 ) along the time steps (n S = 1697, n F = 8, n T = 120)

Figure 16 :Figure 17 :

 1617 Figure 16: Distributions of moisture variations in the panel thickness at various time steps (n S = 1697, n F = 8, n T = 120)

Table 5 :

 5 Cumulative number of iterations n it as a function of the relaxation parameter γ (n S = 1697, n F = 8 and n T = 120)

Ψ 1 ( 1 +

 11 ∆µ) = ρ o f exp M ∆µ RT and B 1 (∆µ) = 1 ρ o AΨ 1 (∆µ)(80)Concerning the dissipative part of the model, the counterpart of (23) reads:J = -ρ o D grad ∆w = -ρ o D ∂∆w ∂∆µ B 1 : σ = Ψ 1 + 1 ρ o Ψ 1 A : σ (82)

Table 2 :

 2 Average number of inner iterations per time step, n N for the Newton solver, n P for the proposed partitioned solver, with respect to the coupling strength (and α = 20 • ).Hooke stiffness C = C ref (but still the same swelling / shrinkage A), such that the solid equilibrium reads Ku -N ∆w = f d

	38:	end if
	39:	until (η	(k) S ≥ S or η	(k) F ≥ F )
	40: end for

Table 4 :

 4 Convergence indicators (average number of iterations per time step and cumulative number of iterations n it ) as functions of the parameter discretizations (number of elements for the diffusion through the thickness n F , number of time steps n T , number of plate elements n S )

	n S	n F n T n it /n T	n it
	509	4	120	16.7	2009
	509	6	120	16.9	2030
	509	8	120	17	2040
	509	10 120	17	2037
	509	8	63	18.2	1166
	509	8	120	17	2040
	509	8	240	15.2	3646
	509	4	120	16.7	2009
	1697	4	120	17	2034
	6788	4	120	16.7	2006
	1697	8	120	17	2039

http://www.woodculther.org

Acknowledgements:

The second author acknowledges fundings from CNRS / Région Languedoc-Roussillon, and from Program Vinci 2007 (Université Franco-Italienne). Both authors wish to thank their collegues L. Uzielli, J. Gril, F. Cherblanc, J.-C. Benet and A. Chrysochoos for fruitful discussions. This work has been supported by European COST Action IE0601 'Wood Science for Conservation of Cultural Heritage' 1 .