Abstract

In the modeling of railway vehicle-track dynamics and wheel-rail damage, simplified tangential contact models based on ellipse assumption are usually used due to strict limitation of computational cost. Since most wheel-rail contact cases appear to be nonelliptic shapes, a fast and accurate tangential model for nonelliptic contact case is in demand. In this paper, two ellipse-based simplified tangential models (i.e., FASTSIM and FaStrip) using three alternative nonelliptic adaptation approaches, together with Kalker’s NORM algorithm, are applied to wheel-rail rolling contact cases. It aims at finding the best approach for dealing with nonelliptic rolling contact. Compared to previous studies, the nonelliptic normal contact solution in the present work is accurately solved rather than simplification. Therefore, it can avoid tangential modeling evaluation affected by inaccurate normal contact solution. By comparing with Kalker’s CONTACT code, it shows both FASTSIM-based and FaStrip-based models can provide accurate global creep force. With regard to local rolling contact solution, only the accuracy of FaStrip-based models is satisfactory. Moreover, Ayasse-Chollet’s local ellipse approach appears to be the best choice for nonelliptic adaptation.

1. Introduction

Wheel-rail damage has been a main concern of railway industry for increasing maintenance cost and potential risk. It is of importance to develop a reasonable maintenance strategy. Therefore, an efficient contact model is needed for computation of both vehicle-track interaction and wheel-rail rolling contact. The term ‘efficient’ here means the contact model is capable of solving nonelliptic rolling contact accurately and fast, because wheel-rail damage simulation such as wear prediction always needs millions of contact cases to account for the actual operating distance. To this end, exact contact models based on finite element method [19] are often not feasible due to their relatively high computational cost. In general, the contact problem can be treated in normal and tangential direction independently since the wheel and rail have similar material properties. This paper focuses on simplified tangential wheel-rail contact modeling.

The tangential wheel-rail interaction depends on rolling contact for low energy dissipation. Rolling contact problem was first studied by Carter [10], who regarded the wheel as a cylinder, thus simplifying this problem as a 2D contact case with longitudinal creepage. Using Carter’s theory, tangential stress, stick-slip division and creep force curve are exact for the contact case whose lateral semiaxis length is much larger than that in longitudinal direction. Johnson [11] extended Carter’s theory to the 3D elliptic contact case containing longitudinal and lateral creepage. Rolling contact considering arbitrary geometry and the combination of creepage and spin is successfully solved by Kalker’s ‘exact theory’ that is implemented in a well-known code named CONTACT [12]. Knothe [13, 14] proposed a more efficient contact model by considering contact zone as several strips along rolling direction instead of elements used in Kalker’s method. Zhao and Vollebreqt [15] accelerated Kalker’s method by using nonlinear conjugate gradient method and fast Fourier transform technique. Similar researches can also be found in [16]. However, computational efficiency still limits the application of these exact models [1215] to on-line wheel-rail damage simulation.

The most widely used tangential contact model for fast calculation is Kalker’s simplified theory and its algorithm FASTSIM [17]. It introduces Winkler’s assumption to relate the displacement and contact stress in linearity and then contact stiffness can be expressed by three groups of springs. Compared to other fast contact models [1820], FASTSIM can also provide detailed rolling contact solution such as tangential stress, stick-slip division and slip velocity. Therefore, FASTSIM is favored by many engineering models [21, 22]. Despite providing satisfactory creep force, FASTSIM is considered less accurate on local tangential stress and slip velocity distribution due to inaccurate assumptions introduced, e.g., linear tangential stress distribution and parabolic traction bound. Recently, an alternative algorithm named FaStrip was proposed at KTH Royal Institute of Technology, which is considered more accurate than FASTSIM in terms of tangential stress distribution [23]. However, FaStrip has not been widely used in the modeling of vehicle-track interaction and wheel-rail rolling contact. In other words, the robust of it should be further verified.

It is noted that both FASTSIM and FaStrip are developed based on the assumption of elliptic contact shape, thus triggering a struggle for many scholars to adapt them to nonelliptic cases. In the literature, three alternative nonelliptic adaptation approaches originally developed for FASTSIM are available [2426]. This is while they have not yet been applied to FaStrip. The core of these three approaches is modifying the flexibility parameters used in FASTSIM by equivalent ellipse [24] or local ellipse [25, 26] to replace the original ones derived from elliptic shape. The tangential performance of these three approaches for nonelliptic contact cases was compared and verified by [2730].

