#### Abstract

Based on the principle of energy variation, an improved Fourier series is introduced as an allowable displacement function. This paper constructs a calculation model that can study the in-plane and out-of-plane free and forced vibrations of curved beam structures under different boundary conditions. Firstly, based on the generalized shell theory, considering the shear and inertial effects of curved beam structures, as well as the coupling effects of displacement components, the kinetic energy and strain potential energy of the curved beam are obtained. Subsequently, an artificial spring system is introduced to satisfy the constraint condition of the displacement at the boundary of the curved beam, obtain its elastic potential energy, and add it to the system energy functional. Any concentrated mass point or concentrated external load can also be added to the energy function of the entire system with a corresponding energy term. In various situations including classical boundary conditions, the accuracy and efficiency of the method in this paper are proved by comparing with the calculation results of FEM. Besides, by accurately calculating the vibration characteristics of common engineering structures like slow curvature (whirl line), the wide application prospects of this method are shown.

#### 1. Introduction

As a common structure in the engineering field, curved beam structure has attracted widespread attention for its dynamic characteristics. Unlike the straight beam structure, due to the influence of deformation, there is considerable coupling between the vibration displacement components of the curved beam structure. To simplify the calculation, the axial tensile deformation, shear deformation, and rotational inertia are sometimes ignored [1–5], but this may result in a large error in calculating the high-order natural frequency of the curved beam [6]. Besides, since the vibration of the plane-curved beam is not coupled between the in-plane and out-of-plane, most scholars discuss it as two relatively independent issues [7, 8], and few documents have studied both at the same time [9].

For the problem of in-plane vibration of curved beams, if all influencing factors such as axial tension are taken into account, discrete Green function combined with numerical integration was used in early research by Kawakami et al. [9]. Chen and Shen [10] used the trigonometric function combined with Galerkin’s method, while Austin and Veletsos [11] used the energy approximation method to deal with this kind of problem. At present, the most commonly used methods by researchers in this field are the dynamic stiffness method [12], transfer matrix method [13], and FE method. As one of the most widely used numerical methods, FEM has been used for the analysis of in-plane vibration of curved beams almost since its birth [14]. To further improve the calculation accuracy and efficiency of this method, many new curved beam elements have been proposed in turn. Based on the Timoshenko beam theory, Raveendranath et al. [15] developed a three-node curved beam element. Based on the extended Hamilton principle, Yang et al. [2] developed a high-order Lagrangian-type element suitable for variable curvature curved beams. Cannarozzi and Molari [16] proposed a stress-mixed element based on the improved Hellinger–Reissner function and proved that the element can be used to calculate arbitrary geometrical curved beam structures through calculation examples. Saffari et al. [17] regard the curved beam structure as a combination of arch elements, each basic arch element introduces shear and axial deformation energy, and the elements are connected by a transformation matrix constructed by curvature. Kim et al. [18] introduced additional node-independent degrees of freedom in the displacement interpolation function, which significantly improved the numerical accuracy of the displacement-stress hybrid element, making it more suitable for predicting higher-order modes. Based on the B-spline interval wavelet (BSWI), Yang et al. [19] transform the displacement field of the element from wavelet space to physical space, instead of the traditional polynomial interpolation finite element. Due to the good approximation performance of BSWI, the new element has better accuracy and efficiency.

Similar to the problem of in-plane vibration, FEM is most widely used to solve the problem of out-of-plane vibration of curved beams. However, it is worth noting that because the coupling term in the out-of-plane vibration is bending-torsion coupling, conventional curved beam finite element elements may have shear, bending, and torsion locking phenomena when dealing with this problem. The locking phenomenon will seriously affect the accuracy and convergence [20, 21]. Ishaquddin et al. [8] proposed a finite element model based on independent field interpolation, which can accurately predict the dynamic characteristics of curved beam structures without generating any false modes under the limit stiffness ratio.

