Normal contact performance of mortise and tenon joint: theoretical analysis and numerical simulation

This article aims to investigate the contact characteristics of mortise and tenon (M&T) joints in the traditional timber structures. In particular, the normal embedded compressive contact between contact surfaces of M&T joint was investigated. Based on basic contact theory and contact characteristics between mortise and tenon, a normal elasto-plastic contact model, which can reflect the real normal contact behavior of M&T joints in traditional wooden structures, was proposed. Coulomb friction contact was utilized to describe the tangential slipping characteristics of the contact surfaces. Micro-morphology scanning tests of wood samples with different roughness were carried out to determine the parameters involved in the normal contact model. The normal contact model subroutine of M&T joint was compiled by FORTRAN language, implemented into ABAQUS through user-defined interface (UINTER). Then the proposed model was verified by shear tests of wood contact surfaces considering different normal pressures. Finally, a finite element model (FEM) of straight tenon joint subjected to cyclic reversed loading, based on the proposed normal elasto-plastic contact model, was developed, and a FEM considering normal “hard contact” between the contact surfaces, was also performed. The simulation results were validated by the experimental results. Results showed that the user-defined normal elasto-plastic contact FEM was more in line with the actual force state and mechanical behavior of M&T joints, which can more accurately predict the failure modes and simulate the hysteretic behavior of M&T joints, compared to the FEM considering normal “hard contact” of the contact surfaces.


Introduction
Wooden structures have been widely used all over the world, notably in China, Korea, Japan, Canadian and USA, for the wide availability, low cost and ease of erection of wood [1][2][3]. In China, mortise and tenon (M&T) joint is the typical connection type of the traditional timber structures, as shown in Fig. 1a, which has a significant effect on performance of structure [4][5][6][7][8][9]. The types of M&T joint mainly include straight tenon joint and dovetail tenon joint (see Fig. 1b).
For M&T joints of the ancient timber structures, analyzing the force mechanism of M&T joints could be very difficult. The mechanical model of M&T joint subjected to moment and shear force is shown in Fig. 2 [10][11][12][13][14][15]. When the joint produces counterclockwise rotation, the contact surface I (CSI) of upper edge on the tenon neck and the contact surface II (CSII) of lower edge on the tenon head will be compressed by the edge of the mortise, which will generate normal contact force F 1 and F 2 , respectively. When the rotation angle further increases, the mortise and the tenon will produce the relative sliding because there is a certain gap between M&T joint, and the mentioned surfaces will transfer tangential frictional force f 3 and f 4 . The mechanical process of M&T joint mainly includes normal embedded compressive contact and tangential frictional contact.
The mechanism of normal compressive contact on rough surfaces of M&T joints is extremely complex. Figure 3 shows normal contact pairs of M&T joint. It can be seen that there are two kinds of normal contact between the mortise and the tenon: one is the contact between the upper and lower surfaces of mortise and tenon (m1-s1, m2-s2, as shown in Fig. 3); the other is the contact between the two sides of tenon and mortise (m3-s3, m4-s4, see Fig. 3). Generally speaking, the elastoplastic contact behavior lies more in the ductile nature and strain stiffening of wood in the perpendicular to grain direction under macroscopic condition. However, the rough surfaces of the mortise and the tenon on the microscopic scale are actually composed of many asperities [16,17], as illustrated in Fig. 4. In the real contact process of rough surfaces, the contact state of asperity will change with the increase of external load [18], which not only experiences elastic state, but also undergoes plastic state. To deeply understand the complex normal compressive contact mechanism of M&T joints, it is of vital importance to investigate the normal contact properties of the rough contact surfaces in the micro-scale.
In recent decades, to investigate the interface contact behavior of M&T joints in traditional timber structures, many researchers have studied the mechanical performance of M&T joints. Chang et al. [19] and Ogawa et al. [20] researched the effects of initial gap between members on the mechanical behavior of M&T joints, respectively. Feio et al. [21] and Ma et al. [22] revealed the compressive strength perpendicular to grain of wood,    The normal contact surfaces of M&T joint at the microscopic scale [16] between normal contact surfaces, had a remarkable influence on the ultimate load of M&T joints. Ma et al. [22] investigated the normal embedded compression mechanism of straight tenon joint under the cyclic loading, and derived the moment and rotation angle relationship of these joints.
Also, theoretical models of normal contact surfaces have been proposed for the investigation of the normal contact characteristics between rough surfaces. Chang et al. [23,24] pointed that the deformation of asperity will transfer from the elasticity to the plasticity with the increase of load, in which classical Hertz elastic contact theory [25], GW model [26,27] and AF model [28] were no longer applicable, in view of this, an elasto-plastic CEB contact model for rough surfaces was proposed. Gao et al. [29] presented a normal damping model of joint contact surfaces, based on the two rough surfaces model [27], taking into account the pre-slip state and asperities in lateral contact. Jamshidi et al. [17] and Ahmadian et al. [30], respectively, put forward a surface contact model considering the couplings between normal and friction forces of contact surface. Most previous studies of contact models of rough surfaces have only considered the elastic contact of two contact surfaces and not taken into account the coupling of three-state changes of the asperity on rough surfaces, i.e., the pure elastic, elasto-plastic, and full plastic states.
Furthermore, many researchers conducted numerical simulation on the contact behavior of rough contact surfaces on M&T and carpentry joints. Villar et al. [31] carried out the finite element (FE) simulation on cogged joints, and found that the friction coefficient, between the contact surfaces of components, had an obvious influence on the seismic performance of the structure. Li et al. [9] carried out numerical simulation analysis on the double span timber frame with the M&T joint by setting the nonlinear semi-rigid (NSR) link element in commercial FE analysis software SAP2000, investigating the hysteretic behavior of timber frame, however, the normal embedded compressive contact characteristics between the rough contact surfaces of M&T joints was ignored. Xie et al. [16] proposed a frictional constitutive model of contact surfaces, which can be utilized to describe the tangential contact performance between the mortise and tenon, the hysteretic performance of Chinese traditional timber frame with M&T joints was simulated by ABAQUS more exactly.
In conclusion, although many scholars have conducted a large number of theoretical research and numerical simulation analysis on the contact characteristics between the rough contact surfaces of M&T joint, the normal contact of rough surfaces was regarded as "hard contact" and the elasto-plastic deformation of asperity on the rough surfaces was neglected in the FE simulation research of M&T joint, which fails to exactly simulate the normal contact characteristics of M&T joint, nor reflect the corresponding normal contact mechanism from microscopic view.
In this article, on the basis of the classical contact theory and the normal contact characteristics of rough contact surfaces, a normal elasto-plastic contact constitutive model was proposed. The proposed model was able to characterize the embedded compressive contact characteristics of M&T joints in ancient wood structures. Constant friction coefficient was used to describe the tangential behavior of contact surfaces. The normal contact constitutive model was implemented into the commercial software ABAQUS as subroutine by using FORTRAN. The proposed normal contact model was firstly verified by rough contact surface shear test of wood, then a three-dimensional FE model of M&T joints was developed based on the proposed model, the modeling results were compared with that of the conventional modeling method (in which the normal "hard contact" characteristics of contact surfaces was considered) and the experiments.