However, it is worth mentioning in [2730] that the normal contact is solved by different simplified models that cannot provide an exact solution for nonelliptic contact. Since tangential traction boundary depends on normal pressure distribution and flexibility parameters are determined based on contact shape, the estimation errors made by simplified contact models in normal modeling may compensate the estimation errors in the tangential solution. It inevitably remains a question whether the performance of tangential contact models is due to tangential algorithms themselves or the combination of simplified normal and tangential modeling. Therefore, the scope of the present work is to perform an investigation on which nonelliptic adaptation approach in [2426] is more accurate for simplified ellipse-based tangential models (i.e., FASTSIM and FaStrip) to use for nonelliptic contact case.

In this paper, both FASTSIM and FaStrip are adapted to nonelliptic contact condition using three available approaches proposed in [2426]. The novelty of this paper is the nonelliptic normal contact solution is accurately solved by Kalker’s NORM algorithm [12] rather than simplification as done in [27, 28]. Thus, it helps making a comparative study not affected by inaccurate normal modeling. The performance of these simplified algorithms is evaluated by Kalker’s CONTACT under typical rolling contact cases using standard wheel-rail profiles. Such a comparative study contributes to a fast and accurate contact model applied to fast calculation of vehicle-track dynamics and wheel-rail rolling contact considering nonelliptic contact.

The outline of this paper is as follows. Section 2 makes a brief introduction of the methods used in this paper and describes our methodology to evaluate the performance of simplified models. Section 3 presents a detailed comparative study to find the best choice for nonelliptic adaptation using simplified tangential models. Section 4 closes with conclusions.

2. Methods Description

2.1. General

A methodology aiming at evaluating the performance of simplified tangential contact models on nonelliptic cases is presented in Figure 1.(i)In previous comparative studies [2730], the normal modeling used simplified contact models, namely, KP, Linder, STRIPES, and ANALYN. Their normal results were only similar to the exact solution and then the inaccurate traction bound (the product of normal pressure and friction coefficient) could affect the prediction of stick-slip division and tangential stress distribution in tangential modeling. To address this issue, our methodology employs Kalker’s NORM algorithm [12] to obtain exact normal solution including contact patch and normal pressure distribution. Note that although NORM is computationally expensive compared to simplified contact models developed in [2426], it is still accepted for on-line vehicle dynamics calculation [31]. In this paper, standard wheel-rail profiles with different lateral displacements are considered to model actual nonelliptic contact conditions.(ii)Both FASTSIM and FaStrip are employed to evaluate their performance on predicting global creep force and local contact solutions such as tangential stress distribution, stick-slip division, slip velocity, and wear distribution. Three alternative approaches proposed in [2426] are used to adapt the two tangential models to nonelliptic contact cases. All of them will be described in Sections 2.22.4.(iii)Kalker’s CONTACT code is employed as a reference to evaluate the simplified contact models under the combination of creepage and spin.

2.2. Contact Models
2.2.1. Kalker’s NORM Algorithm for Normal Wheel-Rail Contact

In Kalker’s NORM algorithm, the wheel-rail contact is locally approximated by elastic half-space. The potential contact zone is discretized by rectangular elements as shown in Figure 2. In every element, the normal contact satisfies (1),where δ is penetration; h is the undeformed distance that can be solved via a rigid geometrical method [32, 33]; u(x, y) is the elastic deformation within the contact patch Ac and can be obtained from half-space equation, where AIJ is the influence coefficient in normal direction, meaning the displacement in element I subjected to a unit load in element J; p is the normal pressure in element J; v and E are material’s Poisson’s ratio and elastic modulus, respectively.

From (2), we can see elastic deformation depends on pressure within contact patch, and the determination of contact boundary also needs to consider elastic deformation. Therefore, there is no analytical solution for (1). The final solution can only be achieved by iteration until pressure of any element within the contact patch meet (3), This iteration process is solved using an optimization procedure by Kalker [12].

2.2.2. Kalker’s FASTSIM Algorithm for Tangential Wheel-Rail Contact