Although the numerical method represented by FEM has made considerable achievements in solving the problem of curved beam vibration, the (semi)analytic method still has its unique advantages in mechanism interpretation and parameter research. As a widely used method, the energy method has many research results in the field of structural vibration calculation [22, 23]. However, the traditional energy method needs to select displacement functions according to different boundaries when dealing with structures under different boundaries, and the applicable range of its calculation model is relatively limited. To solve this problem, Qu [24] proposed a modified variational method, which is based on the modified variational principle and the weighted residual method to solve the internal interface and boundary constraints of the structure, thereby relaxing the requirements for the selection of structural allowable functions. Subsequently, Su et al. [25] used this method to solve the in-plane vibration problem of curved beams and used the advantages of flexible application range to study the influence of concentrated mass points on the vibration characteristics of curved beams. Also, Li [26] proposed an improved Fourier series, which is a combination of traditional Fourier series and auxiliary functions and used it to analyse the lateral vibration characteristics of a straight beam under arbitrary supports. Subsequently, using the series as the allowable displacement function, combining artificial spring system, and basing the Hamilton principle, the vibration characteristics of the elastically supported FGM beam structure [27], the Timoshenko beam structure [28] with arbitrary cross-sectional shape, and the Euler curved beam structure [5] in the plane were also studied.

In this work, all the vibration displacement allowable functions of curved beam structures will be uniformly expressed in the form of improved Fourier series. Different from solving the partial differential equation directly, this paper will be based on the energy variation principle and study the free or forced vibration characteristics of the curved beam system. Firstly, according to the generalized shell theory, the displacement-strain relationship considering the shear and inertial effects is obtained, and the kinetic energy and strain potential energy of the curved beam structure are obtained. Subsequently, the boundary constraint condition or the continuous condition of the curved beam connection will be added as the elastic potential energy via the artificial spring system to the energy functional. Finally, the concentrated mass point or concentrated dynamic load will be expressed in terms of additional kinetic energy or external force work.

The comparison with [5] and the FEM calculation results proves the correctness of the calculation model in this paper, and the improved Fourier series also shows good convergence. Besides, the artificial spring system at the boundary can accurately simulate various classical or elastic boundary conditions. By setting different constraint stiffness values, the boundary conditions of the model can be directly modified without rederived. Compared with the common (semi)analytic method that needs to reselect the displacement function, this paper successfully constructed a unified calculation model with a wide range of applications. Due to the small number of series cutoffs required, the mass and stiffness matrix formed by the model in this paper will be significantly smaller than FEM. Therefore, the model in this paper has more efficiency advantages in the calculation of forced vibration. Examples that consider changes in curvature or section parameters show the broad prospects of the model in practical engineering applications.

#### 2. Theoretical Method

##### 2.1. Model Description

As shown in Figure 1, the research object of this paper is a curved beam with initial deformation in the *xy* plane under the right-handed coordinate system *O-xyz*. Generally, the vibration in the *xy* plane, that is, the vibration in the plane where the deformation is located, is called in-plane vibration, and the plane vibration perpendicular to it is out-of-plane vibration. There are 6 vibration displacement components at any point on the centerline of a curved beam. Although these displacement components may have a coupling relationship with each other, the in-plane and out-of-plane displacements can be decoupled. Specifically, there are three in-plane displacements, which are tangential displacement *u(s)*, radial displacement , and bending rotation displacement , and out-of-plane displacements, which are out-of-plane transverse displacement , bending rotation displacement , and torsional rotation of the cross-section . *R(s)* is the radius of curvature at any point of the curved beam. *L*_{s} is the length of the curved beam. *s* is the beam axis curvilinear path.

##### 2.2. Variational Formulation Model of the Curved Beam

###### 2.2.1. In-Plane Problem

In this paper, an artificial spring system will be introduced at both ends of the curved beam structure to simulate the actual boundary. In the problem of in-plane vibration, both the left and right end faces of the beam will use one rotation spring and two linear displacement springs to restrain the rotation displacement and the tangential (normal) displacement, respectively. The curved beam structure supported by springs is shown in Figure 2.

In the above figure, *k*_{u0}, , and are the constraining spring groups at the left end of the curved beam (at the starting point), which, respectively, constrain the displacement and rotation displacement in *u* and directions. The definitions of *k*_{u1}, , and are similar, acting on the right end (endpoint) of the curved beam. The unit of displacement spring is N/m, and the unit of rotation spring is N/rad. *f*_{u} and are external concentrated dynamic loads, and the relevant research in this paper will adopt unit amplitude excitation force. *m*_{p} is the quality of the concentrated quality point.

