Experimental and Numerical Studies on Repaired Wooden Beam of Traditional Buildings in Huizhou Region, China

: In this study, three different repair methods are proposed for the common broken parts of traditional wooden beams. Four wooden beams without initial damage and six repaired beams in a planar wooden frame are tested and numerical simulated. The test phenomena, bearing capacity, and stiffness of all the wooden beams are analyzed. Compared with the wooden beams without initial damage, the bearing capacity of the beams with upper inlay repair, upper core repair, partial tenon repair, and lower core repair increases by 38.93%, 13.06%, 5.08%, and 3.94%, respectively. Furthermore, the experimental and numerical results indicate that the upper and lower parts of the wooden beam with initial damage can be more effectively repaired by the inlay and core methods, respectively. When the tenon is partially damaged, the mechanical properties of the wooden beam are improved after repair. The simulation of lower inlay repair shows that the bearing capacity of the wooden beam is inversely proportional to the repair height and the distance between the repair position and the span. Based on the experimental results, a dovetail form of inlay repair is proposed, and it is numerically veriﬁed that this method can effectively reduce the stress concentration.


Introduction
As an environment-friendly and sustainable construction material with excellent thermal, acoustic, electrical, mechanical, and aesthetic properties, wood has been commonly used in buildings since ancient times [1]. Over the past centuries, numerous ancient buildings have been seriously damaged due to natural and man-made factors. The damage to the wooden beam can impact the overall stability of timber structures. Notably, some ancient buildings are of great cultural and artistic value to human beings [2]. Therefore, the maintenance of such buildings is extremely important, and it is essential to establish effective repair methods for the damaged wooden beams and evaluate their mechanical properties after repair [3]. For example, as shown in Figure 1, Dunmu hall, a damaged traditional building in the Huizhou area, has been repaired. It can be seen from the Figure 2 that the timber structural forms of China and Europe are different. In the Huizhou region, moon patterns are carved on the surface of wooden beams. And the sparrow brace which is at the junction of wooden beam and column, is one of the characteristic decorative features of ancient Chinese architecture.
Several experimental, theoretical, and numerical studies have been conducted on reinforced wooden beams. Considering the superior performance of fiber materials, some researchers have examined the flexural performance of wooden beams with fiber-reinforced polymer (FRP) and carbon-fiber-reinforced polymer (CFRP) materials [4][5][6][7][8][9][10][11]. For example, Bhat [12] carried out two-point loading tests on different wooden beams reinforced with CFRPs. The results indicated that the CRFP material had the strongest influence on the flexural stiffness of all the reinforced wooden beams when the percentage of CFRP area was greater than 0.50%. Rescalvo et al. [13] summarized the different layouts available for FRP fabrics and pultruded laminates based on the previous studies. Jose et al. and Shi et al. [14,15] investigated the flexural performance and stiffness of wooden beams with glass FRP reinforcement. Further, the reinforcement of wooden beams has been investigated theoretically. Yang et al. [16] derived the relationship between bending moment, axial force, and strain through the micro-beam element model. Combined with the experimental tests, the expressions for the load bearing capacity in the tensile zone of wooden beams reinforced with FRP were derived. Guo et al. [17] proposed a numerical model for tenon joints reinforced by self-tapping screws. They studied the effects of different bolt diameters and yield strengths on the damage modes, ductility, and moment-resisting capacity. Tsalkatidis [18] developed a numerical model for glued laminated timber beams using ANSYS software and compared the flexural performances of steel plate reinforced timber beams and unreinforced timber beams. Baño [19] established finite element models for defective wooden beams and predicted their maximum load under bending. Kotšmíd et al. [20] performed buckling tests on three square beams with different lengths based on the relationship between measured axial displacement and axial load.  Combined structures have been proposed as an effective approach to enhance the loadbearing performance of wooden beams. Fajdiga et al. [21] presented a numerical model for a metal-reinforced wood composite. The model consisted of several interconnected components made of different base materials. It was used to estimate the design bending stiffness and strength of a certain composite material. Furthermore, some researchers have studied the influence of steel strength, shape, and size on the flexural performance of timber-steel composite beams [22][23][24]. Corradi et al. [25] presented a construction method for stainless steel reinforced timber with some engineering examples. Chen et al. [26] experimentally investigated the flexural performance of wood-bamboo composite I-beams to enhance the structural performance of I-beams. Duan et al. [27] proposed a steel-wood combination beam with H-shaped steel beam webs pasted onto timber, and the flexural load capacities of three hot-rolled H-arch combination beams and one pure steel beam were tested. Inspired by the natural bamboo structures with high stiffness and stability when subjected to lateral loads, Zhang et al. [28] conducted flexural tests on a combined structure with hollow sections and internal stiffener ribs. Lukin et al. [29] numerically analyzed the characteristics of wood composite floor beams. The results revealed that compared to the conventional wooden beams, the deformation capacity of glued laminated beams was reduced.
However, some studies on the reinforcement of wooden beams focused on the combination of different materials (steel, bamboo, etc.), which destroyed the original appearance of traditional buildings to some extent. This is not in line with the fundamental principle of conservation for traditional buildings: "repair the old as before" [30]. To address the above issues and maintain the sustainability of traditional buildings, the Scientific Committee of the COST Association has carried out extensive research on the reinforcement of timber structures in Europe from 2011 to 2015. Franke et al. [31,32] summarized several forms of wooden beams and general failure modes, illustrating possible reinforcement techniques for the wooden beams. In terms of reinforcement technology, reinforcing elements such as rods, bars and plates can be used in combination with the original structure. Dietsch et al. [33] applied self-tapping screws and threads as reinforcement for wooden structures in the European region, illustrating their wide range of prospective applications. Schober et al. [34] investigated the timber structures reinforced by FRP. The results verified that FRP is an economical and effective reinforcement for timber structures. Unfortunately, in China, the efficacy of replacement and reinforcement methods for the damaged parts of practical engineering structures has been rarely examined.
Although the Scientific Committee of COST Association has proposed some methods in 2015, the effectiveness of these methods for repairing wooden beams in the Huizhou region of China has not been tested and numerically simulated. Besides, the extensive research was carried out in Europe and America, but the wooden beams are connected with other members to form a wooden frame in practical applications. And the entire structure is jointly stressed, not a single beam. These kinds of simplification usually lack both the consideration of the geometric features of the beams, as well as of the appropriate numerical verification. Therefore, it is more reasonable to study the mechanical properties of wooden beams in timber frames. Finally, there are few studies on the establishment of finite element models for the wooden beams in the ancient buildings of Huizhou region [35][36][37].
To this end, the goal of this study is to evaluate the effectiveness of various repair methods in improving the mechanical properties of wooden beams in a planar timber frame in the Huizhou region of China. Based on different engineering cases, five sets of two-span timber frames, including six repaired wooden beams and four wooden beams without initial damage (as comparison specimens), are tested, and the accuracy of the results is numerically verified. The failure modes, flexural capacity, and stiffness of the wooden beams are investigated. In addition, the influence of different parameters on the wooden beam repaired by lower inlay method is analyzed. Finally, a new form of inlay repair is proposed and verified.