Kalker’s FASTSIM is developed based on his simplified theory [17]. It originally assumed the contact patch to be an ellipse composed of some strips along the rolling direction; see Figure 3(a). The wheel-rail contact is assumed to be replaced by a set of springs with flexibility parameters L1, L2, and L3 in longitudinal, lateral, and normal directions as in (4), where c11, c22, and c23 are Kalker’s creep coefficients depending on semiaxes ratio a0/b0; G is shear modulus; m and n are Hertz coefficients. In this way, linear relation between elastic deformation and tangential stress is supposed. Increment of tangential stress (δqx and δqy) in each strip can be estimated from the leading edge where the tangential stresses (qx and qy) are set to be zero, where υx, υy, and are longitudinal, lateral creepage and spin, respectively. The traction bound is assumed to follow parabolic distribution, where μ is the coefficient of friction; p0 is the maximum pressure. Note that the parabolic traction bound is unrealistic but necessary in FASTSIM to ensure a sufficient exact creep force.

If , the point at (x-δx) locates in the stick area and its tangential stress becomes Otherwise, the point is determined to be in the slip area and its tangential stress isWithin the slip area, the relative slip can be determined from (9), It is noted that the first two terms on the right side are rigid slip and the last part denotes elastic slip.

2.2.3. KTH’s FaStrip Algorithm for Tangential Wheel-Rail Contact

In FaStrip algorithm, the contact patch is also discretized as strips that in FASTSIM; see Figure 3(b). The half-length of slip area in each strip is determined from where ξ, η, and ψ are nondimensional terms as a function of creepage and spin. where p0 is the maximum pressure. The longitudinal and lateral shear stress in stick area can be obtained from where κ, , λ, and are dependent on nondimensional terms ξ, η, and ψ. In slip area, the total tangential stress equals traction bound qt (x, y), The components of tangential stress in x- and y-direction are determined via FASTSIM and relative slip is determined by where ∂ux/∂x and ∂uy/∂x, i.e., the elastic slip, are derivative of longitudinal and lateral tangential stress along running direction, respectively. For more detailed description of FaStrip, the readers are referred to [23].

2.3. Nonelliptic Adaptation Approaches

In this part, it describes three nonelliptic adaptation approaches proposed by Kik-Piotrowski [24], Linder [25], and Ayasse-Chollet [26], which are schematically illustrated in Figure 4.

2.3.1. Kik-Piotrowski’s Approach

In Figure 4(a), the nonelliptic contact patch is regarded as one equivalent ellipse and the three flexibility parameters used in FASTSIM are kept constant throughout the contact patch. The equivalent ellipse is defined, where ae and be are semiaxes length of the equivalent ellipse; L and W are length and width of the nonelliptic contact patch; β is expressed by an empirical formula depending on L/W. Kalker’s three flexibility parameters are then redefined, where A is the area of contact patch; yL and yR are left and right boundary of the contact patch; Rc and h(y) are corrected radius and interpenetration function, respectively.

2.3.2. Linder’s Approach

The idea of Linder’s approach is to apply the tangential model on every strip of contact patch as shown in Figure 4(b). These strips are determined by respective virtual ellipses that share the same width 2be as defined in (20), Then, we define a new coordinate system whose lateral axis is obtained by shifting be and longitudinal axis keeps the same as the original one, From the ellipse equation, the semiaxes length of virtual ellipse isFor each strip, the tangential model is applied locally with respective flexibility parameters, It should be noted in Linder’s approach, the parabolic traction bound used in FASTSIM is changed to

2.3.3. Ayasse-Chollet’s Approach

This approach also regards the nonelliptic contact as several strips along the rolling direction. Each strip is assumed to locate in the center of respective virtual ellipse as shown in Figure 4(c). For each virtual ellipse, the semiaxis length in x-direction is the longitudinal contact boundary of each strip, and lateral semiaxis length is determined by local curvature. Therefore, it needs a reliable smooth procedure to filter lateral combined curvature. Otherwise, negative curvature or discontinuous curvature can significantly affect contact shape and pressure. In [26], it proposed a smoothing strategy based on Boussinesq approach. where l0 is characteristic length; B and Bf are original and filtered lateral curvature. In each y-coordinate, Hertz coefficients m(y) and n(y) can be obtained through a predefined table according to an angle, Thus, the flexibility parameters of this approach are drawn from (27),

2.4. Wear Simulation

