Numerical modelling of 3D dowelled timber joints using advanced fully coupled hydro-mechanical constitutive equations

– In this work, a constitutive law describing the hydro-mechanical behaviour of wood is presented. The main novelty of developed 3D F.E. model, among literature hardening models, is the coupling between plasticity and moisture eﬀect. The hydro-mechanical wood behaviour is analyzed within the framework of non-linear mechanical theory of Mechanical Continuous Media (MMC). The local integration of constitutive laws is carried out with an iterative implicit outline using plastic corrections and the calculus procedure is implemented on Abaqus software and validated in the context of dowelled timber joints. The results of simulation were compared with published experimental results, showing good agreements and demonstrating clearly the eﬀect of the moisture on the plastic behaviour.


Introduction
Wood is composed of polymers which are sensitive to moisture. Therefore to study its hydro-mechanical behaviour is indispensable to develop its usage on construction field. The mechanical behaviour of wood is largely governed by its physical state. Originally, wood contains a certain amount of water that circulates within, which is necessary for its physiological growth. Hence, the moisture content is a fundamental characteristic that affects its mechanical properties. When moisture level lies between 30 to 35%, there is no longer water in liquid phase in the cell cavities. This zone corresponds to the fibre saturation area (FSA).
For several years the hydro-mechanical behaviour of wood has been analysed by many researchers [1][2][3][4][5][6][7][8]. More recently, Saftand et al. [2] proposed a hygro-mechanical coupled macroscopic material model to use within the framework of the FEM. Wood plastic behaviour is described by a linear strain hardening function. In the present work, this behaviour is analysed within the framework of nonlinear mechanics of continuous media. A theoretical model based on the thermodynamics of irreversible processes is used to describe the hydro-mechanical coupling between mechanical properties and moisture effect. The effect of moisture is incorporated to an elastoplastic model through the improved equations proposed by Guitard [3], which reflect the effect of moisture on the orthotropic elastic properties of timber. The effect of moisture on plastic flow is included through the expression of the yield function in terms of moisture content (Eq. (7)).
Timber is made up of polymers susceptible to humidity. Many parameters, such as moisture, influence on its mechanical behaviour. The mechanical behaviour of wood under the mixed action of water and mechanical effects is called hydro-mechanical. Characterization and study of HM behaviour of wood are needed to predict timber behaviour during his life time. To date, no commercial calculus code which is able to predict the behaviour of wood connections taking into account the effect of water content is known. Our work is dedicated to this aspect: the hydro-mechanical behaviour of dowelled timber joints.
Article published by EDP Sciences The new contributions in this paper in relation to our previously published works in references [9][10][11] can be summarized as: -Guitard equations described in [3] are included in the proposed model to take into account the rate of moisture on the orthotropic elastic properties of wood. We generalized these equations in plastic field correcting the form of the yield in our model developed in references [9][10][11]. -The analytical formulation of the proposed model based on an incremental approach is different from those used in references [1][2][3][4][5][6][7][8], it is considered as an approach "advanced" compared to those available in the literature. In addition, the model proposed in this study is rich in state variables; it contains 28 parameters including 9 for elasticity, 9 for plasticity and 2 for moisture.
The proposed finite element model (FEM) was implemented on Abaqus finite element code [12] to explore the hydro-mechanical behaviour of wood in dowelled timber joints. Its development has a particular interest for timber engineering industry, especially with the increase in demand of wood as a construction material. The main characteristics and performances of the model are the purpose of the discussions below.

