This article mainly analyzes the free vibration characteristic of the sandwich piezoelectric beam under elastic boundary conditions and thermal environment. According to the first-order shear deformation theory and Hamilton’s principle, the thermo-electro-elastic coupling equations of the sandwich piezoelectric beam are obtained. Meanwhile, elastic boundary conditions composed of an array of springs are introduced, and the displacement variables and external potential energy of the beam are expressed as wave functions. By using the method of reverberation-ray matrix to integrate and solve the governing equations, a search algorithm based on golden-section search is introduced to calculate the required frequency parameters. A series of numerical results are compared with those reported in literature studies and obtained by simulation software to verify the correctness and versatility of the search algorithm. In addition, three parametric research cases are proposed to investigate the frequency parameters of sandwich piezoelectric beams with elastic restraint conditions, material parameters, thickness ratio, different temperature rises, and external electric potential.

1. Introduction

Piezoelectric ceramics are an electronic ceramic material with positive and negative piezoelectric effects; they have been widely used in sensors, energy harvesting devices, and smart actuators. Given their piezoelectricity, dielectric nature, elasticity, and other characteristics, piezoelectric ceramics are highly transformable in terms of size, external parameters, and lamination methods. Piezoelectric materials and general materials are usually laminated or pasted in engineering to achieve improvement in their performance. The ultrasonic motor widely used in the aerospace and automobile field is composed of two layers of piezoelectric sheets and a base metal material in the middle layer. The structure under such high-temperature and pressure environment is prone to overloading and damage. Therefore, a comprehensive understanding for vibration characteristics of the piezoelectric laminated structure in a thermal environment is necessary to reduce resonance through the rational design of structural parameters.

Given that the aspect ratio of the piezoelectric laminated structure is generally greater than 3.5, it can be suitably simplified into a beam model. The related beam theories mainly include elasticity theory [1, 2], Euler–Bernoulli theory which is also named classical beam theory (CBT) [3, 4], first-order shear deformation theory (FSDT) [58], and high-order shear theory [9, 10]. The early CBT was suitable for studying the characteristics of slender beams. Given that the theory disregards the influence of transverse shear deformation, it is unsuitable for studying medium-thickness beams. As a remedial measure, the Timoshenko beam theory was obtained by adding shear deformation to the Euler–Bernoulli equation. Thus, the rotation of the cross section is caused by deflection and transverse shear. Different cross sections are needed to introduce the shear correction factor to reduce errors. For the higher-order shear theory without using the shear correction factor, different shape functions are introduced for the shear strain, resulting in the more accurate calculation of laminated plate problems.

Given the various beam theories mentioned above, numerous researchers have achieved remarkable results in the analysis of piezoelectric beams. From the perspective of their shape, beams can be divided into bimorph beams [1113], functionally graded piezoelectric beams (FGPM) and sandwich beams [1422], nano-microbeams [2327], and curved beams [28, 29]. The above literature shows that the transition from the traditional composite material lamination theory to the piezoelectric material lamination structure has been relatively complete. Carrera et al. [30] explained in detail the derivation process of piezoelectric and governing equations when smart materials are applied to beam structures and perfected the constitutive equation of the piezoelectric laminated structure under the above deformation theories. On this basis, the analysis of the dynamic characteristics of piezoelectric beams in various coupling environments has become an important research topic in recent years. Annigeri et al. [31] used the in-plate finite element method (FEM) to analyze the vibration characteristics of magneto-electroelastic beams based on CBT. He et al. [32] studied the free vibration of thermo-magneto-electro-nanobeams based on the wave solution method and carried out a series of parametric research including slenderness ratio, thickness, temperature rise, and external potential under different boundary conditions. Yang and Xiang [33] used the differential quadrature method (DQM) to study the FGPM beam under thermoelectric coupling. Lee [34] studied the displacement and stress responses of the FGPM beam with layer-wise finite element formulation and analyzed the effects of different loads on the deflection and stress of bimorph actuators. Benjeddou et al. [35] used the strain gradient theory to study the dynamic responses of the cantilever beam and discussed the effect of structure size and electric potential. Summarizing the results of these researchers, the external conditions for piezoelectric beams generally include the thermal environment, electric field, and magnetic field. However, these studies mostly focused on traditional functionally graded structures or microstructures, and a relatively limited number of research studies centered on sandwich structures under the above conditions. Therefore, the dynamic characteristics of sandwich piezoelectric beams under a coupling environment must be analyzed, and a parameterized study on their geometric properties and external conditions must be conducted.

