#### Abstract

This paper describes some finite element models for analyzing the mechanical behavior of adhesive joints. In these models five layers of solid elements were used across the adhesive layer in order to increase the accuracy of the results. The finite elements were refined gradually in steps from adherends to adhesive layer. In these models, most of the adherends and adhesive were modeled using solid brick elements but some solid triangular prism elements were used for a smooth transition. In some of the models, linear interpolation elements of full or reduced integration and of hybrid formulation were used. In other models, quadratic interpolation elements of full or reduced integration and of hybrid formulation were used. Comparisons are drawn between models with different modeling approaches as well as different types of element combinations in order to find a suitable model to predict the behavior of adhesive joints.

#### 1. Introduction

Fastening techniques are used extensively in different industry fields for joining various materials in the assembly of components and structures. Many efforts have been spent to develop sheet material joining techniques for application into light-weight structures [1–3].

Adhesive bonding has many characteristics comparable with conventional mechanical fastening and welding methods used in structural engineering. It also has many exclusive advantages such as low bonding temperature, light weight, high stiffness, and good fatigue resistance. Consequently it is becoming a widespread candidate technique for joining light-weight structural components. A considerable amount of theoretical and experimental research has been carried out on the static and dynamic behavior of adhesive joints (e.g., [4, 5]).

To design structural adhesive joints, it is necessary to be able to analyse them. This means determining the stresses and strains under a given loading and predicting the probable points of failure. When different boundary conditions are considered by a closed-form analysis, the limitation is how tractable a realistic mathematical model is within an algebraic solution. Usually it is necessary to simplify the models to some extent to make analytical solutions feasible. Many studies have been published all with analytical or experimental simplifications that restrict the usefulness of the results. With finite element (FE) techniques, however, the limiting factor is more likely to be computing power. The FE method now commonly used is well suited to the estimation of stresses in joints of almost any geometric shape [6–9].

During the last four decades, many of the existing adhesively bonding processes have been simulated by FE methods. Woole and Carver’s paper [10] was concerned with the stress analysis by FE method of a bonded single-lap adhesive-bonded joint. A modified version of the well-known Wilson stress analysis program was used in the case of plane stress. The authors used 2 elements’ thickness to model the through-the-thickness behavior of the adhesive layer. Stress concentrations as functions of dimensionless, geometric, and material parameters were presented. However, because of the sharp discontinuity between the mechanical properties of the adherend and the adhesive, the use of 2 elements is not sufficient. Smooth transitions between the adherends and the adhesive are necessary in order to obtain accurate results. In later work by Adams and Peppiatt [11], stresses in a standard metal-to-metal adhesive-bonded lap joint were analyzed using a two-dimensional FE method and comparisons were made with previous analyses. In the paper, particular attention was paid to the stresses at the ends of the adhesive layer. Unlike previous work, which assumes the adhesive to have a square edge, the adhesive spew was treated as a triangular fillet. The results show that the highest stresses exist at the adherend corner within the spew. This model is closer to realistic adhesive joints. Carpenter and Barsoum [12] modeled the adherends as two-node beam elements and the adhesive layer as a linear plane element with offset nodes. The number of degrees of freedom is reduced appreciably by this approach because the adherends and the adhesive use the same node.

Anyfantis and Tsouvalis’s work [13] was focused on the numerical simulation of single-lap bonded joints, based on cohesive zone modelling techniques. The models were built in a 3D FE space. The adherends were modelled with continuum elements whereas the entire adhesive layer was modelled with cohesive elements. A mixed-mode cohesive model was used as the constitutive relationship between the cohesive elements. The traction increase part of the cohesive laws is given by an exponential function, which describes the elastoplastic adhesive response, and the traction decrease part is given by a linear function, which describes the initiation and propagation of damage. By using this model, it was possible to calculate the developed peel and in-plane and out-of-plane shear stresses over the adhesive area. Hybrid-adhesive joints are an alternative technique for stress reduction in adhesive joints. The joints have two types of adhesives in the overlap region. The stiff adhesive should be located in the middle and the flexible adhesive at the ends. The effect of the hybrid-adhesive bond line on the shear and peeling stresses of a double lap joint was investigated by Özer and Öz [14]. A 3D FE model of the double lap joint has been created based on solid and contact elements. The contact problem was considered by modelling the interface as two surfaces belonging to adherend and adhesive. The results show that the stress components can be optimized using appropriate bond-length ratios.