Analytical modelling of hydro-mechanical behaviour
The plasticity model presented is based on the thermodynamic approach of irreversible processes with internal variables. The state relations are given as [9,10,[13][14][15][16][17][18][19]: σ is the Cauchy stress tensor, ε e the elastic strain tensor, Λ the fourth-order operator of elastic properties, R the scalar variable representing the isotropic hardening stresses, r the parameter of the isotropic hardening and and Q the isotropic hardening modulus. Λ is a function of , and Poisson's ratios ν 12 , ν 13 and ν 23 in the directions (1), (2) and (3) as shown in Figure 1.
The operator Λ is the inverse of the compliance operator C −1 : Guitard [3] adopted the following corrections for the fluctuations in moisture occurring around 12%: The plastic flow is described by a plastic potential, which is written as: with: where σ represents the equivalent stress, σ yield the yield stress, σ h yield represents the yield strength of a humid material, α h the parameter of relative humidity, h the moisture content and H is Hill's operator, which is a function of six parameters F , G, H, L, M and N [9][10][11]20]. It is given as: Using the normality rule, the complementary relations are derived as follows [9,10,[13][14][15][16][17][18][19]: where ε p is the plastic strain tensor; n the normal to the loading surface; b the coefficient of isotropic hardening andλ is the plastic multiplier.

Numerical implementation
Using a dynamic explicit scheme (DE), the constitutive equations are integrated over a time increment to 501-page 2 obtain the quantities σ n+1 and R n+1 . The plastic strains tensor εp n+1 is obtained using the θ-method while the isotropic hardening variable is r n+1 is obtained with the asymptotic method; that is: The plastic criterion f n+1 which should satisfy the consistency condition at time t n+1 is given as: Given a total strain increment Δε at time t n , assuming the response to be initially elastic, the trial stress at time t n+1 is obtained as: which is then substituted in the expression of the failure criteria, Equation (13), to yield: with If f * n+1 < 0, the obtained response is elastic, and the forces variables at time t n+1 are obtained as: Otherwise, if f * n+1 0, it becomes necessary to apply a correction to the trial stresses σ * n+1 to obtain σ n+1 and R n+1 such that: The final system of equations to solve then becomes: The model is implemented in user subroutine Vumat in Abaqus/Explicit [12]. In an explicit dynamic analysis, the global equilibrium equation becomes: where [M n ] is the mass matrix, and {R n } the residual vector.
The explicit algorithm consists in obtaining solutions at the time t n based on solutions obtained at the previous times: where {Ü n } is the acceleration vector, {U n+1/2 } is the velocity vector, and {U n+1 } is the displacement vector. The stability and precision of an explicit scheme depends very much on the chosen time increment.

Results and discussion
In this section, Vumat subroutine is applied, implemented on code Abaqus/Explicit, to simulate steel-totimber joints behaviour. The main objective is to show the abilities of this coupled "hydro-mechanical" approach to predict the strength of timber joints. The validation is done by the simulation of a single and a multiple-fastener timber joints tests.

Embedding tests
In this section the model is used to simulate the embedding tests of spruce timber as it is shown in Figure 2. Embedding tests parallel to the grain orientation were performed to determine the spruce timber strength with a dowel presence [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. The tests were made according to ASTM D5767. The loading is applied as it is shown in Figure 3, and only vertical displacement is allowed. Firstly, the specimen is loaded up to 40% of its strength, then it is unloaded to 10% with a rate of 0.35 mm.min, and finally loaded up to failure. The details of the experiments can be found in reference [7].
In the finite element model, steel rod is modelled with 550 shell elements type S4R available in Abaqus/Explicit. On the other hand, timber specimen is modelled with 10 900 brick elements type C3D8R with reduced integration. The mesh is shown in Figure 2b and material properties and all model parameters are given in Table 1. Contact between rod and sample is simulated using classical isotropic Coulomb friction model with a friction coefficient of μ = 0.23.
The model parameters leading to simulations results are summarized in Table 1.   Figure 3 shows the global responses on terms of load versus displacements plotted for h = 13%. Predicted results match closely to the experimental ones. For h = 13%, the overall response is linear up to a displacement of 0.07 mm corresponding to a force of 4.7 kN. Beyond this point, the response becomes nonlinear until the force reaches 12.2 kN at a displacement of 1.1 mm. Figure 3 shows the evolution of the equivalent stress at different displacements. At a displacement of 0.05 mm, the formation of stress concentration zones with circle shapes is observed on the contact region between rod and specimen. These zones expand as the rod penetrates into the sample. Other zones of stress concentration appear on sample sides at a displacement of 0.35 mm (F = 9.8 kN). These two zones progress and join each other at a displacement of 1.1 mm. Figure 4 shows the global responses in terms of load versus displacements plots for h = 13% and h = 26%. The predicted results match closely the experimental ones. For a moisture content of 26%, the behaviour is linear up to a displacement of 0.07 mm. Beyond this point, the behaviour becomes non-linear. The corresponding force is equal to 2.1 kN. For a displacement higher than 0.7 mm, the force is constant and equal to 5.1 kN, which is well predicted by the model. For a moisture content of 13%, the load-carrying ability response doubles to reach 10.8 kN at a displacement of 1 mm. This example illustrates very clearly the influence of moisture on the overall behaviour of timber.

Multiple-fastener connection
For this section, test results are taken from [36]. Specimen geometry and boundary conditions are given in For the simulation, a fourth of connection is discretized. The finite element model used is composed by a mesh of 54 800 elements of type C3D8R, it is shown in Figure 5. Friction coefficient between all components (wood specimen, plate and fasteners) is μ = 0.21.
The material parameters of the wood species are given in Table 2.
To model steel, an isotropic elastic behaviour coupled with perfect plasticity is considered. Steel material characteristics are summarized in Table 3.
Simulation results are compared with those of the results obtained by Laplanche in [36] (Fig. 6) Fig. 7a), which is consistent with experimental observations. Figure 7c shows the failure modes of fasteners by plastic coating. Fasteners plasticize with "V " shape, which corresponds to the experimental observations performed on connections of two shear planes.

Conclusion and perspectives
In this work, a constitutive 3D F.E. model for hydro-mechanical behaviour of wood is presented and successfully implemented in the commercial Abaqus software. A theoretical model based on the thermodynamics of irreversible processes is presented to describe the hydro-mechanical coupling between the elastic-plastic anisotropic properties and moisture for timber. Theoretical and numerical modelling were discussed in detail. The numerical integration of the 3D constitutive equations of the model was performed by implicit methods.
To validate the model, two different steel-to-timber joints tests were used: an embedding test and a multiplefastener connection. First simulation studies the embedding of a metal rod into a timber specimen, which is widely used in timber structures to simulate the behaviour of connections. Numerical results show a good agreement with the published experimental ones. Simulation results of multiple-fasteners joints approach successfully the connection embedding and the fasteners failure modes. The predicted results confirm that the model is capable of describing the hydro-mechanical behaviour of timber.
Research is now in progress and further improvements and developments, such as failure modes in tension, material temperature dependence, will be addressed in a future research work.