Theoretical analysis on normal contact model of M&T joint
In this section, a normal elasto-plastic contact model of M&T joint was presented based on Hertz elastic contact theory [25] and GW model [26]. The three deformation stages of the asperity were considered in normal contact model, which are as following: (1) the Hertz elastic contact theory was employed in the complete elastic stage of the asperity; (2) the Abbott plastic contact theory was adopted in the full plastic stage of the asperity; (3) the power-exponent function was used to model and analyze in the elasto-plastic stage of the asperity. Finally, the method of mathematical probability and statistics was utilized to make the contact research between a single asperity extend to that of two rough surfaces [32,33], and the normal elasto-plastic contact model of rough surfaces in M&T joint was further established.

Normal contact model of single asperity
To simplify the deduction process of normal contact model of single asperity, the following assumptions are made: 1. The top of the asperity is spherical; 2. The curvature radius of all asperities on the rough surface is equal; 3. The interaction between the adjacent asperities on the same surface was not considered; 4. The height of asperity obeys Gaussian distribution [26,27].
The contact state between a single asperity and a rigid smooth plane in two-dimensional coordinates is illustrated in Fig. 5, where z is the asperity height, w stands for the normal deformation of the asperity, d is the distance between the average height line of the surface peak and the rigid smooth plane, and R stands for the asperity radius.
The deformation of asperity increases as the increase of normal contact load N, and the deformation process experiences three contact states: elastic, elasto-plastic and plastic. The simulation results investigated by Kogut et al. [34] showed that the full plastic contact will begin when the average contact pressure increased to a constant value equivalent to the hardness of the contact material, meanwhile, the normal deformation w will increase to the critical value of full plastic deformation accordingly, namely w = w 2 = 110w 1 , where w 1 is limit of elastic deformation. Finally, the conclusions were as follows: when w < w 1 , the asperity was in complete elastic deformation; when w ≥ w 2 , it went into full plastic deformation; when w 1 ≤ w < w 2 , it is in the elasto-plastic deformation state. The normal contact models of asperity under different contact states are analyzed in the following sections, respectively.