In the case of analysis of adhesive joints, the thickness of adhesive is much smaller than that of the adherends. FE meshes must accommodate both the small dimension of the adhesive thickness and the larger dimension of the remainder of the whole model. Moreover, the failures of adhesive joints usually occur inside the adhesive layer. In other words, the strength of adhesives is usually lower than that of adherends. It is thus essential to model the adhesive layer by a FE mesh which is smaller than the adhesive thickness. The result is that the FE mesh must be several orders of magnitude more refined in a very small region than is needed in the rest of the joint. It is also important that a smooth transition between the adherends and the adhesive be provided. To determine the physical nature of adhesive joints, many researchers have limited their investigations to single-lap joints because they involve relatively simple and convenient test geometries. However, most other joints may be obtained through some combination or repetition of this basic type.

This paper describes some FE models for analyzing the behavior of single-lap adhesive joints. To overcome the limitations described above, five layers of solid elements were used across the 0.05 mm thick adhesive layer. The main objective of this treatment was to increase the accuracy of the results. The FE models were refined gradually in steps from adherends to adhesive layer. Most of the adherends and the adhesive were modeled using solid brick elements but some solid triangular prism elements were used for a smooth transition. Comparisons are drawn between models with different modeling approaches as well as different types of element combinations in order to find a suitable model to predict the behavior of adhesive-bonded single-lap joints.

#### 2. Configurations and Material Properties

The single-lap adhesive joint studied in the present paper includes the lower adherend, adhesive layer, and upper adherend, as shown in Figure 1. The two adherends used were 2024-T3 aluminium alloy plates of dimensions 200 mm long × 25 mm wide × 4 mm thickness. The elastic material constants of the adherends were as follows: Poisson’s ratio and Young’s modulus GPa. The elastic material constants of the adhesive investigated were Poisson’s ratio and Young’s modulus GPa.

#### 3. FE Models

##### 3.1. Model 1 (Coarse Transition Mesh Design)

The FE mesh in model 1 was created using the PATRAN menu-driven FE pre- and postprocessing program operating in an X-window environment. Input into the program was the description of 14920 linear brick and triangular prism elements by indicating the material properties for the elements. The locating of nodal points was accomplished by dividing the configuration into 81 solid models. The original FE mesh of model 1 is shown in Figure 2.

Most of the geometry of the adherends and the adhesive was modeled using the 8-node solid elements. But at the transition zones from the adherends to the adhesive, where the mesh density is very high, some 6-node transition elements were used. Furthermore, the adhesive layer was divided into 40 equal parts along its length ( direction) and 20 equal parts along its width ( direction) in order to obtain an accurate indication of the variation of stresses in the lengthwise () and breadthwise () directions. Nodal points were located automatically by the PATRAN software as a function of the length and width of the adhesive layer, that is, in accordance with the geometric parameters of the model. Also the material parameters of the adhesive and adherends were input via the PATRAN software.

It can be seen from model 1 that the ratio of the thickness of the adherend elements to the thickness of the adhesive elements is 12.5. This is an abrupt transition in thickness. Also, the ratio of the lengths of the adhesive elements in the and directions to their thickness is 62.5 and 125, respectively. The adhesive elements are therefore very long and thin. This mesh is therefore regarded as a coarse mesh.

##### 3.2. Model 2 (Smooth Transition Mesh Design)

The FE mesh in model 2 was created using the ABAQUS FE analysis preprocessing program operating in an X-window environment. It was necessary to define the coordinates of the key nodes and the node number of the key elements in this case. Input into the program was the description of 2700 elements by indicating the material properties for the elements. The original FE mesh of model 2 is shown in Figure 3.

The geometry of the adherends and adhesive was modeled mainly using the 20-node solid elements. At the transition zones from the adherends to the adhesive, some 15-node transition elements were used. These transition elements were used only in the sections of the adherends which were outside the lap jointed section. The adhesive layer was divided into 10 equal parts along its length ( direction) and 10 equal parts along its width ( direction). Nodal points were located by the ABAQUS input file as a function of the length and width of the adhesive layer. The location of these points was in accordance with the geometric parameters of the model. Also the material parameters of the adhesive and adherends were input via the ABAQUS input file.