Establishment of Repair Methods
Based on a site survey of the damaged parts of wooden beams in the traditional building Dunmu Hall, three repair methods have been designed to stably connect the original structure and replacement part. These repair methods are shown in Figure 3. In the core repair method, bolts are used to connect the original structure with the replacement part. In the inlay repair method, the original structure and replacement part are bonded with an adhesive layer. In the tenon repair method, the wooden beam is connected to the completely repaired tenon with steel plates and bolts and the 1/3 repaired tenon with bolts.
These methods are designed based on the principle of "maintaining original appearance after restoration". A Q235 steel plate with a yield strength of 235 MPa and thickness of 10 mm is used. The length of the M8 bolt is 190 mm. All the plates and bolts are placed inside the structure to retain the original appearance of the old building. The wooden beams are divided into five groups: 1-5. Among them, two groups (group 1 and 2) are without initial damage and the other three (group 3-5) are repaired by different methods. The letters Y and Z refer to the right and left side of the timber frame respectively.

Experimental Design
With reference to the traditional building Dunmu Hall in Xu Village, Huizhou region, five 1:2.75-scale wooden frames have been designed for the vertical loading test. The details were presented in Table 1. Figures 4 and 5 show the assembly process and overview of wooden beams after repair, respectively. The planar wooden frame has been assembled from wooden beams, wooden columns, and sparrow brace. The sparrow brace in Figure 4c is a load-bearing and decorative component. A rectangular block named "Guanjian" is used to connect the beam to the column. The dimensions of the wooden beams are 2340 mm × 210 mm × 200 mm, and their surface is chamfered with a radius of 6 mm along the cross-sectional direction. The radius and height of the wooden column are 90 and 180 mm, respectively.