Considering the effects of shear and inertia, this paper will construct the energy functional in representing the system’s in-plane vibration based on the generalized elastic theory, which is composed of the total potential energy , total kinetic energy , and external force work of the curved beam structure:

The total potential energy of the structure includes the strain potential energy stored in the structure and the elastic potential energy of the artificial spring.

The relationship between strain and displacement [7] can be expressed aswhere and are the membrane and transverse shear strain on the neutral plane, respectively. is the curvature strain.

The strain potential energy stored in the curved beam structure can be expressed as

In the formula, *E* and *G* represent Young’s modulus and shear modulus, respectively. *A* is the cross-sectional area of the curved beam, is the moment of inertia to the *z*-axis, and *κ* is the shear correction coefficient.

The elastic potential energy stored in the artificial spring group at the boundary of the curved beam is

The total kinetic energy of the structure is composed of the kinetic energy stored in the beam structure and the kinetic energy of the additional mass points that may exist.

Structure’s kinetic energy is as follows:where is the material density of the curved beam.

Additional mass point kinetic energy is as follows:where is the position coordinate of the concentrated mass point on the curved beam.

External force work is done by concentrated load where is the position coordinate of concentrated load on the curved beam.

Before using the variational principle to obtain the final vibration equation of the system, the displacement tolerance function needs to be clarified.

The differential equation of motion for the in-plane vibration of a curved beam [7] is as follows:

It can be seen from the above differential equations that the tangential displacement , the normal displacement , and the rotational displacement are required to have a continuous first derivative on the entire beam (including the ends). In this paper, a form of the improved Fourier series proposed by Su et al. [27] will be used to construct the displacement tolerance function. The specific form is as follows:where *λ*_{m} = *mπ/L*_{s}. *ω* is the natural circular frequency of the beam structure. , , , , , and are the unknown coefficients of the improved Fourier series. The supplementary auxiliary function has the following form:

Combine the above formulas and find the extreme value of the energy functional

The generalized coordinate vector is composed of the following form [*A*_{1}, …, *A*_{M,}*a*_{1}, *a*_{2}, …, *c*_{1}, *c*_{2}]^{T}.

The final matrix in-plane vibration equation is as follows:where is the mass matrix and is the stiffness matrix. is the generalized force vector of external force. For free vibration, the external force vector *F*^{in} is zero. The determinant of the coefficient matrix needs to be zero to get the solution of equation (14). Then, the natural frequencies can be determined by calculating the roots of the determinant and the corresponding mode shapes can be identified.

###### 2.2.2. Out-of-Plane Problem

Similar to the in-plane problem, when solving the out-of-plane problem, the spring system shown in Figure 3 will be used to simulate the constraint of the boundary for the three vibration displacement components.

In the figure above, *K*_{u0}, , and are the constraint spring groups at the left end of the curved beam (at the starting point), which, respectively, constrain *u*, rotational displacement, and direction displacement. *K*_{u1}, , and have similar definitions and act on the right end (endpoint) of the curved beam. is the external concentrated dynamic load.

The energy function of out-of-plane vibration Π^{out} is composed of the total potential energy , the total kinetic energy *T*^{out}, and the external force work term *W*_{e}^{out}.

The total potential energy of the structure includes the strain potential energy stored in the structure and the elastic potential energy of the artificial springs.

The out-of-plane bending, shear, and torsional strain components of a curved beam under the curvilinear coordinate system [8] are

The strain potential energy stored in the curved beam structure is

In the formula, and *J*, respectively, represent the *y*-axis moment of inertia and the torsional constant of the cross-section.

The elastic potential energy of the artificial spring group at the boundary of the curved beam is

The total kinetic energy *T*^{out} of the structure is

Structure’s kinetic energy is

Additional mass point kinetic energy is

External force work is done by concentrated load .

The out-of-plane vibration differential equation [29] is

Similarly, an improved Fourier series is also used to construct the displacement allowable function of out-of-plane vibration. The specific form is similar to the previous one

Subsequently, the out-of-plane vibration equation can also be written in the following matrix form:

is also composed of unknown Fourier coefficients in the form of [*D*_{1}, …, *D*_{M,}*d*_{1}, *d*_{2}, …, *h*_{1}, *h*_{2}]^{T}.