##### 3.3. Model 3 (Smoother Transition Mesh Design)

The FE mesh in model 3 was created also using the ABAQUS FE analysis preprocessing program operating in an X-window environment. Input into the program was the description of 16160 elements by indicating the material properties for the elements. The original FE mesh of model 3 is shown in Figure 4.

Most of the adherends and adhesive were also modeled using the 20-node solid elements. But at the transition zones from the adherends to the adhesive, some 15-node transition elements were used. In this case, the transition elements were used in both the lap jointed section of the adherends and the section of the adherends which were outside the lap jointed section. The adhesive layer was divided into 64 equal parts along its length ( direction) and 20 equal parts along its width ( direction) in order to obtain an accurate indication of the variation of stresses in the lengthwise and breadthwise directions. The nodal points were located by the ABAQUS input file as a function of the length and width of the adhesive layer, that is, in accordance with the geometric parameters of the model.

##### 3.4. Model 4 (Smoothest Transition Mesh Design)

The FE mesh in model 4 was created using the ABAQUS FE analysis preprocessing program operating in an X-window environment. Input into the program was the description of 57440 elements by indicating the material properties for the elements. The original FE mesh of model 4 is shown in Figure 5.

Again, the adherends and adhesive were mostly modeled using the 20-node solid elements. However, at the transition zones from the adherends to the adhesive, some 15-node transition elements were used. The adhesive layer was divided into 256 equal parts along its length ( direction) and 20 equal parts along its width ( direction) in order to obtain an accurate indication of the variation of stresses in the direction normal to the bond line. As in the previous cases, the nodal points and the material parameters of the adhesive and adherends were input via the ABAQUS input file.

##### 3.5. Comparison of FE Models

In order to illustrate the influence of the choice of FE model on the prediction of the mechanical behavior of adhesive joints, comparisons were performed between models with different modeling approaches. Since the failures of adhesive joints usually occur inside the adhesive layer, then only the lap jointed section is of interest.

It is easy to create a FE mesh using the PATRAN software as the nodal points are located automatically. The number of nodal points was reduced by this approach because part of the surface of the adherends and the surface of the adhesive uses the same node. In model 1, for example, input into the program was the location of 26922 nodal points but output of the analysis result was only 16968 nodal points. As a result, unfortunately, the number of nodal points was arranged discontinuously by the program. This makes postprocessing of the FE analysis results difficult. In fact it was found that it is better to use the ABAQUS FE preprocessing program to create the FE mesh of the single-lap adhesive joint.

Unlike model 1, model 2 was created using the ABAQUS FE analysis preprocessing program. This model has a limited number of elements and nodes. In the lap joint section, smooth transitions between the adherends and the adhesive in the direction were provided. Obviously the adhesive layer needs to be divided into more equal parts in order to obtain an accurate indication of the variation of stresses in the and directions.

It can be seen from Figure 5 that, in model 4, the FE model was refined gradually in steps from the adherends to the adhesive layer. That means smooth transitions were provided between the adherends and the adhesive in both the and directions. In the direction, the model was divided into more equal parts than in model 2. Therefore, model 4 was expected to provide more accurate analysis results. However, the disadvantage is that model 4 needs more computing time as it has a larger number of elements and nodes.

In the case of model 3, smooth transitions between the adherends and the adhesive were also provided in both and directions. In addition, model 3 has a moderate number of elements and nodes. Thus, model 3 was expected to be the most cost-effective of the 4 models studied.

In order to confirm this point, the stress distributions of the 4 models under tension were investigated. A distributed load of 1000 N was applied at the right end face of the upper adherend in the direction. This distributed load does not refer to any load condition in particular and is used simply as an example for comparisons between different modeling approaches and different combinations of elements. The boundary conditions of the joint are shown in Figure 6. Because five layers of solid elements were used across the adhesive thickness, a total of six interfaces were obtained. The lower interface, which is between the lower adherend and the adhesive, is denoted by interface 1. Similarly, the upper interface, which is between the adhesive and upper adherend, is denoted by interface 6. The intermediate interfaces are denoted by interfaces 2 to 5.