Mechanical Characteristics of Materials
The wooden beams were made of Malaysian Merbau and Chinese fir. Figure 6 shows the test blocks prepared to examine the mechanical properties of wood based on the relevant standards [38][39][40]. The mechanical properties of wood are listed in Table 2.

Experimental and Loading Devices
The experiments were performed in the structural test hall at the Anhui University of Architecture. A 50 t hydraulic jack was used to apply vertical load to the distribution beam and transfer it to the wooden beam. To prevent out-of-plane displacement of the wooden frame, lateral support restraint was set. Considering the effect of actual constant load of roof on the mechanical properties of the wooden frame, the wooden column was loaded to 20 kN in the first stage of loading, which was kept constant during the test. For failure loading, a preliminary value of 20 kN was maintained for 5 min to eliminate the installation error of the specimen. When the test was officially started, the vertical load was applied at a rate of 2 kN/min until beam failure. When the load value dropped to 80% of the ultimate load, it was considered that the specimen is destroyed. The loading test setup is shown in Figure 7.

Failure Modes
The failure modes of the B0 and B1 groups are shown in Figure 8a,b. The wooden beam moves downward with the increase in load during the initial loading process. When the load reaches 47% of the final load, transverse cracks appear in the compression zone along the down-grain direction of the beam. When the load reached 80% of the final load, the pulling and cracking sound of wooden fibers started to become intense. Eventually, the tensile zone at the bottom one-third and two-third of the length cracked with a loud sound. Figure 8c,d present the failure modes of the B2 group. When the load of the specimen B2Z reaches nearly 45% of the ultimate load, a preliminary crack is developed in the upper-pressure zone of the specimen. As the load increases, folded wood fibers appear at the pressure area. When the load reaches nearly 80% of the maximum bearing capacity, transverse cracks appear at the tensile zone along the left one-third of the specimen. At the same time, a horizontal crack at one-third height of the span gradually expands to the ends. Further, the specimen cracks at the span center. Eventually, the crack extends toward the bottom of the tensile zone of the specimen and the tensile zone suddenly breaks. The reinforced part of the B2Y specimen is completely peeled off from the structure.
Figure 8e,f demonstrate the failure mode of the B3 group. Under loading, the two specimens of this group exhibit significantly different damage patterns. For specimen B3Z, extrusion wrinkling gradually appears in the repaired part of the pressure zone when the load reaches approximately 63% of the maximum bearing capacity. As the load increases, it extends along the edge of the repaired part toward the loading points on both the sides. For the B3Y specimen, the initial failure occurs at the repaired part of the vertical contact with the structure. This phenomenon is primarily attributed to the considerable transverse stress at this location and the initial formation of transverse cracks. Finally, the wooden fibers gradually pull apart laterally along the top of the repaired section. Figure 8g,h show the failure mode of the B4 group. For specimen B4Z, as the load increases, the repaired part of the pressure area is horizontally extruded from the structure, and a large gap is formed between the reinforced contact area and the structure. When the load reaches nearly 58% of the maximum bearing capacity, a clattering sound continuously comes from the end of specimen B4Y. Eventually, the reinforced part splits along the bending direction of the bolts, and the bottom of the tensioned area is detached from the original structure.