Pure elastic deformation stage
When the normal deformation w < w 1 , the asperity is in the stage of pure elastic deformation. On the basis of Hertz theory [25], the actual contact area A e is equal to the truncation of the undeformed surface profile at the intersection of a single asperity and a rigid smooth plane, which is expressed as Normal contact load N e is calculated by the following equation The average contact pressure P e is defined as where where E* represents the equivalent elastic modulus; E 1 and E 2 stand for the elastic modulus of the two contact objects, respectively; ν 1 and ν 2 represent the Poisson's ratios of two contact objects, respectively.
The contact pressure at the beginning of plastic deformation P p, initial is associated with hardness of material [23,24], among, the maximum elastic limit contact pressure P e, max is equivalent to the contact pressure at the onset of plastic deformation P p, initial , which can be expressed as With Eqs. (3)-(5), the critical normal deformation w 1 corresponding to the initial yield point of plasticity can be derived as in which H stands for the hardness of softer material; k represents the mean contact pressure coefficient. Tabor [35] introduced an empirical relationship between P p, initial and the hardness H for various materials. After calculation, the asperity began to yield when the mean contact pressure factor k is equal to 0.4.
With Eqs. (1), (2), (3) and (6), the normal contact load N 1 , actual contact area A 1 and average contact pressure P 1 corresponding to the critical normal deformation w 1 can be, respectively, derived as

Full plastic deformation stage
When the normal deformation w ≥ w 2 , the asperity is in the state of full plastic deformation, and the mean P e, max = P p, initial = kH .

Fig. 5
Normal contact model between a single asperity and a rigid smooth plane contact pressure value P p is equal to the ratio of contact load N p to contact area A p , which is also equivalent to the hardness value H of the contact material [28,35], namely The critical deformation of complete plasticity w 2 is defined as [34] According to theory deduced by [28], the actual contact area A P is expressed as The normal contact load N P is calculated as

Elasto-plastic deformation stage
When the normal deformation w 1 ≤ w < w 2 , the asperity is in the stage of elasto-plastic deformation. The whole deformation process of asperity on rough surfaces presents continuous and smooth characteristics. Therefore, it can be assumed that the actual contact area and contact pressure, at the two normal critical deformations w 1 and w 2 , are continuously and differentiable. In this section, the form of power-exponential function was adopted to carry out the analysis of elasto-plastic deformation [36], which can be expressed as where P * ep , A * ep and N * ep , respectively, represent dimensionally normalized average contact pressure, contact area and the contact load in the state of elasto-plastic deformation; w * is the ratio of normal deformation w to the critical normal deformation w 1 corresponding to the initial yield point.
According to considering the continuity condition at the critical deformation of initial yield w 1 and the critical deformation of full plastic w 2 , the following equations can be obtained: Combining with the equations from Eqs. (14) and (18), unknown parameters can be obtained Substituting Eqs. (7)-(9) and (19) into Eqs. (14)(15)(16), the expressions of the contact load N ep , the contact area A ep and the average contact pressure P ep in the stage of elasto-plastic contact deformation can be expressed as Eqs. (20)-(22), respectively:

Normal elasto-plastic contact model between two rough surfaces
For the contact between two rough surfaces, which can be equivalent to the contact between equal rough surface and rigid smooth surface [27], as shown in Fig. 6, where d represents the distance between the rigid smooth surface and the equivalent rough surface. The contact condition between rough surface and rigid plane is z > d. In addition, the rough surface is composed of uneven asperities randomly, and the contact between the two rough surfaces can be regarded as the contact between the asperities. According to assumption(4), the peak height of the asperities

Fig. 6
Normal contact between two rough surfaces [27] obeys the Gaussian distribution (standard normal distribution) (see Fig. 6), and the expression is as follows where σ is the standard deviation of height distribution function.
On the basis of the contact theory of a single asperity and the method of probability and statistics, the total contact area and the total contact load between rough surfaces can be obtained.

Total contact area of rough surfaces
Assuming that there are N 0 asperities between the two rough surfaces, the total number of asperities n in actual contact can be described as Eq. (24), based on Eq. (23): where η represents the number of asperities per unit area; A n stands for the nominal contact area.
Because there are three contact states between asperities on rough surfaces, the total contact area of rough surfaces is equivalent to the sum of the actual contact area of all asperities under the three contact states. With Eqs. (8), (12), (21) and (24), the total contact area of rough surfaces A can be calculated as where A e , A ep and A p are the contact area of asperities in elastic, elasto-plastic and full plastic states, respectively.

Total normal contact load on rough surfaces
Similarly, the total normal contact load of rough surfaces is equivalent to the sum of the actual contact load of all asperities in the three contact states. With Eqs. (7), (13), (20) and (24), the total contact load N of rough surface can be obtained, which is described as Eq. (26): Furthermore, the Coulomb friction model is adopted to describe the tangential slipping contact characteristics of contact surfaces, namely, the friction is proportional to the normal load between the contact surfaces, the expression of the tangential friction f t can be described as Eq. (27): where μ represents the friction coefficient; N is the normal load between the contact area.

Determination of model parameter Material test of wood
Combined with the normal contact model of rough surfaces in M&T joint proposed in "Normal elasto-plastic contact model between two rough surfaces" section, the wood property parameters to be measured are as follows: the elastic modulus parallel to grain E L , the elastic modulus perpendicular to grain along radial direction E R , Poisson's ratios of transverse and tangential section ν LT and ν RT , hardness of transverse and tangential section H RT and H LT .
The relevant wood property tests had been presented previously [37,38] and are not be repeated for brevity. Larix gmelinii was used to make the small clear specimens, the material test results of which are shown in Table 1.