Figures 7 and 8 show the distributions of normal stress S_{11} predicted by the 4 models at interface 1 of the adhesive layer. The dimensions in the and directions are displayed in nondimensional form as and where is the width of joint and c is the length of the bonded section. Figure 7 shows the stress distributions at the front edge of interface 1. From the symmetry of the direction, it is clear that the stress distributions at the rear edge of interface 1 are the same as that at the front edge. Figure 8 shows the stress distributions at the centre line of interface 1.

From Figures 7 and 8, it can be seen that in the direction the left hand region is subjected to much higher stresses than the right hand region. In the case of model 1, the stress distribution curve does not extend as expected, possibly because the transition mesh is coarse. In the case of model 2, there is a severe oscillation in the stress distribution curves. The stress distribution curves of models 3 and 4, however, extend smoothly. Comparing the predicted direct stress for these models, it can be seen that the stress distributions of models 3 and 4 are similar. Therefore it can be inferred that the adhesive joint is better represented by models 3 and 4. The anomalies observed in stress distributions obtained using models 1 and 2 show that these two models are not adequate for modeling the adhesive joint. The figures show that the predicted results improve in accuracy as the model size increases. The worst result is obtained using model 1 which has a coarse FE mesh, whereas the best result is obtained using model 4 which has the finest FE mesh. Of course model 4 is more accurate in predicting the results than model 3, but it needs much longer CPU time. We conclude, therefore, that model 3 is the most cost-effective choice. Moreover, the results support a preference for choosing finer elements, for example, 20-node elements, to improve the accuracy of prediction rather than choosing more complex models which need much longer CPU time.

#### 4. Element Types and Topology

Stress/displacement elements were chosen for this study because they are suitable for modeling linear or complex nonlinear mechanical analyses that may involve contact, plasticity, and large deformations [15]. There are different types of stress/displacement elements, which are appropriate for different types of analysis. In the case of analysis of adhesive joints, the elements must accommodate the material properties and joint dimensions of both the adherend and the adhesive. The element combinations are defined in pairs to denote the element types used to model the adherend and adhesive, respectively. The first-order element combinations investigated in the present study include C3D8-C3D8, C3D8-C3D8H, C3D8R-C3D8H, C3D8R-C3D8R, and C3D8R-C3D8RH element combinations. The second-order element combinations investigated in the present study include C3D20-C3D20, C3D20-C3D20H, C3D20R-C3D20H, C3D20R-C3D20R, and C3D20R-C3D20RH element combinations.

#### 5. Discussion of Results

A recent study by the present author [16] showed that the spatial distributions of all the 6 components of stress are similar for different interfaces even though the stress values are slightly different. Since the failure of single-lap bonded joints initiates where high stresses occur, we are only interested in the maximum stresses. The overwhelming majority of maximum stresses occur at interface 1 though a few occur between interfaces 1 and 2 and at interface 6. Furthermore, the maximum stresses at interface 1 are much bigger. This section describes the predicted stress distributions obtained using some 3D stress/displacement element combinations to model a single-lap adhesive joint under tension.

##### 5.1. Stress Distributions Using First-Order Element Combinations

As stated previously, the first-order element combinations investigated in the present study include C3D8-C3D8, C3D8-C3D8H, C3D8R-C3D8H, C3D8R-C3D8R, and C3D8R-C3D8RH element combinations. The stress distributions corresponding to different element combinations were obtained. However, only a few typical distributions will be discussed here. Figure 9 shows the distributions of the 6 components of stress for the C3D8R-C3D8R element combinations at interface 1 as an illustration of the typical 3D stress distribution in the adhesive layer of the first-order element combinations. To enable easy comparison of these stress distributions, all the 6 components of stress are drawn using the same coordinate scales.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The figure shows that the highest stresses are concentrated near the left edge of the adhesive layer. The S_{33} has the highest magnitude of stress while the S_{12} has the least magnitude. The stress distributions of other first-order element combinations were omitted because they look similar to Figure 9, though there are some distinctions between them. These distinctions are discussed by using two-dimensional plots.