The analysis of the free vibration of this type of structure can ultimately be attributed to the problem of solving differential equations. In addition to various beam theories mentioned above, various analysis methods are used to solve the dynamic characteristics of composite material laminated beams under coupling fields; such methods include the FEM [3640] and DQM [4143], modal superposition method [44, 45], transfer matrix method [4648], Galerkin method [49], Ritz method [50], and pseudo-arclength continuation method [51]. However, the modal superposition method has the problem of a large amount of calculation and low accuracy, the transfer matrix method also has the problem of insufficient high-order frequency accuracy, and the differential quadrature method is not convenient when calculating general boundary conditions. The research on structural boundary conditions also stays in traditional clamped, simply supported, and free. In actual engineering, there are more elastic connection ways such as bolt, rivet, and glue connections instead of the above methods. Therefore, an artificial boundary realized by changing the spring stiffness coefficient is introduced to obtain more accurate calculations. In addition, an analysis technique based on the frequency domain named the method of reverberation-ray matrix (MRRM) was proposed by Howard and Pao [52], which was used to solve the problem of elastic wave propagation in frame structures. The basic idea is to establish a pair of coordinate systems on the structure, derive the phase and scattering relationships of the structure through the dual relationship under two coordinate systems, and obtain the overall MRRM equation containing a series of dynamic properties of the structure. Guo et al. [53] discussed the MRRM in the functionally graded layered medium formula and equation-solving method. Miao et al. [54] used the MRRM to establish the frequency response function of the displacement in the frequency domain. The frequency is determined by the frequency response curve, and the mode shape can be calculated by the joint matrix from the system control equation. Thus, the MRRM is highly suitable for solving the dynamic characteristics of related structures. The scattering matrix is an all-constant matrix resulting from the dual coordinate system, its numerical value is relatively stable, and the high-order frequency can be accurately solved. However, limited research studied the dynamic characteristics of sandwich piezoelectric structures with arbitrary boundary conditions based on the MRRM.

As in previous summaries, this article uses the MRRM to analyze the free vibration characteristics of the sandwich piezoelectric beam with elastic boundary conditions and coupling field. In accordance with the FSDT and Hamilton’s principle, the thermo-electric-elastic coupling equation and the boundary coordination relationship are obtained. Three sets of mechanical springs are used to implement any combination of the boundary conditions. On these bases, the frequency results are compared with the reported literature studies and software. For enriched practicability and flexibility of this theory, a series of parameters that affect the vibration characteristics, including the elastic modulus, thickness ratio, ply orientation of the middle metal layer, elastic boundary conditions, temperature rise, and external electric potential, is investigated.

2. Mathematical Model

2.1. Model Description

Figure 1 shows a sandwich piezoelectric beam with piezoelectric sheets as the top and bottom layers and a general metal material layer in the middle. The length and thickness are set as L and h, respectively. The piezoelectric and metal layers have a thickness of hp and hm, respectively. Each layer is assumed to be perfectly pasted, and the influence of the adhesive layer is disregarded. For simulation of the boundary force at both ends of the beam, there are three sets of mechanical springs, including two sets of displacement springs and a set of rotary springs, which are set in each group. Furthermore, a set of electric displacement constraints is considered. In the following derivation, the piezoelectric material equation will be described first. For general metal layers, only the piezoelectric term in the equation needs to be removed. Moreover, zk+1 and zk indicate the distance from the top and bottom of the kth layer to the middle face of the beam, respectively.

2.2. Kinematic Equation

The external electric potential and temperature rise are introduced to describe the thermo-electro-elastic coupling equation of the sandwich piezoelectric beam. According to the assumption of the FSDT, the displacement of any point can be written as follows:

In (1), are the displacements and rotation angle of a node in the middle plane. The variables of external electric potential are given by the following [32]:where and Φ0 is the initial external electric potential.

According to (1) and (2), the mechanical strains and electric field are, respectively, obtained as follows:

Based on the stress-strain relationship of orthotropic materials and thermal stress, the constitutive equations of the kth layer from the sandwich piezoelectric beam can be written aswhere the overline indicates the transformed coefficient. is the elastic constants; and are piezoelectric and dielectric constants, respectively. Notably, and are equal to zero when the kth layer is a general metal material; these constants can be specifically expressed as follows [23]:where and represent the normal and shear stresses, respectively. and represent the electric displacement. and are neglected, and . and correspond to the thermal expansion and pyroelectric constants, respectively. Thus, the thermal strain can be written as

By integrating the mechanical stress in (4) in the thickness direction, the resultant expression of stretching stress, moment, and transverse force can be calculated with the following:

Similarly, the corresponding stiffness coefficients can be obtained in the same integration way:where represents the shear correction factor that is determined by a cross-sectional shape. , , and refer to the mechanical stiffnesses, and and denote the electromechanical stiffnesses. stands for the dielectric stiffness [30]. Therefore, we can obtain the constitutive equations of the piezoelectric laminated beam as follows:

2.3. Governing Equations and Geometric Boundary Conditions

The coupling equation under thermal and electric fields of the sandwich piezoelectric beam is obtained by Hamilton’s principle defined by the following:where represents the variation and indicates the virtual strain energy and virtual electric internal energy calculated by the following [30]:

Substituting equations (3), (4), and (7) into (11) and performing variational operations on the above equation, the expression of total virtual internal work is shown as

and represent the external thermal potential energy and electric potential energy, respectively. In this article, only the axial strain component should be considered, and it can be expressed as

Thus, the external thermal and electric potential energies can be written with the following equations:where and are the normal forces as

From the beam displacement field in (1), represents the kinetic energy of the sandwich piezoelectric beam:where the material density of the kth layer is expressed as , and the superscript indicates the integral over time.where are the kth layer’s inertia coefficients defined by

Equations (12), (14), and (17) are substituted into equation (10) to derive the governing equations of the sandwich piezoelectric beam under the coupling environment as follows:

According to the first term of equation (19) and the basic principle of variations, the coefficients of , , , and should be set to zero and can be written as

The coordinate boundary conditions of the sandwich piezoelectric beam are determined from the second term of (19). The ends of the beam combined with the corresponding spring constraint can be expressed as :

Therefore, setting different spring stiffness values is equivalent to applying different forces at both ends. Clamped (C) and free (F) correspond to a large value of k and 0, respectively.

2.4. Wave Solutions

The Laplace transformation of (19) is carried out to promote the MRRM in Section 4, and the governing equations in the frequency domain are shown as follows:where the expressions of can be calculated with the following:

In general, the displacement variables that satisfy the governing equations can be set aswhere is the characteristic wave number and , , , and are the wave contribution factors. Submitting equation (24) into equation (22), a set of homogeneous functions about wave amplitude is obtained as follows:where the expressions of can be written as follows:

For the general solution of (25), we can obtain the characteristic wave number by setting as zero:where are the coefficients corresponding to the wave numbers determined by matrix T. The four pairs of wave number and corresponding eight eigenvectors can be obtained by equation (27). According to the general method of solving matrix equations, the elements in can be solved by making one element equal to 1. Thus, the basic solution for the displacement variables can be defined as follows:where the coefficients are given by the following:

Based on equation (28) and the superposition principle, the displacement variables can be written as follows:where are the amplitude of the arriving waves and are the departing waves.

3. Solving Method

Before using the MRRM to formulate the structure, the displacement vector δ and force vector f are introduced first, and equation (30) is rewritten in the following form:whererepresent the arriving wave vectors, departing wave vectors, and phase propagation matrix. and are the generalized displacement coefficients and generalized force coefficient matrix, respectively:where

3.1. Generalized Phase Relationship

As the first step in constructing the MRRM formulation, a pair of local dual coordinate systems should be established for each unit on the structure (Figure 2). Nodes 1 and 2 are considered the two ends of the structure, and the positive direction is set from node 1 to node 2 as x12. The coordinate systems at nodes 1 and 2 are defined as and , respectively. If the first coordinate system is defined to be parallel with the positive direction, then the stiffness coefficient and inertia moment of the structure in the “21” coordinate system will change accordingly. The related coefficients in the “21” coordinate system can be written as follows:

Therefore, the relations of the generalized displacement coefficients and in the dual coordinate system can be obtained as follows:where matrix and . Although two wave solutions can be obtained, the physical quantity at the same point on the beam is constant, and the vectors and in the dual coordinate system are expressed with the following:

Before deriving the phase relationship of the beam structure, the wave solutions need to be classified first. When the wave direction is the same or opposite with , it is called the departing or arriving wave, respectively. Thus, substituting equation (37) into equation (31), the relationship between and can be simplified and classified as follows:where is the phase matrix of the sandwich piezoelectric beam.

3.2. Scattering Relation

Using only the phase relationship of the structure is insufficient to obtain the dynamic characteristics of the sandwich piezoelectric beams. The scattering relationship at each node can be obtained through geometric equations and force balance conditions. Given that the spring boundary exists at nodes 1 and 2 on both ends, equation (31) can be used to establish the relationship between generalized displacements and forces through the spring stiffness matrix:

Thus, we can obtain the following:

Table 1 provides the stiffness matrix including the C, S, F, and elastic support (E) for common boundary conditions in engineering. The value of 1018 represents a very large stiffness coefficient to achieve the C condition because an infinite value cannot be calculated. The scattering relation can be obtained by classifying vectors aij and dij in the same manner:where S represents the scattering matrix of the sandwich piezoelectric beam.

3.3. Reverberation-Ray Matrix and Free Vibration

According to the relation of equations (39) and (42), the overall MRRM equation of the piezoelectric laminated beam without considering external excitation can be obtained as follows:where and represents the identity matrix. is called the reverberation-ray matrix of the present sandwich piezoelectric beam. For free vibration analysis, the determinant of equation (43) is set to zero:where the solution of equation (44) is the natural vibration parameter of the sandwich piezoelectric beam.

Equation (44) is a complex transcendental equation about frequency ω and is not easily solved by general analytical methods. However, for this complex equation, as ω gradually increases, the modulus of the complex determinant always assumes a minimum value, which is zero when it passes through its zero point. Therefore, the golden-section search (GSS) method can be used to search for natural frequencies. The evaluation function is defined as follows:

The specific process is performed as follows:Step 1: define the initial trial frequency , frequency search step size , and maximum frequency  Step 2: generating two search frequencies and , calculate the evaluation function values corresponding to the three frequencies  Step 3: if and , proceed to Step 4; else, return to Step 2 with resetting Step 4: execute the GSS method to search the undetermined frequency , and calculate the value of function  Step 5: when the maximum frequency is reached, stop the search; else, return to Step 2 to reset

Figure 3 shows the specific flowchart of the GSS algorithm.

4. Result Verification and Discussion

4.1. Result Verification

Currently, several types of numerical examples are used to study the free vibration characteristics of sandwich piezoelectric beams. First, the calculated results are compared with the findings of literature studies and software simulations to verify the calculation correctness. On this basis, a series of parameter coupling studies, including elastic modulus, ply orientation of materials, thickness ratio, external electric potential, and temperature rise, is considered. The effect of single and coupling factors on vibration characteristics will be investigated separately. Table 2 lists the relevant material parameters for verification and subsequent parameterization studies. The nondimensional results will be used to express the frequency parameters if no extra instructions are given in the following examples.

First, a double-layer composite beam composed of materials A and B is introduced without the thermal environment and piezoelectric effect. As mentioned earlier, the general material layer only needs to set the piezoelectric coefficients in the equation to zero. He et al. [32] analyzed this example based on the DQM and considered the boundary conditions and slenderness ratio. Table 3 shows the comparison results of the first six dimensionless frequencies under different boundary conditions and slenderness ratios, where ω is defined as . The numerical errors are in a relatively small range, and the poorest result error is at 1.7%.

Next, the material in the above numerical example is replaced with piezoelectric material, a double-layer beam with equal thickness composed of PZT-5A analyzed by Fernandes and Pouget [11]; in consideration of the Bernoulli–Euler beam theory, a purely mechanical beam theory and an electromechanical coupling model under an electric field environment are established, and the related material parameters can be found in Table 2 (dimensions: L = 0.4 m and h = 0.01 m). Table 4 lists the first six dimensionless frequencies in different mechanical boundary conditions, which show that our frequency parameters are notably close to those of Michael Krommer and FEM simulation results. In the table, “NE” stands for the layers not equipped with electrodes, which are consistent with the model in this example.