Table 1 Material properties of Larix gmelinii
Six samples were utilized for each test. E L , E R , E T are the elastic modulus in longitudinal, radial, and tangential directions, respectively; G LR , G LT , G RT are the shear modulus under different directions; ν is the Poisson's ratio under different directions; M represents the moisture content; ρ is the density of wood; f t,0 and f c,0 are the tension and compression strength parallel to grain, respectively; f c,90 stands for the compression strength perpendicular to grain.

Surface morphology parameters of wood
On the basis of the proposed normal contact model of rough surfaces existing in M&T joint, as can be seen the primary parameters used to describe the rough surface morphology include the height distribution function of asperity z(x, y), the peak radius of asperity R, and the per area density of asperity η.
For the sake of gaining the above parameters, the surface roughness parameters to be solved are as follows: arithmetic average height of asperity S a , root average square height S q , peak vertex density of asperity S pd , peak arithmetical mean curvature of asperity S pc , kurtosis S ku . The specific meaning of the above parameters can be found in Chinese National Standard GB/T 33,523.2-2017 (2017) [39] and Franco's work [40]. Chen et al. [41] and Majumdar et al. [42] pointed that the asperity height of rough surface conformed to the standard normal distribution (see Eq. (23)) when kurtosis S ku was equal to 3, and the standard deviation σ was equal to the root mean square height S q . Then the computational expressions of the surface morphology parameters are shown as follows:

Surface micro-morphology test of wood
Northeast Larix gmelinii of China was used to make small specimens; the moisture content and density of wood were 13.3% and 0.65 g/cm 3 , respectively. According to GB/ T 14,495-2009 (2009) [43], 30 small cube samples with a 50 mm × 50 mm section size and with 50 mm in height, were fabricated to determine the surface morphology parameters. Each sample was divided into three sections: transverse section, radial section and tangential section, as illustrated in Fig. 7. Then different surface roughness can be obtained by sanding with sandpaper. The grinding methods were divided into no-grinding, 120 mesh sandpaper grinding and 400 mesh sandpaper grinding.
The scanning system of Leica three-dimensional (3D) confocal microscope (DCM3D), produced by Leica microsystems Co., Ltd., was utilized to observe the surface topography of the specimens, as shown in Fig. 8, which consisted of confocal microscope, data acquisition equipment and data analysis system, and whose  The scanning system of Leica 3D confocal microscope advantage is that it can obtain clear and fine surface micro-morphology without contacting the sample surface. The sampling parameters of the instrument were set as follows: the horizontal scanning range was 50 mm × 50 mm, the vertical measuring range was 15 mm, the sampling interval was 100 nm, and the scanning speed was 0.5 mm/s. For brevity, the scanning results of some samples were used for analysis, the 3D profile of three different sections under different processing methods are shown in Figs. 9,    10,11. It can be seen that the color of the three sections gradually transformed from warm color to cold color with the increase of processing accuracy, and the color was gradually consistent, indicating that the height of profile decreased and the surface became more and more smooth. It is also clear that the color distribution of cross section was deeper and presented uneven, under the same polishing condition, compared to that of tangential and radial section, indicating that the cross section was rougher than the tangential section and radial section. The micro-morphology parameters of wood surfaces are illustrated in Table 2, in which the values of the root average square height S q , the peak vertex density of asperity S pd , the peak arithmetical mean curvature of asperity S pc , the kurtosis S ku are automatically calculated by 3D scanning data analysis equipment (see Fig. 8), respectively, and the values of the peak radius of asperity R, the per area density of asperity η and the standard deviation σ are obtained based on the calculation expressions (Eqs. (28)- (30)). Furthermore, it is worth noting that the value of S ku is close to 3, which proves that it is reasonable to assume that the peak height of asperity obeys Gaussian distribution in the process of theoretical derivation in "Normal contact model of single asperity" section.

Subroutine development
Based on the normal elasto-plastic contact model of rough contact surfaces formulated in "Normal elastoplastic contact model between two rough surfaces" section, the user-defined interface (UINTER) subroutine was compiled by FORTRAN language according to the corresponding rules [44], as shown in Fig. 12, which can be implemented in ABAQUS, where the stress and displacement contact characteristics and the coupling relationships, between the normal and tangential, can be defined.