Figures 10 and 11 show two-dimensional plots of the maximum values for the 6 stress components of C3D8-C3D8, C3D8R-C3D8H, and C3D8R-C3D8R element combinations against the nondimensional distances and , respectively. In order to make the figures clear, the following codes are used to denote the element combinations and interfaces in Figures 10 and 11: 8-8: C3D8-C3D8 element combinations, 8R-8H: C3D8R-C3D8H element combinations, 8R-8R: C3D8R-C3D8R element combinations, Int1: interface 1, Int5: interface 5, Int6: interface 6.

The C3D8 element is an 8-node linear brick element. In the case of the C3D8-C3D8 element combinations, most of the maximum stresses (S_{11}, S_{33}, S_{12}, and S_{13}) occur at interface 1. The maximum stresses of S_{22} and S_{23}, however, occur at interface 6. In addition, all 6 maximum stresses occur at the left end of the adhesive layer which is closer to the clamped end.

Similarly, the C3D8R element is an 8-node linear brick, reduced integration with an hourglass control element, while the C3D8H element is an 8-node linear brick, hybrid with a constant pressure element. In the case of the C3D8R-C3D8H element combinations, , , , , and occur near the center of the left edge of interface 1, while occurs at interface 5.

In the case of the C3D8R-C3D8R element combinations, most of the maximum values of the 6 components of stress occur at the center of the left end of interface 1 except which occurs near the left-rear corner of the interface 1 and which occurs at the left-rear corner of the interface 6.

From Figure 10, it is obvious that the distributions of the direct or normal stresses S_{11}, S_{22}, and S_{33} are similar. The magnitudes of S_{11} and S_{22} are almost identical but the magnitude of S_{33} is more than double the magnitudes of S_{11} and S_{22}. While the distributions of the shear stresses S_{12}, S_{13}, and S_{23} are fairly similar, their magnitudes are widely different. Furthermore, the figures show that the stress distributions predicted by the 8R-8H and the 8R-8R element combinations are very closely correlated. But the stress distribution predicted by the 8-8 element combinations deviates significantly from predictions for the other element combinations. Similarly, it can be seen from Figure 11 that the stress distributions predicted by the 8-8 element combinations are quite different from those predicted by the 8R-8H and the 8R-8R element combinations, the predictions for which are fairly similar. We conclude, therefore, that the 8-8 element combination is not suitable for the analysis of single-lap adhesive joints.

##### 5.2. Stress Distributions Using Second-Order Element Combinations

The second-order element combinations investigated in the present study include C3D20-C3D20, C3D20-C3D20H, C3D20R-C3D20H, C3D20R-C3D20R, and C3D20R-C3D20RH element combinations. Figure 12 shows the distributions of the 6 components of stress for the C3D20R-C3D20R element combinations at interface 1 as an illustration of the typical 3D stress distribution in the adhesive layer of the second-order element combinations. It can be seen that the highest stresses are concentrated near the left edge of the adhesive layer. Also, S_{33} has the highest magnitude of stress whereas S_{12} has the least. These observations are identical with those made previously from Figure 9.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figures 13 and 14 show two-dimensional plots of the maximum values of the 6 stresses components of the C3D20-C3D20, C3D20R-C3D20H, and C3D20R-C3D20R element combinations against the nondimensional distances and , respectively. The stress distributions for the C3D20-C3D20H and C3D20R-C3D20HR element combinations are not included in these figures because they are very close to those of the C3D20-C3D20 and C3D20R-C3D20H element combinations, respectively. The following codes are used to denote the element combinations and interfaces in Figures 13 and 14: 20-20: C3D20-C3D20 element combinations, 20R-20H: C3D20R-C3D20H element combinations, 20R-20R: C3D20R-C3D20R element combinations, Int1: interface 1, Int6: interface 6.

The C3D20 element is a 20-node linear brick element. In the case of the C3D20-C3D20 element combinations, it can be seen that most maximum stresses occur at the interface 1 except the which occurs at the interface 6.

The C3D20R element is a 20-node quadratic brick, reduced integration element, while the C3D20H element is a 20-node quadratic brick hybrid with a linear pressure element. In the case of the C3D20R-C3D20H element combinations, the occurs near the center of the left region. The does not occur at an interface but between interfaces 1 and 2. In addition, there is a severe oscillation in the S_{13} curve.

In the case of the C3D20R-C3D20R element combinations, most of the maximum values of the 6 components of stresses occur in the center of the left end of interface 1. The shear stress (, ) occurs at the left-rear corner of interface 1 while (, ) occurs at the left-rear corner of interface 6. It is also clear from figures that the stress distribution curves extend smoothly.