For wear calculation, the wear model proposed by KTH Royal Institute of Technology [34] is used. It assumes wheel-rail rolling contact process keeps constant. The wear depth at one position after one wheel passing over can be calculated by integrating the wear contribution in x-direction. where δ(x, y) is the wear depth at each element and can be expressed bywhere k(x, y) is the wear rate associated with pressure p(x, y), running speed v, relative slip s(x, y), and hardness of wheel-rail material H as shown in Figure 5. For a detailed explanation of this method, the readers are referred to [28, 34].

3. Results and Discussion

3.1. Normal Contact Modeling and Cases Study Illustration

In order to evaluate the performance of aforementioned six simplified models, it will present a comparative study of wheel-rail rolling contact under the combination of creepage and spin. Standard wheel-rail profile of S1002CN-CN60 is considered, with prescribed lateral displacement, yaw angle, friction coefficient, and rail cant. A further study for worn profiles will be investigated in the next work. Note that before performing a contact simulation; it should carry out a contact geometry calculation [32, 33]. The input parameters for contact geometry and wheel-rail’s material property are listed in Table 1.

In every lateral displacement case, the input penetration is obtained using Hertz theory under the prescribed load of 83.3 kN [28]. Such an input mode is used to resemble on-line calculation of wheel-rail rolling contact. Figure 6(a) illustrates the variation of contact shape and pressure under 10 shift cases (ranging from -3 mm to 6 mm at 1 mm increment) calculated by NORM. For visualization, the central position of each contact patch is shifted along the longitudinal direction while the lateral position still remains its true location. With the increase of lateral shift, the contact patch locates towards rail gauge and wheel flange root. The detailed normal solutions for the 10 shift cases are graphically presented in Figure 6(b).

For including lateral creep behavior, specified yaw angles are applied to the ten shift cases. It means the yaw angle (being defined as 1.2 times the lateral displacement) increases from -3.6 mrad to 7.2 mrad in 1.2 mrad increments. Since the yaw angle is not enough large to change contact shape obviously, we only consider its effect on creepage. For a detailed investigation of yaw angle’s role on contact shape and creepage, readers are referred to [33]. At each case, the creepage (cx and cy) and spin () are calculated by (30), where ψ, , and δ are yaw angle, roll angle, and contact angle, respectively; R and R0 are rolling radius and nominal radius, respectively.

The variation of longitudinal, lateral creepage, and spin with lateral shift are illustrated in Figure 7 (detailed data are available in Table 2). The longitudinal creepage decreases with the increase of lateral shift. This variation can be explained from (30) that the longitudinal creepage is mainly affected by rolling radius difference (i.e., the first term of right side in cx). For lateral creepage, it changes linearly since lateral creepage can be considered as yaw angle when it is small. The variation of spin is similar to that of longitudinal creepage because the determination of spin mainly depends on contact angle whose variation is close to rolling radius difference.

For the sake of simplicity, three typical nonelliptic contact cases under different shifts (0, -2 and 6 mm), creepage, and spin are selected to evaluate detailed tangential solutions including ‘tangential stress and stick-slip division’ in Section 3.2.1 and ‘slip velocity and wear distribution’ in Section 3.2.2. Figure 8 shows the tangential stress distribution and stick-slip division (the slip zone is encircled by red lines) using CONTACT for the typical cases, which are selected as a reference for Section 3.2.1. In these results, the arrows point the direction of tangential stress and their lengths are proportional to tangential stress’ magnitude. Furthermore, all shift cases in Figure 6 are considered to estimate ‘creep force’ in Section 3.2.3.

3.2. Comparative Study
3.2.1. Tangential Stress and Stick-Slip Division

Figure 9 compares tangential stress distribution and stick-slip division predicted by different models for the case wheelset locating on the central position and the exact solution is referred to Figure 8(a). Detailed tangential stress distribution occurring at y = 0 is shown in Figure 10. For simplicity in plotting results, the three nonelliptic approaches stated in Section 2.3 are termed as KP, Linder, and AC, respectively.

Compared to the reference, it is found FASTSIM-based models result in different tangential stress distribution. This is because the tangential stress in FASTSIM increases linearly in the adhesion area until it saturates in the slip area and follows a parabolic distribution. However, in FaStrip-based models, the stress distributions in the adhesion area are nonlinear and follow the elliptic traction bound in the slip area, which are similar to that in CONTACT. Among them, AC’s approach performs the best.