Validation of model Shear tests for rough contact surfaces of wood
A verification of the wood normal contact model subprogram was carried out by means of comparison to shear tests results of wood contact surfaces. Under the action of the normal stress, direct shear tests were carried out on two pairs of wood normal contact interfaces existing in the M&T joint, namely, transverse-tangential section and tangential-tangential section. The experimental program included 12 small cylindrical specimens having dimensions of 61.8 mm diameter and 10 mm height, as shown in Fig. 13, in which there were three groups of transverse-tangential contact surfaces and three sets of tangential-tangential contact surfaces.
To form normal stress of 0.05 MPa, 0.1 MPa and 0.2 MPa, respectively, three different weights were applied to each group. The direct shear apparatus controlled by strain was adopted as the loading device of the test, as shown in Fig. 14, which mainly consisted of upper and lower shear boxes, lever pressure equipment, rotatable handwheel as well as force and displacement gauge. After installing the shear box, the vertical pressure was applied by adding weights at the lever, and the horizontal displacement was applied to the lower shear box by turning the handwheel. The loading rate was controlled at 1.2 mm/min. The comparisons of the experimental results and the numerical results will be analyzed in the following section.

Finite element implementation
A 3D finite element analysis (FEA) for transverse-tangential contact surfaces of wood were carried out with geometry and load according to the direct shear tests described in "Shear tests for rough contact surfaces of wood" section. To represent the timber species, a local coordinate system was used to define the orthotropic properties of wood. C3D8R-type element was utilized, and the mesh size was set as 2.5 mm. The normal  13 Physical diagram of contact surface of samples: a transverse-tangential contact surface; b tangential-tangential contact surface pressure of 0.05 MPa was applied to the top surface of the specimen, the shearing process between contact surfaces was controlled by displacement; the FEA model of transverse-tangential contact surface of wood is illustrated in Fig. 15. The elasticity constants were used to characterize the orthotropic anisotropy of wood under the elastic state, the plastic compressive stress and strain relationship of timber parallel to and perpendicular to grain were used to simulate the plastic behavior of transverse section and tangential section, respectively, as shown in Fig. 16. All the material properties adopted in the model are summarized in Table 1. Further details can be discovered in the research carried out by Zheng [37] and Xie et al. [38].
The interaction between two normal contact surfaces in the model was also defined and modeled by the userdefined normal contact constitutive subroutine. Based on the surface micro-morphology parameters and wood material parameters measured by the experiments, the parameters of normal contact model subroutine selected are shown in Table 3, it is noteworthy that the friction coefficient between normal contact surfaces was still set as 0.5 to describe the tangential friction slipping characteristics.   Figure 17 draws the normal contact stress and analysis step time relation curve of wood contact surfaces obtained from numerical simulation. It can be found that the normal contact stress and analysis step time relationship curve of the contact surfaces showed strong nonlinear characteristics, which can be roughly divided into three states, that is elastic phase (OA phase), elasto-plastic phase (AB phase) and plastic phase (BC phase). This indicates that the FE simulation of wood rough surfaces embedded with user-defined subroutine can realize the real contact nonlinear characteristics of rough surfaces, and the normal elasto-plastic contact characteristics were well simulated. A comparison between FEA considering normal contact subroutine and tests regarding shear displacement and shear stress is presented in Fig. 18. It was found that a good consistency between the numerically simulated results and test results, and the overall behavior of the wood contact surfaces was well simulated. As shown in Fig. 18, the shear stress rapidly increased with the increase of shear displacement in the initial stage, then gradually tended to moderate, ultimately, there were no obvious horizontal segment, which is due to the contact area being variable and elastic slip deformation was also introduced, so the actual shear stress-shear displacement curve did not show the two-stage behavior (in general, the contact behavior consists of 2 stages: (1) rapid increment of shear stress value without increment of displacement; (2) constant shear stress value with increment of displacement after first stage). It can also be clearly observed that the shear stress-shear displacement relation curve obtained from FEA showed strong nonlinear contact characteristics due to the Coulomb model was used in the process of numerical simulation. Based on the above analysis, it can be obtained that the UINTER normal contact subroutine, embedded in the finite element model (FEM) of wood rough contact surfaces, can better describe the normal contact characteristics of wood rough surfaces, and the proposed normal contact model was suitable for analyzing the problem of wood normal contact surfaces.

Application example-FE model of straight tenon joint
Based on the UNITER normal elasto-plastic contact subroutine model verified in "Validation of model" section. In this section, it was applied to the FE simulation of the straight tenon joint in the historical timber structures, by means of comparing with the experimental results of the straight tenon joint under low-cyclic reversed loading, the universal applicability of UNITER normal contact subroutine in numerical simulation of M&T joints was further verified.