Force-Displacement Curves
The experimentally obtained mechanical properties of the test groups are compared in Table 3, where F u represents the ultimate load capacity, D u represents the maximum displacement, and K represents the stiffness of the specimen. The stiffness is the slope of the curve between 10% and 40% of the ultimate load [41]. According to Table 3, the wooden beams with lower inlay repair exhibit reduced ultimate load capacity (−38.75%). As the load increases, the stress between the repaired part and the bottom of the wooden beam increases, which leads to a horizontal crack corresponding to adhesive degumming. The wooden beams with complete tenon repair exhibit significantly lower ultimate load capacity (−65.37%) than the undamaged beams because the bolts bend under tension and squeeze the wood, and then the repaired part is gradually separated from the wooden beam.  Figure 9 shows a comparison between the force-displacement curves of the five specimen groups. The ultimate load and ultimate displacement of each group of specimens are not equal, but the load and displacement linearly increase at the beginning of loading. In the wooden beam without initial damage, the ultimate load capacity of group 2 is nearly 3.71 times higher than that of group 1. Due to the large wood dispersion, the wooden fibers of the unrepaired beam are easily pulled off, resulting in smaller plastic deformation of specimens B1Y and B1Z. However, the variation trends of stiffness and ultimate bearing capacity of both the specimens are approximately similar. In the group 3, the specimen B2Y is less stiff than the specimen B2Z, which is due to the gap formed between the reinforced part and the wooden beam at the bottom of specimen B2Y as the load increases. When the specimen enters the yielding stage, the reinforcement contact area at the tensile zone is co-tensioned and the gap gradually increases. Eventually, the specimen peels off laterally along the core-planted reinforcement part, and brittle damage occurs. The ultimate bearing capacity of specimen B2Y is increased by 3.96%. For the bending specimens in group 4, the ultimate bearing capacity of specimen B3Y is only 0.41 times that of specimen B3Z. When the bottom tensile stress of specimen B3Y is large, stress concentration occurs in the reinforcement part, inducing a transverse tensile gap after gradual detachment. The reinforcement in the compression zone is better than that in the tensile zone because the plastic deformation capacity is enhanced after the replacement and reinforcement of damaged parts in the compression zone. However, because stress concentration occurs in the reinforced part of the tensile zone, it is more likely to cause brittle damage. For the B4 group, the ultimate load of specimen B4Z is 0.3 times that of specimen B4Y, while the displacement of specimen B4Z is 0.73 times that of specimen B4Y. This is mainly because with the increase in the load, cracks are formed in the contact area between the screw and beam in the repair section (as shown in Figure 8g), which prevent the repaired section from being stressed together with the wooden beam.  Figure 9f compares the load-span displacement curves of each group of specimens. It is clear that the curve can be divided into of two stages: elastic and plastic stages. When the span displacement of each group of repaired specimens made of Merbau is about 10 mm, the specimens enter the plastic stage from the elastic stage. The stiffness of the reinforced specimen B3Z in the plastic stage is larger than that of the B2 and B4 groups specimens. Within the load range of 152-225 kN, the stiffness of repaired specimen B4Y is smaller than that of group B1 specimens, but the deformation of specimen B4Y is improved to some extent. Compared with the group B1 specimens, the bearing capacity of repaired specimens B2Z, B3Z, and B4Y increase by 13.06%, 38.93%, and 5.08%, respectively.

Force-Strain Curves
Taking specimens B1Z, B3Z and B3Y as examples, seven measurement points are set at the mid-span section height of the wooden beams, and four points are set along the one-third and two-third span section heights of the timber beam (Figure 10a). The span of the inlay repair group contains three measurement points, and additional measurement points (4, 5, 6, 7, and 8) are on the repaired area. The specific arrangement of measurement points is shown in Figure 10b,c. According to the measured strain at each measurement point with equally spaced distribution in the span section of the specimen, the maximum tensile stress occurs at one-third and two-third height of the specimen at the bottom. It can be seen in Figure 11 that the force-strain curves of the specimen are almost linear. The maximum compressive and tensile strains of the two groups of specimens at one-third and two-third height are larger than that at the corresponding height in the span. For the specimen B3Z, when the load value is up to 260.32 kN, the reinforcement part of the compression zone exhibits extrusion folding. As shown in Figure 11g, the strain value of measurement points 8 and 4 increases at a higher rate, and the maximum compressive strain is 0.00081. The strain values of the measurement points 5 and 7 arranged at the bottom of the reinforced part change consistently. The strain variation is relatively stable, and the maximum compressive strain is 0.0003. This suggests that the stress value on both sides of the upper inlay part changes significantly with the vertical contact part of the original structure, which should be given more attention in the analysis. By contrast, for the specimen B3Y, an increase in the load has no effect on the strain at measurement points 4 and 8 on either end. The rate of strain change at the measurement point 6 is 44.4% of that at the points 5 and 7, but the strain value is 0.0005 greater than at the other measurement points, as shown in Figure 11 h.