For the case Δy = -2 mm under a moderate creep level, the tangential stress distribution by simplified models is presented in Figures 11 and 12 and the reference can be found in Figure 8(b). Unlike case Δy = 0 mm where the adhesion zone dominates, slip zone in this case is almost the half of the contact patch. Figures 11 and 12 clearly show a good agreement between FaStrip-based models and the reference. For FASTSIM-based models, their maximum tangential stresses are found to coincide with CONTACT. It is owing to the combination of smaller predicted slip zones and sharper increase of tangential stress compared to the reference (see the first row of Figure 12).

Full slip is observed in the case Δy = 6 mm according to the reference in Figure 8(c) and this phenomenon is successfully captured by simplified models except two approaches of KP and Linder using FASTSIM, see Figures 13 and 14. Higher maximum tangential stress in FASTSIM-based models is found due to unrealistic parabolic traction boundary employed. It is mentioned that the traction bound in Linder’s approach is lower than the ones in KP’s and AC’s approaches, see (24) and (6), respectively.

3.2.2. Slip Velocity and Wear Distribution

Figure 15 shows relative slip distributions along longitudinal direction (at y = 0) for the three shift cases in Section 3.2.1. In the case Δy = 0 mm, results predicted by FaStrip-based models agree with the reference while those using FASTSIM significantly overestimate the relative slip. Moreover, it is found the relative slip calculated by FASTSIM-based models increase linearly while it shows a nonlinear change in the reference and FaStrip-based models. This phenomenon also exists for the case Δy = -2 mm although the errors of FASTSIM-based models are reduced. In the last case, only result using both FaStrip and AC’s approach coincides with the reference.

It is mentioned that the rigid slip contained in the relative slip for all models are identical while the elastic slip contributes to the difference. The discrepancy of elastic slip is influenced by different flexibility parameters. This influence is evident for case Δy = 6 mm where derivative of deformation is higher than those of the other two cases.

Then, the capability of simplified models on simulating wear distribution is explored. From (28) and (29), relative slip is found to play an important role on determining wear rate and then wear depth. Therefore, wear distribution also imply the performance of relative slip occurring at all strips along the longitudinal direction. This is while Figure 15 only presents the result belonging to one strip.

The wear distributions under 200 km/h predicted by different models are illustrated in Figure 16. Here, the wear depth is normalized by the maximum value predicted by CONTACT and such normalization facilitates comparison. Significant errors caused by FASTSIM-based models are observed in the case Δy = 0 mm. It is attributed to the obvious discrepancy of relative slip as shown in Figure 15(a). The higher relative slips calculated by FASTSIM-based models make the wear rate (locating at k2 zone in Figure 5) much larger than those in the other methods (whose wear rate is k1 in Figure 5) and hence a more severe wear can appear. By contrast, FaStrip-based models coincide with the reference due to good agreement in relative slip.

The difference of wear depth using different models is reduced in the case Δy = -2 mm. It can be explained wear rate locates at the same zone in Figure 5 for similar relative slip and then the different wear depth is solely affected by the discrepancy of relative slip. For the case Δy = 6 mm, wear depth by all models agree with each other very well because their average relative slip are almost the same.

3.2.3. Creep Force

In this part, we evaluate the performance of simplified models on predicting creep force, which is very important for vehicle dynamics simulation such as vehicle negotiating a curved track [29] and turnout [35, 36] as well as modeling track vibration [37, 38].

Figure 17 illustrates creep force estimation by different models under lateral shift and mixed creepage and spin. Here, the creep force is normalized by the normal force and friction coefficient. Generally, all models agree with the reference and FaStrip-based models among them tend to be more accurate than FASTSIM-based models.

In order to quantify the performance of these models for the calculation of creep force, we define an average error in (31),where i (= 1, …, m) means the number of shift cases; j (= 1, 2, 3) means normalized longitudinal, lateral, and resultant direction; Fijs and Fijc are predicted normalized creep force by simplified models and CONTACT for every shift case.

The average errors caused by the six simplified models are compared in Figure 18. These errors range from 6.27% to 11.99%, showing a satisfactory capability for predicting creep force. The main contribution to average error is the case Δy = 0 mm where the longitudinal creep force is very small and thus the error is exaggerated. It is believed such an error can hardly affect vehicle dynamic behaviors simulation. As expected, AC’s approach is the best selection for nonelliptic adaptation from a creep force viewpoint.