FE simulation of the straight tenon joint
The test setup of straight tenon joint is shown in Fig. 19, the specific details of the test for straight tenon joint were Table 3 The parameters used for normal contact model subroutine g crit is the abbreviation of the elastic slip deformation [45]; Props stands for the abbreviation of material properties in the normal contact model subroutine; Props (1)- (8) Fig. 18 Comparison between experimental and numerical shear stress-shear displacement response of wood normal contact surfaces presented in [37] and [38] and will not be repeated here. The detailed construction and geometry size of the FEM are illustrated in Fig. 20. Similar to the FE simulation of direct shear tests of wood conducted in "Finite element implementation" section, a local coordinate system was also defined to characterize the orthotropic properties of wood. The constitutive model of wood selected is shown in Fig. 16. Material properties used for straight tenon joint model are listed in Table 1. The C3D8R-type element was utilized. Considering the calculation cost, the grid size was set as 30 mm and the meshes consisted of 7144 elements. Typical FE meshes are presented in Fig. 21.
The boundary condition and the loading system were consistent with those of the tested specimens, where the bottom of the column was fully fixed, two horizontal displacements and three rotational degrees of freedom were constrained on the top surface of the column, a vertical load of 54 kN was exerted on the top of the column simultaneously. The loading point was set at the end of the beam to simulate the loading of the actuator.
The established corresponding FE model, with boundary and loading conditions, is illustrated in Fig. 21. The low-cyclic reversed loading under the control of displacement was adopted. During the loading process, firstly, the amplitudes of ± 0.5 mm, ± 1 mm, ± 2 mm and ± 4 mm were utilized for one cycle, respectively, and then the amplitudes of ± 8 mm, ± 16 mm, ± 24 mm , ± 32 mm and ± 48 mm were applied for three cycles, respectively. The loading procedure is illustrated in Fig. 22.
As mentioned in "Determination of model parameter" section, it is worth pointing out, that there are two forms of normal contact between M&T joints, namely, transverse-tangential and tangential-tangential contact. In respect of the contact setting between the contact surfaces of the FEM of straight tenon joint, two different forms of interface contact were, respectively,  defined to simulate the hysteretic behavior and contact status of straight tenon joint, the purpose of which, was to illustrate the applicability and superiority of the established normal contact model by comparison. One was the ordinary contact FEM (abbreviated as OC-FEM), where the "hard contact" was adopted for normal contact between contact surfaces, and the frictional coefficient µ, between tangential contact surfaces, was set as a constant, which was equal to 0.5. The other was the varying contact FEM (abbreviated as VC-FEM) considering the normal elasto-plastic contact between the contact surfaces, in other words, embed the developed and validated UNITER normal contact model subroutine into ABAQUS. Coulomb model was still utilized for tangential contact of contact surfaces, and the frictional coefficient µ was set as 0.5. The normal contact model parameters adopted are shown in Table 4.

Failure patterns
Under the low-cyclic reversed loading, the typical failure modes of the VC-FEM of straight tenon joint, such as the pulling-out of tenon from mortise (Fig. 23a) and the irreversible extruded plastic deformation of mortise and tenon neck (Fig. 23b, c), are shown in Fig. 23. The failure patterns of the straight tenon joint, obtained from the test performed by Zheng [37], are illustrated in Fig. 24, it can be seen that the failure behaviors of the straight tenon joint were consistent with that of VC-FEM, which indicated that the VC-FEM can accurately predict the failure behaviors of straight tenon joints.

Hysteretic curves
The hysteretic curves of straight tenon joint obtained from numerical simulations and experiment are compared in Fig. 25. In general, it was found that the hysteresis curves of the experiment and FEM all exhibited anti-"S" shape, and the pinching effects were obvious, demonstrating the successive slippage features between the contact surfaces of joint. In order to distinguish more clearly the comparison results between VC-FEM, OC-FEM and the test, the hysteretic loops of those were analyzed under the loading amplitude of 48 mm, as illustrated in Fig. 26. It is clear, that numerically simulated  Table 4 Parameter used for wood contact subroutine g crit is the abbreviation of the elastic slip deformation [45]; Props stands for the abbreviation of property in the normal contact model subroutine; Props (1)- (8) are the 8 parameters in the normal contact model subroutine, respectively

Parameter
Symbol Transversetangential contact results of VC-FEM were more consistent with the test results in respect of the hysteresis loop shape, the slippage features and peak moment, and hysteresis loop of VC-FEM was fuller, compared to that of OC-FEM. In addition, a comparison of the experimental and simulated hysteresis loops at the first four small displacement cyclic amplitudes was obtained, as shown in Fig. 27, it is clear that there was a big difference between the simulated and experimental results when the displacement amplitude was less than 4 mm, the distinctions between the two FEMs were very close, the above differences can be explained as: the clearance between the mortise and the tenon was too small due to the manufacturing and installation errors of the tested model, making the initial slope of the hysteresis loops obtained from the test larger than that of the FEM.  Figure 28 displays the skeleton curves for the straight tenon joint, obtained from the FE simulations and the experimental work. It is evident that the change trend of envelop curves, between FEM and the test, was consistent as a whole, in which the bending capacity of the joint gradually increased as the increase of the rotation angle, the growth trend of curve tended to be flat finally, exhibiting an obvious yielding behavior. It is also clear that the skeleton curve of FEM showed obvious symmetry, however, the bending behavior of the tested joint during negative loading were larger than those in the positive loading process, which can be interpreted as: it is difficult to guarantee the complete symmetry of the loading boundary conditions during the loading process; in addition, the initial imperfection, such as wood knots and local cracks, and fabrication errors existed in the tested model. The results of yielding moment, ultimate moment and bending stiffness, obtained from Fig. 28, are summarized in Table 5, in which only the results of the negative loading cycles were showed for simplicity. As shown in Table 5, in terms of yielding moment, ultimate moment and bending stiffness, compared with that of OC-FEM, the simulation results of VC-FEM were closer to the experimental results and the errors were smaller than those of OC-FEM, which were 7.8%, 10.4% and 9.5%, respectively, which were far less than the errors of OC-FEM. This is mainly due to the elasto-plastic deformation of the contact surfaces was considered in the UNITER subroutine, the mechanical mechanism of FEM of the straight tenon joint was more consistent with those of the tested straight tenon joint.