Since the in-plane and out-of-plane vibrations of a curved beam are not coupled with each other, the two can be combined in the following form to obtain the matrix vibration equation of the entire system:

#### 3. Numerical Results and Discussion

##### 3.1. Analysis of the Convergence and Accuracy of the Calculation Model

The verification model takes a curved beam in [5], which is an arc of 1/6 part of a full circle. The material and geometric parameters of the arc are as follows: *E* = 2.1 × 10^{11} Pa, the radius of curvature *R* = 1 m, *ρ* = 7800 kg/m^{3}, *L*_{s} = 1.0472 m, the section is circular, the shear correction factor *κ* = 9/10, and the section radius *R*_{b} = 0.01 m. Since the reference only considers in-plane vibration, the calculation results in this paper will also be compared with the FEM calculation results. The number of FEM units is 210.

The order intercepted after the expansion of the improved Fourier series used in this paper has a direct effect on the accuracy of the calculation result. The more the expansion order, the closer to the accurate value, but it will affect the efficiency of the solution. Generally speaking, the calculation result can stabilize after the selected order reaches a certain value. Table 1 shows the calculation results of the verification model under different modified Fourier series cutoff numbers *M* and compares them with the calculation results of reference and FEM. At this time, the verification model adopts the free-free boundary condition, which corresponds to the value of 0 for the boundary spring groups in the model. The comparison content is the natural frequencies after ignoring rigid modes.

It can be seen from Table 1 that the vibration performance calculation model in this paper has good convergence and numerical stability. When the value of *M* is 16 or greater, there is no significant change in the natural frequencies. Therefore, in the subsequent calculations, the number of truncated terms of the displacement function series is *M* = 16. Besides, since both the paper and the FEM element consider the shear effect and [5] uses the Euler beam theory, the results of this paper and the FEM calculation are smaller than those of [5], and the results of the former two are closer to each other.

As mentioned earlier, this paper uses a spring system to simulate the boundary conditions of the curved beam. Therefore, the influence of the value of the springs on the natural frequency of the curved beam structure needs further research. Firstly, define the dimensionless natural frequency of the curved beam as Ω = 100*ωL*_{s}. The research object is still the previous arc, but the values of the 6 springs *k*_{u1}, , , , *K*_{u1}, and will increase from 0 to 10^{14} in turn, and the first three-order dimensionless natural frequencies of the curved beam structure are shown in Figure 4.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

It can be seen from Figure 4 that the greater the spring stiffness, the greater the stored elastic potential energy and the stronger the corresponding displacement constraint at the boundary of the curved beam, but the boundary constraint effect produced by the increase in stiffness will inevitably tend to be completely constrained. Therefore, a spring group with a certain value can be used to simulate various classical boundaries. For example, the fixed-support boundary should be sufficiently large for all spring groups. Secondly, if the value of a certain spring stiffness is not 0, the rigid body mode of the curved beam structure to the corresponding displacement will disappear, which is reflected in Figure 4 that all the 1st modes increase from 0. Besides, for the mode with a dimensionless frequency of 10.28 is an in-plane vibration mode, it will only be affected by the springs *k*_{u1}, , and of the in-plane model; similarly, the mode with a dimensionless frequency of 27.79 is an out-of-plane vibration modal; it will only be affected by out-of-plane springs *K*_{u1}, , and . Reflected in the three pictures *a*, *b*, and *c* in Figure 4, the 3rd mode (out-of-plane mode) remains unchanged when the in-plane springs *k*_{u1}, , and change, and the 2nd mode (in-plane mode) has an obvious trend of “slow increase–rapid increase–steady change” when the in-plane spring changes. The reason why the 2nd mode remains unchanged in the three pictures *d*, *e*, and *f* in Figure 4 is also the same. In particular, as shown in Table 2, even if it is the corresponding out-of-plane vibration, the change of *K*_{u1} has little effect on most natural frequencies of the curved beam structure, so the 3rd mode curve in Figure 4(e) does not change significantly.

In general, although the curves of the natural frequencies of the curved beams are slightly different when the parameters of different springs change, they can be stabilized when the stiffness value increases to about 1 × 10^{12}. Therefore, the rest of the work in this paper considers that when the stiffness value of the spring group is 1 × 10^{12}, the displacement or rotation in the corresponding direction is completely restrained. Table 3 will give several cases of spring simulation boundary.