3.3. Discussion

The present work is an extension of our previous study [28], answering the question proposed in that paper, i.e., which approach is the best for nonelliptic adaptation. From the comparative study of tangential stress, stick-slip division, slip velocity, wear distribution, and creep force presented in Section 3.2, we can find AC’s approach is the best choice. It is obviously reflected in the usage of FASTSIM because it uses flexibility parameters (a function of semiaxes ratio) to express the linear increase of shear stress in the adhesion zone.

For instance, in Figures 9 and 10, the contact patch shape is far from an ellipse as shown in Figure 19(a). The flexibility parameters used in Linder’s approach is higher than the other two approaches. Therefore, the increase rate of its linear tangential stress in the adhesion area is lower. This slow rate of stress build-up delays the saturation and hence the formation of slip area. By contrast, AC’s approach of varying flexibility parameters throughout the contact patch is more close to the nature of contact.

On the other hand, in the case Δy = -2 mm, the predicted results using three nonelliptic adaptation approaches are similar. It is due to the fact that, in the case Δy = -2 mm, the contact patch is similar to an ellipse as shown in Figure 19(b). To this end, KP’s approach using constant flexibility parameters can also contribute to a good result.

By contrast, FaStrip is not as sensitive to nonelliptic adaptation approaches as FASTSIM. This is because in FaStrip, the expressions of stick zone and tangential stress are nonlinear, and they cannot significantly be affected by nonelliptic adaptation approaches as shown in (10)-(15).

It is also found though the total tangential stress using FaStrip is very accurate, the component of it in x- and y-direction deviates a bit from the reference as presented in the second rows of Figures 10, 12, and 14. This is because their components ratio are determined by FASTSIM that are inaccurate (see the first rows of Figures 10, 12, and 14). One improvement approach may be the introduction of influence coefficient combined with the available stick-slip division by FaStrip. This improvement is within the scope of our going on work named INFCON.

Even so, the simplified models in this paper are considered as promising solutions to be used in vehicle-track dynamics and wheel-rail damage simulation, because they are both computationally efficient. FaStrip is a bit slower than FASTSIM since it uses FASTSIM to identify the component of analytical total tangential stress in slip zone. In addition, the nonelliptic adaptation approaches will not markedly increase their computational cost since they can be solved explicitly. In our DELL Precision T7910 workstation with the processor of Intel Xeon E5-2630 v3 @ 2.4 GHz, the average computational time is about 2.4 ms and 1.7 ms for FaStrip and FASTSIM, which are about 1000 times faster than CONTACT.

4. Conclusions

In this paper, two simplified tangential models and three nonelliptic adaptation approaches are firstly recalled. They are then combined with Kalker’s NORM to carry out a comparative study of wheel-rail rolling contact and evaluated by Kalker’s CONTACT.

All simplified models can provide satisfactory global creep force, meaning all of them are applicable to vehicle-track dynamics modeling. This is while only FaStrip-based models can capture local nonlinear distribution of tangential stress and slip velocity, contributing to more exact prediction of wear distribution.

By using Ayasse-Chollet’s local ellipse approach, it results in the best prediction compared to the reference. It implies varying flexibility parameters throughout the contact patch are more close to the nature of contact rather than constant flexibility parameter employed in Kik-Piotrowski’s approach. This approach can be further applied to the extension of other ellipse-based method, such as thermal wheel-rail contact analysis [39].

FaStrip is found not very sensitive to the selection of nonelliptic approach. This characteristic shows Kik-Piotrowski’s and Linder’s approaches are also alternative choice and they are expected to be more suitable for worn wheel-rail profile cases, in which the nonsmoothing curvature may affect the accuracy of Ayasse-Chollet’s approach [40]. This point should be further investigated and clarified in the future work.

Data Availability

All data of the current study are available from the corresponding author on reasonable request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The present work is supported by National Key R&D Plan of China (grant no. 2016YFB1102600), National Natural Science Foundation of China (51605318, 51608459, 51778542, and U1734207), Fundamental Research Funds for the Central Universities (2682018CX01), and Jiangsu Provincial Natural Fund project (16KJB580008).