Stiffness degradation curves
The stiffness degeneration curves of the joint obtained from FEM and test are compared in Fig. 29, in which the stiffness of the joint was considered as the mean values of the absolute values of the forward and negative stiffness for the sake of eliminating the error caused by the manufacturing process of the joint and the asymmetry of loading. It can be found that the stiffness of the joint in both test and FEM degenerated gradually with the increase of rotation angle. The stiffness dropped drastically with the rotation angle increasing from 0.001 rad to 0.016 rad. And when the rotation angle was larger than 0.016 rad, the stiffness degradation tended to be gentle. It is clear that the initial stiffness obtained from OC-FEM and VC-FEM was smaller than those of the tested joint, which can be explained as: due to the manufacturing error, the initial clearance between the tenon and the mortise was too small in the test, which led to the larger initial stiffness.  Table 5 Comparisons of yielding moment, ultimate moment and bending stiffness between FEM and test M yield is the yielding moment of the negative loading cycles, and is calculated by the equivalent energy elasto-plastic (EEEP) method [8,46,47]; M u represents ultimate moment of the negative loading cycles; K stands for the bending stiffness, which is equal to the ratio of ultimate moment to ultimate rotation angle

Fig. 29 Stiffness degradation curves
However, when the rotation angle greater than 0.016 rad, compared with that of OC-FEM, the stiffness degradation tendency of VC-FEM was more in line with that of the tested joint, showing VC-FEM can better reflect the plastic deformation behavior of joint under larger displacement amplitude.

Energy dissipation curves
The envelop area of hysteretic loops reveals the amount of energy dissipated [9]. Figure 30 shows the comparisons of the energy dissipation curves between FEM and the experiment. As illustrated in Fig. 30, it was found that the cumulative energy dissipation obtained from VC-FEM was more consistent with that of the experimental joint when the rotation angle ranging from 0.001 rad to 0.048 rad, compared to that in OC-FEM. When the rotation angle changed in the range of 0.048 rad to 0.096 rad, the cumulative energy dissipation obtained from FE simulation was larger than the experimental value. The main reasons of above phenomena are as follows: firstly, the assumption of the contact state adopted in the normal elasto-plastic contact theory was still different from actual contact state of the contact surface of M&T joint during the experiment, there were still great difficulties in reflecting the complicated contact characteristics of the contact surfaces of M&T subjected to cyclic loading, truly and comprehensively. Secondly, the fabrication and installation error existed in tested joint, having an influence on experimental results. Thirdly, the cyclic loading caused irreversible plasticity deformation of the joint in the test, which made the contact surfaces become more and more smooth and reduced energy consumption under large loading amplitude.

Comparison between VC-FEM and OC-FEM
According to the analysis in "Comparison between FE simulation results and experimental results" section, it can be clearly noted that the simulation results of the VC-FEM embedded in the UINTER subroutine were more consistent with the experimental results. In order to further make a clear distinction between VC-FEM and OC-FEM, the contact behaviors of VC-FEM and OC-FEM of straight tenon joint, in terms of contact state, normal contact stress and normal contact force of rough contact surfaces, were analyzed, respectively.   Figures 31 and 32, respectively, compare the interface contact state for reverse loading and reverse unloading of the two FEMs at a displacement amplitude of 48 mm, in which different nephogram colors indicate different contact states: red represents "sticking" state, which means that the surfaces of the tenon and the mortise are in contact, but there is no slippage between the contact surfaces; green represents "sliding" state; blue stands for "open" state.

Analysis of contact state
It can be clearly noted that the contact state between VC-FEM and OC-FEM was quite different under the same loading amplitude. Compared with OC-FEM, VC-FEM of straight tenon joint had more sticking and slipping phenomena in the same region. This is because that the normal elasto-plastic contact properties of rough contact surfaces were considered in VC-FEM, which was more conducive to the energy dissipation of the joint. This was also more consistent with the plump hysteresis loop of VC-FEM observed from Fig. 26, namely, the larger the hysteresis loop area, the greater the energy consumption.