Finally, a sandwich piezoelectric laminated beam with L = 0.4 m and h = 0.03 m with equal thickness of each layer composed of PZT-5A/aluminum/PZT-5A is considered. This verification example and the following parameter studies will all use this model unless special instructions are specified. In the comparison with Ansys v18.2 software, 80 solid 5-grid elements and 40 solid 45-grid elements are used to simulate the piezoelectric and metal layers, respectively. The external coefficients Φ0 and △T are set to 10 kV and 25–200°C. Table 5 shows the first six dimensionless natural frequencies compared with Ansys v18.2 software under the clamped-clamped boundary condition. The frequency values of the two variables are extremely close, and the largest error is 0.47%.

Based on the above comparisons and verifications, the proposed method is a simple and effective technique for the analysis of the vibration characteristics of piezoelectric composite beams, and it provides a convenient adjustment interface for various boundary conditions and coupling of different external conditions or parameters. In this manner, the following parameterization studies require minor changes to improve the computational efficiency and accuracy.

4.2. Case 1: Elastic Modulus and Boundary Conditions in the Thermal Environment

Different from general composite materials or functionally graded beams, the sandwich piezoelectric laminated structure has a good design capability because it is composed of different kinds of materials, and various parameters have varied effects on its dynamic performance. This paper mainly studies five aspects of variables, including temperature rise, external electric potential, material and geometric properties of the middle layer, elastic boundary conditions, and the free vibration characteristic under coupled fields.

For the detailed investigation of the effect of spring restraint constant on frequency characteristics, three elastic boundary conditions (E1, E2, and E3) are additionally defined, and they correspond to the z-direction spring constant, rotational spring constant, and condition when both constants are considered. Table 1 provides the detailed essential conditions and spring coefficients. With the calculated natural frequency under different boundary combinations, the changing rules in Figure 4 show that the frequency parameters of each order increase within a certain elastic coefficient interval that is generally between 105 and 1010. The interval gradually moves to the right as the frequency order increases. Furthermore, different boundary combinations will result in different frequency characteristics. Figures 4(a) and 4(b) show the independent influence of and kθ when one end is C or S, respectively. Figure 4(d) shows the effect of and kθ acting at the same time when one end is C or S. Figures 4(c) and 4(e) show that the effect of both ends is elastic restraint.

Next, the effects of elastic boundary and traditional boundary conditions on frequency parameters are studied by giving a specific spring constant. Based on the spring constants corresponding to the frequency variation interval in Figure 4, two defined elastic boundary conditions Ea and Eb are considered besides the conventional S, C, and F conditions, where Ea indicates that only the z-direction is elastically restrained (), and the spring constant kw is set as 106; Eb indicates that the rotation direction and z-direction are elastically restrained (, θ≠0), and kw and kθ are set as 106. Figures 5(a)5(d) show the frequency parameters under different boundary conditions. Figures 5(c) and 5(d) illustrate that the frequency parameters can be further adjusted by changing the spring constant.

On the basis of the above studies, we further consider the effect of the properties of the middle metal layer under the thermal environment on the frequency parameters. First, the elastic modulus of the middle layer is set from 30 GPa to 350 GPa. However, the density and Poisson’s ratio remain unchanged, and the range of temperature rise is 0°C–200°C. Figure 5 also shows the change trend of frequency parameters under different elastic modulus and temperature rise conditions. The thermal environment has a negative effect on the free vibration of the piezoelectric beam because the stiffness is weakened when the temperature rise causes internal stress in the beam. From the horizontal view in Figure 5, when Poisson’s ratio and density of the core layer are constant, the free vibration parameters will increase with the increase of its elastic modulus, and high-order frequencies are more evident than the low-order ones.

4.3. Case 2: Ply Orientation and the Number of Layers under Different Electric Potentials