Mid-Span Height-Strain Curves
Taking specimens B0Z, B1Z, B2Z, B3Z, and B4Y as examples, Figure 12 dem the strain variation with load at each measurement point along the mid-spa height of the specimens. It can be seen that the height of each group of specim the mid-span section is linearly related to the wood fiber strain during the ear loading. With the increase in the load, the height of the neutral axis of the speci out initial damage remains constant. Among other groups, the specimens B2Z move toward the top of the wooden beam in the neutral axis due to the incre strain of the tensile zone with the increase in load. Since the wooden fibers in pression zone are folded, the neutral axis of specimen B3Y moves toward the the wooden beam, and elastoplastic deformation occurs.

Mid-Span Height-Strain Curves
Taking specimens B0Z, B1Z, B2Z, B3Z, and B4Y as examples, Figure 12 demonstrates the strain variation with load at each measurement point along the mid-span section height of the specimens. It can be seen that the height of each group of specimens along the midspan section is linearly related to the wood fiber strain during the early stage of loading. With the increase in the load, the height of the neutral axis of the specimen without initial damage remains constant. Among other groups, the specimens B2Z and B4Y move toward the top of the wooden beam in the neutral axis due to the increase in the strain of the tensile zone with the increase in load. Since the wooden fibers in the compression zone are folded, the neutral axis of specimen B3Y moves toward the bottom of the wooden beam, and elastoplastic deformation occurs.

Finite Element Analysis
To effectively reflect the loading mechanism of wooden beams in the wooden frames, the mechanical performance of different specimens was compared through parametric analysis. The simulations were conducted using the ABAQUS numerical analysis software.

Repair Methods
For the partial damage of wooden beams, the three repair methods were numerically examined with the same structural dimensions as the test specimens. The surface of the wooden beam had a rounded chamfer with a radius of 6 mm in the section direction. The replacement part was of the same size as the damaged part. The contact surfaces between the plate and bolt were set to "tie" constraint. The coefficients of friction of wood to wood and steel were 0.35 and 0.3, respectively [42]. The stiffness constant of adhesive was 200 N/mm 3 . Figure 13 shows the numerical models of the three different repair methods.

Finite Element Model
The structure and dimensions of the five groups of single-story, two-span timber frame models are the same as those of the test model. The wooden column model has a height and radius of 1800 and 90 mm, respectively. The bottom end of the column is set to a fixed restraint with the stone foundation. Figure 14 presents the mesh generation scheme of the planar wooden frame. An elastic-plastic model of wood anisotropy is defined by using the nine elastic engineering constants (as shown in Table 4) and the generalized Hill yielding criterion. Due to the anisotropy of wood, the directions of wooden beam, key, and alternate finch material must be assigned separately. The force surface is set at the top of the wooden column and the third equivalent point of the wooden beam, and the motion coupling constraint is established with the loading point. A constant load of 20 kN is applied to the top of three wooden columns in "Step1," and the wooden beam load is applied in "Step2". To consider the adhesion between the reinforced part in specimen B3Z and the original structure, an adhesive contact is used to connect the two parts. Considering the slip effect, the wood-wood and wood-steel plate contact surfaces are set to hard contact. The connection of Guanjian to the wooden beam and column is set as "tie" constraint. Eight-node linear hexahedral cells (C3D8R) are employed in the simulated wooden beam, spar, and key, and four-node linear tetrahedral cells (C3D4) are adopted in the wooden columns, core reinforcement part, and tenon reinforcement part. To improve the calculation efficiency, the grid division size is set to 20 mm.   Figure 15 shows the "S33" stress cloud parallel to the grain direction of the wooden beam. It can be seen that the stress concentration area at the bottom mainly appears at 1/3 and 2/3 of the length. The maximum compressive stress of group B1 is 63.21 MPa, and the maximum tensile stress is 45.29 MPa. In the group 3, the maximum stress of specimen B2Z at the top of the repaired part is 42.42 MPa, and the maximum stress at the bottom of the repaired part of specimen B2Y is 101.1 MPa. The stress value along the 1/2 section height of specimen B2Y is larger than that of the unreinforced group. For specimen B3Y, the adhesive layer is about to reach the plastic, and the top of the corresponding repaired part is open to the bonding surface of the structure. It can be seen from Figure 14f that stress concentration occurs at the top bonding position. With the increase in the load, horizontal and transverse tensile cracks appear on the bonding surface. For specimen B4Z, the stress at the transverse contact part between the repaired part and the original structure gradually increases. In general, the numerically obtained damage morphology of the specimens is consistent with the experimental results.