Normal contact stress analysis of contact surface
Under the displacement amplitude of 48 mm, the comparison of the normal contact compressive stress between two simulation methods is shown in Fig. 33. It can be seen that there will be extrusion and friction between the end and neck of tenon and the edge of mortise when the joint was subjected to the low-cyclic reversed loading, resulting in obvious stress concentration phenomenon. This failure mode was consistent with those observed in the tests presented in [37] and [38]. It can be also found that the stress distribution and stress value of the two models, under the same displacement amplitude, were different in the same region, and the difference of maximum contact compressive stress between VC-FEM and OC-FEM was 15.5%, which must be considered for the sake of more exactly simulating hysteresis behavior and normal contact characteristics of M&T joints. This implies that VC-FEM of M&T joint can better reflect and exactly simulate the real normal contact behaviors of rough contact surfaces.

Normal contact force analysis of contact surfaces
The normal contact forces of the upper and lower contact surfaces (m1-s1, m2-s2, see Fig. 3) of tenon and mortise were selected for analysis, in which the normal contact force is equal to the product of the contact area and the contact stress for all units on the contact surface of the model. A comparison of contact force for the normal contact surfaces between the two simulation methods is shown in Fig. 34. It can be seen that there was little difference for the contact forces between the two models in the early stage of loading, due to the model is in the elastic state; at the late stage of loading, the normal contact force of VC-FEM was gradually greater than that of OC-FEM, which is due to the mutual extrusion of the tenon and mortise, causing the irreversible plastic deformation of contact surfaces, and the embedded compressive behavior between the contact surfaces was more obvious. In addition, a comparison of the frictional energy dissipation of the two simulation methods was obtained, as shown in Fig. 35. It was found that the frictional energy dissipation of VC-FEM was greater than that of OC-FEM at the late stage of loading, which is due to that the tangential friction between the contact surfaces increased with the increase of the normal contact force, making the frictional energy dissipation of the joint increased as well. This is consistent with the phenomena shown in Figs. 26, 31, 32, and it can be concluded that VC-FEM can dissipate more energy through frictional slipping, compared with OC-FEM, especially in the later stage of loading, which was more in line with the experimental results and better reflected the frictional energy-consuming mechanism of M&T joints.

Conclusions
This paper investigated the normal contact performance of M&T joints in traditional timber structures, a normal elasto-plastic contact model of rough contact surfaces was proposed and derived. The normal contact model subroutine of M&T joints was compiled by FORTRAN, which can be implemented into ABAQUS through UINTER interface. The shear tests of transverse-tangential and tangential-tangential contact surfaces and the corresponding FE simulation embedded wood normal elasto-plastic contact model were carried out, respectively. The VC-FEM and OC-FEM of the straight tenon joint subjected to cyclic loading were developed, respectively. The following conclusions can be obtained: (1) The proposed normal elasto-plastic contact model of rough contact surfaces can be utilized to characterize the normal contact behavior of the rough surface of M&T joints, and describe the relationships between the actual contact area, the normal contact load and wood mechanical properties and surface micro-morphology parameters. (2) The wood property parameters and surface micromorphology parameters of the normal elastoplastic contact model of the rough surfaces were determined. The variation laws of wood micromorphology for different surface grain directions under different roughness were analyzed. The results indicated that the cross section was rougher than the tangential section and radial section. (3) The user-defined normal elasto-plastic contact subroutine was validated against the results of the shear tests of wood rough surfaces, showing that the test results agreed well with the simulation results, and the normal elasto-plastic contact model could well describe the normal nonlinear contact characteristics between rough contact surfaces. (4) The universal applicability of UNITER normal elasto-plastic contact subroutine of wood contact surfaces was verified by the comparison between the simulation results of VC-FEM and OC-FEM and the experimental results. The results of VC-FEM for straight tenon joint were more consistent with the experimental results in terms of hysteresis curve, envelop curve, stiffness degradation and energy dissipation behavior, compared to that of OC-FEM, which also showed that the correctness and applicability of the normal elasto-plastic contact model to simulate the hysteretic performance of M&T joints accurately. The VC-FEM can be used to predict the failure behaviors of M&T joints. (5) The contact behavior of VC-FEM of straight tenon joint, in terms of contact state and normal contact stress and contact force of normal rough contact surfaces, had a difference in comparison with those of OC-FEM. The VC-FEM of straight tenon joint had more sticking and slipping phenomena and good energy dissipation capacity. The maximum contact compressive stress of VC-FEM was 15.5% greater than that of OC-FEM. The VC-FEM can dissipate more energy through frictional slipping of contact surfaces, compared to OC-FEM. The VC-FEM of M&T joint was capable of simulating the real normal contact behavior of rough contact surfaces exactly.
Due to limited resources, the material properties (e.g., new timber) utilized in our present research will differ somewhat from those (e.g., old wood) of traditional M&T joints. In order to reduce this difference, future studies can assess the material properties of traditional joints as closely as possible by artificially simulating material behavior degradation (e.g., establishing the degradation model of wood material properties to obtain the relationship between the properties of ancient wood materials and the properties of new wood materials). Finally, based on the normal contact constitutive model of wood rough surfaces proposed in this paper, future studies will use the degraded material properties to more realistically evaluate the normal contact properties of the ancient wood material.