##### 5.3. Comparison of Maximum Stresses Predicated by Linear and Quadratic Elements

In this section, comparisons are performed between the 8-node element groups and the 20-node element groups. Table 1 shows the maximum values of the 6 components of stress in the 8-node element combinations. It is clear that the stress state in this case is mainly dominated by the normal stress component S_{33} and then the shear stress component S_{13}. Surprisingly, the maximum values of the stress components S_{13} and S_{33} are higher than the maximum value of the stress component S_{11}. The latter would have been expected to be the most dominant since the joint is subjected to tensile loading. The departure from expected behavior is due to the effect of bending at the bonded section of the lap joint. In the 8-node element groups, not only the C3D8-C3D8 and C3D8-C3D8H element combinations, but also the C3D8R-C3D8R and C3D8R-C3D8RH element combinations have exactly the same stress distribution. This observation seems to suggest that the 8-node hybrid elements do not work in the analysis of the single-lap adhesive joints. In addition, the magnitudes of stresses of the 8-node element combinations oscillate in values along the lengthwise direction.

Table 2 shows the maximum values of the 6 components of stress of the 20-node element combinations. The stress components S_{13} and S_{33} mainly dominate the stress state as for the 8-node element groups, but the stress values are larger than those of the 8-node element groups. Unlike the 8-node element groups, however, every combination has a different stress distribution. For example, Table 1 shows that the maximum values of the 6 stress components predicated by the C3D8-C3D8 and C3D8R-C3D8R element combinations are identical to those predicated by the C3D8-C3D8H and C3D8R-C3D8RH element combinations, respectively. However, Table 2 shows that the maximum values of the 6 stress components predicated by the second-order element combinations are all different. Also, for the C3D20-C3D20H, C3D20R-C3D20H, and C3D20R-C3D20RH element combinations, the maximum values of normal stress component S_{22} do not occur at an interface but between interfaces 1 and 2. This observation seems to suggest that the 20-node hybrid elements do not work well in this study. This is not very surprising because the adhesive used is very stiff and therefore has a Poisson ratio less than 0.5. Thus the adhesive is not incompressible and the hybrid elements formulated for incompressible or nearly incompressible materials will not work well.

Also, it was shown previously that the C3D8-C3D8 element combination is not suitable. Thus, the reasonable choice should be between the following three types of element combinations: C3D8R-C3D8R, C3D20-C3D20, and C3D20R-C3D20R.

As mentioned before, second-order elements provide higher accuracy than first-order elements. They capture stress concentrations more effectively and are better for modeling geometric features. In addition, second-order elements are very effective in bending-dominated problems.

Reduced integration uses a lower order integration to form the element stiffness and so reduces the run time, especially in three dimensions. For example, element type C3D20R has 8 integration points while C3D20 has 27. Therefore, element assembly is roughly 3.5 times less costly for C3D20R than for C3D20. In addition, second-order reduced integration elements generally yield more accurate results than the corresponding fully integrated elements.

Finally, C3D20R-C3D20R element combinations would be the best element combinations for the analysis of single-lap adhesive joints.

#### 6. Conclusions

Some FE models for analyzing the behavior of adhesive joints were described in this paper. In these models five layers of solid elements were used across the adhesive layer which was only 0.05 mm thick, in order to obtain accurate results. The FE models were refined gradually in steps from adherends to adhesive layer. Most of the adherends and the adhesive were modeled using quadratic solid elements but some triangular solid elements were used to give a smooth transition. Comparisons were performed between models with different modeling approaches as well as different types of element combinations.

From the comparisons between the FE models, it is clear that of the 4 models presented in this study model 3 is the most cost-effective. This is because it has a moderate number of elements and nodes and a smooth transition between the adherends and the adhesive in both and directions.

The results of the analysis also show that the linear, fully integrated, and hybrid elements are not suitable for the analysis of single-lap adhesive joints. In addition, second-order reduced integration elements provide higher efficiency and accuracy than the corresponding first-order elements.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This study is partially supported by National Natural Science Foundation of China (Grant no. 50965009) and Special Program of the Ministry of Science and Technology, China (Grant no. 2012ZX04012-031).