In this section, the middle metal layer is regarded as a composite layer composed of anisotropic materials (Figure 6(a)). The effects of the lamination method, ply orientation, the number of layers, and external electric potential on the frequency parameters are discussed. Table 2 gives the parameters of material-D, and there are five types of lamination method, denoted as θ, 0/θ, θ/0/θ, 0/θ/0/θ, and θ/0/θ/0/θ, where θ represents the ply orientation to the x-direction which is set as 0°–180°. The external electric potential is set from −20 kV to 20 kV. Figure 7 shows the frequency parameters of the clamped-clamped (C-C) beam under the five types of lamination methods. Figure 7(a) shows that the frequency parameters are symmetrical at about x = 90°. For different lamination methods, the frequency range of beams with more layers is significantly larger than those with few layers. However, the overall trend is similar. From the horizontal view, a ply orientation with the x-axis has a negative effect on free vibration parameters, and the slope is large between 30°–60° and 120°–150°. Figure 7(b) shows the effect of ply orientation under different lamination methods. When the size of the laminated beam is constant and the angle is zero, the effect of the number of middle layers on free vibration can be obtained. It is obvious that more middle layers indicate high-frequency parameters. When the ply orientation becomes larger, the change mode of different lamination methods varies although the overall trend of the frequency is still increasing. The comparison reveals that it is related to the proportion of the zero-degree layers to the total layers, and the range of frequency changes decreases as this ratio increases. In Figure 7(c), external electric potential is similar with temperature rise that has a negative effect on free vibration parameters. However, the overall effect is notably small.

4.4. Case 3: Thickness Ratio and Composition Method in the Thermal Environment

Next, we study the frequency characteristics of the thickness ratio under different material combinations in two cases. Similar to Case 1, the middle layer is composed of isotropic materials, and the beam size is constant. As shown in Figure 3(b), six thickness ratios (hm/hp = 0.5, 1, 2, 5, 20, and 50) of the metal layer to the piezoelectric layer are considered by calculating the frequency parameters under different thickness ratios and the number of layers under C-C boundary conditions. Figure 8 illustrates that a piezoelectric laminated beam with a large thickness ratio has a high frequency at any number of layers. However, the growth rate decreases as the thickness ratio increases and becomes more evident with the increase in middle layers. Furthermore, if the middle layer is composed of isotropic and anisotropic layers with hm and hc for each layer as shown in Figure 3(c), then the effect of different numbers of layers and thickness ratios (hm/hc) on the frequency parameters should be reconsidered, and the temperature rise is set from 0°C to 200°C. Figure 9(a) shows that when piezoelectric, anisotropic material, and isotropic material layers together form a symmetrical composite structure, the frequency characteristics differ from those of previous examples. For each order frequency, the thickness ratio has a positive effect on free vibration parameters, and the value is between 1 and 5. If the thickness ratio is less than 1 and greater than 5, the frequency will change in the opposite direction. In addition, the change trends in Figures 9(c) and 9(d) are similar with the rules in Case 1. The temperature rise still has a negative but negligible effect on the frequency parameters in this example.

5. Conclusion

This article mainly studied the free vibration of thermo-electro-elastic piezoelectric beams with the sandwich structure. To obtain more general characteristics and simplify the study of related parameters, we introduced the spring boundary to simulate the boundary force at both ends and integrated the phase and scattering relationships of the structure to achieve the total MRRM equation. By performing an algorithm based on the GSS for the MRRM equation, the natural frequency of the structure was calculated and compared with the FEM results and existing literature to verify the correctness of the proposed method. To further explore the effects of geometric factors and material parameters on the vibration of the structure, we carried out a series of analyzes of the geometric and material characteristics of the structure, especially the middle metal layer. These parameters included elastic modulus, the number of layers, thickness ratio, ply orientation of anisotropic material, different boundary conditions, each spring constant, temperature rise, and external electric potential of the beam. Through the coupling of external environmental parameters and geometric parameters, we arrived at the following conclusions:(1)The free vibration parameters change obviously when the spring coefficient range is generally between 105 and 1010.(2)For sandwich piezoelectric beams, the elastic modulus, thickness, and the number of layers of the middle metal layer have a positive effect on the free vibration.(3)When the middle metal layer is composed of an anisotropic material, the ply orientation and lamination method influence the free vibration significantly. The more layers that change the ply orientation, the larger the range of frequency fluctuation.(4)If the middle layer is composed of isotropic and anisotropic materials, then its thickness ratios have different effects on the frequency parameters in different numerical ranges, and additional layers will increase the frequency parameters.(5)Factors in the external coupled environment, such as temperature and electric potential, have a negative effect on free vibration parameters of sandwich piezoelectric laminated beams. The temperature rise and external potential increase will cause a negligible decrease in frequency.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (Grant no. 51905511).