To further prove the correctness of the simulated boundary hypothesis, the above boundary is combined in pairs. For example, F-E means that one end is free and the other end is elastically supported, and the natural frequency calculation values are compared with the finite element calculation results. The results are shown in Table 4.

The results in Table 4 show that the calculation model in this paper can be applied when springs are used to simulate various classical and elastic boundaries. The calculation results are accurate and reliable and have a wide range of applications.

Then, consider the influence of external dynamic load and mass point. It is still a 1/6 segment arc model, the material and geometric parameters remain unchanged, and the boundary condition is C-F boundary. Add a 5 kg concentrated mass point at the free boundary end, and apply the excitation force of unit amplitude in *u*, *v,* or direction successively there. Compare the method in this paper with the FEM’s forced vibration response under excitation in the same direction. The frequency sweep range is 0.1∼200 Hz, and the step length is 0.1 Hz. Extract the displacement response of the structure in *u*, , and directions at the mass point and the midpoint of the curved beam. Take the reference displacement *D*_{0} = 1 × 10^{−12}m, *D* is the effective value of the displacement response of the extracted point, and then the vibration response displacement level *L*_{D} = 20log_{10}(*D*/*D*_{0}). The comparison of vibration displacement level is shown in Figure 5.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figures 5(a)–5(c) are the displacement responses at the midpoint of the curved beam, and Figures 5(d)–5(f) are the displacement responses at the mass point of the curved beam. First of all, the results obtained by the method in this paper are in good agreement with the results of the finite element calculation under excitation in all directions, which further verifies the correctness of the calculation model in this paper. Secondly, since the vibrations in *u* and directions are both in-plane vibrations and the two are coupled with each other, they can stimulate responses in both *u* and directions when they are excited by a single excitation in *u* or direction. The vibration in the direction belongs to out-of-plane vibration, so the excitation force in this direction will only stimulate the displacement response in the corresponding direction.

Under normal circumstances, it takes a long time to calculate the forced response of the structure using the FEM method. However, the method in this paper requires a small number of truncated terms, and the generated system stiffness and mass matrix are also small, so the calculation cost is lower. Specifically, on the same hardware (a processor equipped with 3.2 GHz Intel Core i5-6500 CPU and 16 GB memory), it takes about 4.2 seconds to calculate 2000 frequency steps using this method, while the time consumed by the finite element calculation is more than 10 times. Besides, the calculation model in this paper does not need to change the number of mesh elements or remesh like the FEM model when the size parameters such as the length of the curved beam structure are changed. A large number of calculations may be needed for parametric analysis of curved beam structures. The method presented in this paper has obvious advantages because of its low cost and wide application range.

##### 3.2. Engineering Application Examples

The above discussion and most of the research objects in the references assume that the curvature, cross-section, and other parameters are constant, but this assumption may not be appropriate in actual engineering applications. The calculation model in this paper is still applicable when facing such problems. Here are a few examples to introduce how to deal with these problems.

###### 3.2.1. Curvature Continuously Changing Structure-Transition Curve (Whirl Line) Structure

In engineering, it is often necessary to connect a linear structure with a circular arc structure. At this time, the transition curve often chooses a relaxation curve whose radius of curvature gradually changes from infinity. The curve requires the product of the arc length under the curve coordinates and the radius of curvature of the point to be a constant, then the whirl line is mostly used as the easement curve.

The basic formula of the whirl line is

In the above formula, *As* is the whirl line parameter. It is usually defined that the starting point of the whirl line *s*_{0} = 0, where the radius of curvature *R*_{0} is infinite; the endpoint *s*_{1} = *L*_{s}, where the radius of curvature *R*_{1} is arbitrary constant.

To prove that the model in this paper can be applied to structures with varying curvatures, this paper will take a whirl line structure with a length *L*_{s} of 1.4849 m and a radius of curvature *R*_{1} of 1 m at the endpoint as an example to compare the natural frequency calculation results of this method and FEM under different boundary conditions. The relevant parameters such as the cross-sectional shape and material of the structure are still the same as the previous model.

From the comparison of the results in Table 5, it can be seen that the calculation model in this paper also has good accuracy when dealing with slowly varying curvature structures.

###### 3.2.2. Cross-Section Mutation Structure

In actual engineering, the structure is often not composed of a single material or subunits with characteristic geometric dimensions. When dealing with this kind of problem, even if the curvature between the subunits remains continuous, it is necessary to model the subunits separately and then build a coupling model based on the continuous displacement conditions between the subunits. This paper will take the arc structure with abrupt cross-sectional dimensions as an example to introduce how to use this method to solve such problems. An example is a semicircular structure composed of two 1/4 arcs as shown in Figure 6, both of which keep the circular cross-section unchanged, and the cross-section radii *R*_{b1} = 0.01 m and *R*_{b2} = 0.02 m. The material parameters are consistent with the previous model, and the radius of curvature is 0.9549 m.

As shown in Figure 6, when processing the continuous conditions between beam 1 and beam 2, this paper will also introduce an artificial spring system similar to the simulated boundary conditions. The elastic potential energy *V*_{c} of the coupled artificial spring group can be expressed as

In the above formula, is the spring that connects the tangential displacements *u*_{1} and *u*_{2} at the coupling point, and the definition of the rest springs is similar to this, so we will not repeat them. *L*_{s1} is the arc length of beam 1.

If it is considered that beam 1 and beam 2 are welded together, that is, the corresponding displacements of the two at the coupling point should be kept strictly continuous, then the coupling springs can be set to a larger stiffness value of 10^{12} at this time. Compare the natural frequency calculation results of this method and FEM under different boundary conditions.

It can be seen from Table 6 that, after introducing the coupled artificial spring group, the model in this paper can accurately calculate the vibration characteristics of the structure with abrupt cross-sections. Besides, if the value of the rotation restraint spring in the coupling spring is set to 0, the articulation in the actual project can be simulated, and the smaller stiffness value of each spring can be used to simulate the elastic connection.

###### 3.2.3. Curvature Mutation Structure: S-Shaped Structure

In actual engineering, structures such as locks and bridges often take a form of sudden curvature, S-shape. Generally, the curvature of the subunits of this structure is equal but the direction of curvature is opposite. This paper will also take the S-shaped structure shown in Figure 7 composed of two 1/4 arcs as an example to show how to use this method to deal with this type of problem. At this time, the radius of the arc section is *R*_{b1}*=R*_{b2}*=*0.01 *m*, and the other parameters remain unchanged.

When dealing with the continuous condition of the S-shaped structure, the elastic potential energy of the coupling spring will become

Also, compare the natural frequency calculation results of the method in this paper and the FEM under different boundary conditions for the welding connection.

As shown in Table 7, the calculation model in this paper can also well calculate the curvature mutation (S-type) structure. Similarly, different spring values can also simulate different connection conditions.

#### 4. Conclusion

Based on the principle of energy variation, this paper establishes a vibration analysis model of curved beam by introducing an improved Fourier series as an allowable function of vibration displacement. Based on the generalized shell theory, considering shear and inertial effects and coupling between displacement components, the kinetic energy and strain potential energy of the curved beam structure are obtained. Subsequently, the artificial spring system is used to simulate the possible constraint conditions at the boundary of the curved beam, and the constraint effect is converted into the corresponding elastic potential energy. The spring system is suitable not only for elastic boundaries but also for simulating various classical boundary conditions. At the same time, by introducing the corresponding energy term into the energy function of the entire curved beam system, it is convenient to consider the concentrated mass point or the concentrated point excitation force at any position. Compared with the common (semi)analytic method, the calculation model constructed by the improved Fourier series and the artificial spring system is a unified model with wider applicability.

Then, the comparison with the FEM calculation results of the curved beam model under different boundary conditions shows that the calculation model in this paper has good accuracy. Besides, due to the good convergence of the improved Fourier series, the dimensions of the stiffness matrix and the mass matrix formed by the method in this paper are small. The method in this paper has higher efficiency when calculating the forced vibration, which makes it easy to use in the parameters’ study.

Finally, the common structures in engineering such as the curvature easement curve (whirl line) structure, the cross-section abrupt structure, and the S-shaped structure are used as examples to show the wide application prospect of the method in this paper.

#### Data Availability

The data used to support the findings of this study are included within the article and available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (contract nos. 51839005 and 51879113) and the Chinese Fundamental Research Funds for the Central Universities (no. 2019kfyXMBZ048).