Comparison between the Simulation and Test Results
A comparison between the experimental and simulated force-displacement curves is presented in Figure 16. The experimental and numerical results for the ultimate load and stiffness are listed in Table 5. Except for the specimens B3Y and B4Z, the error between the experimental and numerical results for each specimen is within 20%. For specimen B3Y, the tensile fracture strength and shear yield strength of the epoxy resin decrease with time, and this strength reduction problem is not taken into account in the viscous contact setup used in the simulation model. For specimen B4Z, as the load increases, the steel plate gets detached from the original structure, and the repaired part cannot be jointly stressed with the original structure. The comparison results in Table 5 indicate that the performance of the repaired specimens is in the following order: B3Z > B2Z > B4Y > B2Y > B3Y > B4Z. A comparison between test observations and finite element simulation is presented in Figure 17.

Influence of Repair Height in Lower Inlay Method
To investigate the effect of repair height on the load-bearing performance of wooden beams reinforced by the lower inlay method, 1/5, 2/5, and 3/5 of the span section are considered as the repair height. Figure 18 shows that the bearing capacity of wooden beam with 1/5 height is 1.11 times that of the beam with 2/5 height and 1.32 times of beam with 3/5 height. The simulated values of ultimate load under different repair heights are listed in Table 6. It is clear that as the repair height of the wooden beam continues to increase, the decrease in its bearing capacity of the wood beam gradually increases.  To study the effect of repair position on the load-bearing performance of wooden beams reinforced by the lower inlay method, the strengthened part with center distance of 0, 200, and 400 mm are considered. Figure 19 shows the load-displacement curves of the wooden beam with different repair positions. The simulated values of ultimate load under different positions are listed in Table 7. It can be seen that the farther the repair position is from the center of the span, the worse the bending performance of the wooden beam.  According to the test results, a novel form of dovetail inlay named B3N is proposed with reference to the dovetail joints, as shown in Figure 20. The effect of inlay repair form on the load-bearing performance is evaluated by B3N modeling. The stress cloud diagram is shown in Figure 21. Figure 22 shows the force-displacement curves of the wooden beam with different repair forms. The maximum load capacity and maximum stress of specimens B3N and B3Y are compared in Table 8. It can be seen that the stress value of the upper contact part in B3N is smaller compared than that in B3Y, which effectively prevents the stress concentration. The bearing capacity of B3N is 27.09% higher than that of B3Y, which proves the effectiveness of the proposed dovetail inlay form for repair.

Conclusions
The mechanical properties of wooden beams in traditional timber-frame structures (in Huizhou, China) were examined with three different repair methods. The wooden beams with initial damage were repaired by the core repair, inlay repair, and tenon repair methods, and four wooden beams without initial damage were used for comparison. The main results of the study are summarized as follows: 1.
In the flexural test, the bearing capacity of Merbau was 3.73 times higher than that of Chinese fir in wooden beams without initial damage. The wooden beams exhibited brittle damage, and the maximum stress values were observed at the bottom 1/3 and 2/3 height of the tensile zone. The stress value on both sides of the upper inlay part varied greatly with the vertical contact part of the original structure, which should be given more attention in the analysis. For the lower inlay part, the top part of the reinforcement in contact with the original structure should be emphasized.

2.
The experimental and finite element results were combined to compare the bearing capacity and stiffness of wooden beams, and the results revealed the following performance order of the repair methods: complete tenon repair < lower inlay repair < lower 1/2 core repair < 1/3 tenon repair < upper 1/2 core repair < upper inlay repair. It was observed that for engineering applications, the upper and lower parts of the wooden beam with initial damage could be more effectively repaired by the inlay method and core method, respectively. Compared to the overall damage of the tenon, the proposed tenon repair method was more effective for partial tenon damage.

3.
The parametric analysis revealed that the bearing capacity of the beam with lower inlay repair was inversely proportional to the repair height as well as the distance between the repair position and the span. Compared with the rectangular inlay form, the dovetail inlay form provided better performance indicators and effectively prevented stress concentration in the contact area.
Author Contributions: Conceptualization, methodology, investigation, validation, writing-original draft, data curation, software, Y.J.; writing-review and editing, resources, funding acquisition, supervision, Q.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by National Science and Technology Support Program, grant number 2012BAJ08B02.