Effects of total pressures and equivalence ratios on kerosene/air rotating detonation engines using a paralleling CE/SE method

Fang Wang, Chun-sheng Weng, Yu-wen Wu, Qiao-dong Bai, Quan Zheng, Han Xu

National Key Laboratory of Transient Physics, Nanjing University of Science and Technology, Nanjing, 210094, China

Keywords:Kerosene/air rotating detonation waves Total pressure Equivalence ratio Instabilities MPI+OpenMP

ABSTRACT In this paper, the kerosene/air rotating detonation engines(RDE) are numerically investigated, and the emphasis is laid on the effects of total pressures and equivalence ratios on the operation characteristics of RDE including the initiation, instabilities, and propulsive performance. A hybrid MPI + OpenMP parallel computing model is applied and it is proved to be able to obtain a more effective parallel performance on high performance computing(HPC) systems. A series of cases with the total pressure of 1 MPa,1.5 MPa,2 MPa, and the equivalence ratio of 0.9, 1, 1.4 are simulated. On one hand, the total pressure shows a significant impact on the instabilities of rotating detonation waves. The instability phenomenon is observed in cases with low total pressure(1 MPa)and weakened with the increase of the total pressure.The total pressure has a small impact on the detonation wave velocity and the specific impulse. On the other hand, the equivalence ratio shows a negligible influence on the instabilities, while it affects the ignition process and accounts for the detonation velocity deficit. It is more difficult to initiate rotating detonation waves directly in the lean fuel operation condition. Little difference was observed in the thrust with different equivalence ratios of 0.9,1,and 1.4.The highest specific impulse was obtained in the lean fuel cases,which is around 2700 s.The findings could provide insights into the understanding of the operation characteristics of kerosene/air RDE.

1. Introduction

Rotating detonation engine (RDE) is a novel engine based on detonation combustion.Compared with deflagration,detonation is characterized by lower entropy production and faster propagation speed,which theoretically provides higher thermal efficiency[1].In recent years, RDE has been studied by numerical simulations and experiments in different aspects such as the thermodynamic performance [1-3], detonation instabilities [4,5], ignition process[6-11],fuel injection[12-17],hollow combustor[18],propagation mode[19,20],turbine combined RDE[21-28],ramjet test[29-31].During this period, the numerical research has demonstrated its ability to guide RDE design and reveal the internal mechanism behind the complex detonation phenomenon.

Based on the aforementioned researches, it can be found that most of the numerical work is focused on hydrogen-fueled RDE.However, due to the merits of convenience and low cost, the kerosene-fueled rotating detonation engine has recently attracted more and more attention from both academic and engineering applications. Bykovskii et al. [32] first experimentally initiated the kerosene-air-oxygen mixture RDE and obtained typical photographic records of continuous detonation waves with the addition of oxygen. They pointed out the initiation of the kerosene-air mixture remains urgent and needs to be further studied. Their following research [33] focused on adding some high-activity fuel such as hydrogen or syngas to obtain the self-sustained rotating detonation waves. The results show bubbling of the gaseous fuel(hydrogen or syngas) through liquid kerosene in the fuel injection system makes it possible to reduce the mass fraction of the gas in the two-phase fuel down to 8.4%for hydrogen and 47%for syngas.Naour et al.[34]proposed the method of pre-vaporized kerosene to improve the reactivity and mixing of reactants. Zhong Yepan et al.[35,36] conducted the experiments on the RDE fueled by cracked kerosene.The mass flow rate of oxygen is decreased to 30%and the effect of the equivalence ratio and oxidizer slot width are analyzed.Kindracki [37] experimentally investigated the initiation and propagation of rotating detonation waves for liquid kerosene and gaseous air mixtures including isopropyl nitrate,gaseous hydrogen,and gaseous air. The above literature presents the latest research results on kerosene/air RDE and illustrates the urgent need for numerical research on ramjet or turbojet kerosene/air RDE. For a propulsive system, the working condition including the total pressure and the equivalence ratio is always important.Wang and his co-workers [38] investigated the effect of injection conditions on the instabilities of H2/Air rotating detonation waves.Zhao et al.[39]reported that under constant pressure conditions the rotating detonation is more prone to instability at low pressure.Meng et al.[40]pointed out that the equivalence ratio plays a very critical role in the formation process and propagation mode by numerical work on a non-premixed RDE.From the above research on hydrogen/air RDE,the total pressure and equivalence ratio significantly affect the operation characteristics of RDE. However, no such work has been done on the kerosene/air RDE. Due to the low-activity nature of kerosene, the effects of these factors on the kerosene/air RDE remain to be studied.

Additionally, one could see that the numerical simulations of RDE have greatly improved contributing to the development of modern high performance computer(HPC) systems. A large scale parallel programming on the HPC systems become more and more important in the numerical simulations.The hybrid nature of these systems - distributed memory across nodes and shared memory with non-uniform memory access within each node - poses a challenge to obtain the high-performance parallel efficiency [41].The dominant parallel strategies using Massage Passing Interface(MPI) alone could not fully utilize the advantage of HPC systems.Previous research[41-44]on parallel efficiency predicts that a hybrid MPI+OpenMP strategy could further improve the parallel performance and lighten the low convergence problem(or even non-convergence) [41] caused by a large amount of domain decomposition. The MPI + OpenMP strategy employs the MPI for communication across nodes at the first level and the OpenMP for parallelization of the computing solver at the second level. This hybrid parallelism model combines shared memory and massagepassing programming better and matches the hierarchical architecture of computer clusters[45].Therefore,it is possible to obtain an advanced paralleling performance to apply the hybrid strategy to the numerical study of RDE.

Overall,we aim to investigate the effects of inlet total pressures and equivalence ratios on kerosene/air RDE by a paralleling CE/SE method to understand the operation characteristics of kerosene/air RDE including the ignition, instabilities, and propulsive performance.This paper has several contributions to existing research.(1)The effects of the inlet total pressures and the equivalence ratios on the operation characteristics of kerosene/air rotating detonation engines are studied; (2) the insights into ignition limits and instabilities of the kerosene/air rotating detonation waves are analyzed; (3) an advanced parallel model performed on HPC systems is implemented and evaluated in RDE simulations.The rest of the paper is organized as follows:In Section 2, the physical model and parallel strategy are presented. The numerical results are presented and discussed in Section 3. Finally, the main conclusions drawn from our work are summarized in Section 4.

2. Physical model and parallel strategy

2.1. Physical model

The two-dimensional Euler equations associated with chemical reaction source terms are used as the governing equations,which is written as

where hiis the enthalpy calculated by the thermochemical relation[46].

For a detailed chemical reaction model, the ith specie chemical reaction production rate ωiis calculated by

In the simulation,a two-step kerosene/air chemical scheme[48]is employed as the source terms and solved by a third-order total variation diminishing(TVD) Runge-Kutta method [49]. The isentropic flow boundary condition [50] is used as the inlet boundary condition and a non-reflecting outlet boundary[51]corresponding to Mach number is employed as the outlet boundary condition.The κ-CNI scheme space-time conservation element and solution element (CE/SE) method [52] is employed to solve the governing equations.

2.2. MPI + OpenMP parallelization

The mesh data is divided and then assigned to each MPI process.The MPI communication takes place at the splitting boundary between two neighboring subdomains. Due to the irregularity and inhomogeneity of unstructured triangular grids,the determination of the communication grids, including the topological relations of the communication grids between two neighboring subdomains,is very difficult. Since one triangle has three neighboring triangles,the grid at the splitting boundary in mesh data j may have one or two neighboring grids located at subdomain j+1.This uncertainty is solved by a judgment sentence and a tag regarding the number and index of neighboring grids located in the neighboring subdomain.For a grid at the splitting boundary of subdomain j,the index of its neighboring grids located at the subdomain j+1 will be assigned to the process j+1 to recognize and send the flow parameters to process j when communication starts. On the other hand, a tag containing the number of neighboring triangles and information of their topological relations will be used by process j to arrange the flow parameters received from process j+1.

Fig. 1 illustrates the two-level parallelization scheme. Each process reads one mesh data and conducts the computing within subdomain.In each process,several OpenMP threads are employed to provide an inner level of parallelization.The CE/SE solver is fully parallelized in each loop. Communication and boundary updating are included in each iteration. Data output is conducted by each process and a reconstruction code of mesh data is developed. This could avoid the long execution time of gathering flow parameters and outputting data by the main process.

2.3. Parallel performance

To evaluate the parallel performance of the hybrid MPI + OpenMP strategy, numerical tests are conducted using different combinations of processes and threads. The tested numerical domain is 0.3 m × 0.1 m containing 1 million triangular meshes. The execution time per iteration is measured by the average time of 1000 iterations using the same code and meshes.The execution time is collected by a build-in Fortran function named system_clock.The employed software is Intel/19.0.0.117 and the platform is ASPIRE1 of National Supercomputing Centre(NSCC)located in Singapore.Each node of the supercomputer contains 24 CPUs.To make the most of each node,the product of the numbers of processes N1and threads N2in each node remains 24. In the N1×N2notation,N1indicates the number of MPI processes in each node, and N2indicates the number of OpenMP threads employed by each process.

The speedup of parallel computing is calculated by

where l, T1, and Tlare the CPU number, execution time by single CPU, and execution time by parallel CPU number of l. The parallel efficiency Elis calculated by

Fig. 2(a) shows the speedup with different CPU numbers. The ideal speedup is also presented to make a comparison.As shown in Fig. 2(a), the speedup shows a big difference between the combinations. The pure OpenMP parallel (1 × 24) obtains the lowest speedup. The combination (4 × 6) obtains the highest speedup,where each node employs 4 processes and each process employs 6 threads. It could maintain the best speedup with the increase of CPU numbers. The speedup of 4 × 6 with 384 CPUs (16 nodes) is 335 times. The mesh independence was validated with 3 million grids and the result shows the combination(4×6)still obtains the best performance with a speedup of 348.54 and an efficiency of 90.77%using 384 CPUs.

Fig.1. Hybrid MPI + OpenMP parallel model.

Fig. 2. Speedup and efficiency of the hybrid parallel model.

Fig. 2(b) shows the parallel efficiency with different CPU numbers.With the increase of CPU numbers,the pure MPI parallel(24×1)shows an obvious decrease in parallel efficiency while the hybrid strategy shows a slight decrease in parallel efficiency. It could be found that the combinations of processes and threads significantly influence parallel efficiency.The combination(1×24)obtains the lowest efficiency. The combination (4 × 6) obtains the best parallel efficiency and could maintain a high efficiency around 90% with the increase of CPU numbers. It performs even better when a large number of CPUs are used compared with other combinations. This phenomenon is a balanced result of the two parallel strategies. Using more OpenMP threads reduces the corresponding number of MPI processes, which reduces the total number of communication messages and increases the total cache size for the same amount of data. Therefore, the overall performance is improved. However, as the number of OpenMP threads increases, the overhead of the OpenMP constructs and remote memory accesses across sockets in a node increases.As a result,we see a relative performance increase up to 6 threads after which we see a performance drop from 6 threads to 24 threads. A similar phenomenon is also observed by Haoqiang Jin et al.[41]in testing a standard CFD problem of Overflow by the hybrid MPI + OpenMP approach.However,they obtained the best performance when the number of OpenMP threads per process is 4 on a supercomputing cluster with 16 cores in each node.

2.4. Grid independence validation and C-J validation

The numerical method and the chemical reaction model are validated with one-dimension detonation waves. The premixed stoichiometric kerosene/air mixture at 1 atm and 300 K is initiated by a hot jet with p = 2 MPa, T = 3000 K, and u = 1000 m/s. The length of the 1D tube is 1 m and the grid sizes contain 0.3 mm,0.2 mm, and 0.1 mm. To remove the influence of the initial speed,the detonation velocity is measured by a mean detonation velocity from 300 μs to 500 μs.Table 1 shows the comparison between the calculated detonation velocity at different grid sizes and the theoretical C-J value of Jet-A computed by Chemical Equilibrium with application(CEA) software [53]. Since this kerosene fuel is not included in any detonation calculation software,the theoretical C-J value of the kerosene-type fuel Jet-A is used here to make a comparison. As shown in Table 1, all three cases obtain a very similar detonation velocity, which verifies the grid independence of the numerical method. The relative error between the theoretical C-Jvelocity of Jet-A and the calculated KERO detonation velocity is roughly 2%. This error may result from the difference between the two fuels.

Table 1 Comparison of calculation and theoretical C-J value.

As reported by Hua Shen [54], the triangle mesh results in an improved resolution in each of the CE/SE schemes because the average area of the triangle elements is less than that of quadrilateral elements.The grid number of triangles is twice the number of 0.2 mm rectangles.From the perspective of the area,the area of a 0.2 mm uniform triangles is equal to that of a 0.141 mm rectangle.Hence, considering the calculation efficiency and accuracy, the triangle grid size 0.2 mm is employed in the simulation.This grid is sufficient to capture rotating detonation waves.

Fig. 3. Detonation velocity validation with different equivalence ratios.

Table 2 Information of simulated cases.

Fig. 4. Ignition results.

Since the effects of the equivalence ratio are considered in this study, the validation between the calculated detonation velocity and the C-J velocity is also conducted. As shown in Fig. 3, the calculated detonation velocity shows a good agreement with the CJ velocity of Jet-A.The deviation is around 2%with the equivalence ratio ranging from 0.8 to 1.4.

3. Results and discussion

3.1. Cases simulated and ignition analysis

The computation domain is 0.3 m×0.1 m containing 1.5 million triangles. A series of numerical simulations are conducted and the emphasis is laid on the effects of inlet total pressures and equivalence ratios on the kerosene/air rotating detonation engines. The information about simulated cases is shown in Table 2. Overall, 9 cases are simulated, and 8 cases succeed to obtain self-sustaining rotating detonation waves. All the cases are performed with a total temperature of 500 K.

The region below y = 0.03 is filled with the stoichiometric kerosene/air mixture with a pressure of 0.3 MPa. To initiate the detonation waves directly,the ignition region[0-0.01]×[0-0.03]is filled with the oxygen-enriched (42% mass fraction) stoichiometric kerosene/air mixture with p = 2.8 MPa, T = 3400 K, and u = 1000 m/s. In this way, the detonation wave could be initiated directly.

Fig. 4 shows the initiation results using the present ignition method. With the increase of inlet total pressures, the required equivalence ratio for successful initiation increases.The cases with the equivalence ratio of 0.9 such as cases 1 and 4 obtained a selfsustained rotating detonation wave, but it failed in case 7 whose inlet total pressure is 2 MPa. A new initiation method where the reactant in the initiation region was replaced with the stoichiometric kerosene/oxygen was tested to initiate case 7. However, it still failed to initiate the rotating detonation waves.Zhao et al.[39]obtained the hydrogen/air rotating detonation waves even with a total pressure of 4 MPa.From the comparison,it could be found that the initiation of kerosene/air rotating detonation waves is much harder.According to the experiments on cracked kerosene/oxygenenriched air RDE conducted by Zhong [35,36], the lean fuel boundary under 30% oxygen-enriched air is around 0.7-1.0 corresponding to the mass rates. Due to the addition of oxygen in experiments,the lean fuel boundary of kerosene/oxygen-enriched air rotating detonation waves would be lower than kerosene/air mixture.

Fig. 5. Initiation process of case 7 and case 8.

Fig. 6. Pressure and temperature history.

The reason accounting for this phenomenon is the low-activity nature of kerosene/air mixture. The lean fuel operation condition will lead to a lower reaction rate, which is even worse when the reactant has a high speed in the axial direction.Besides,the higher density and pressure in the injected fresh mixture region result in the higher compressibility, increasing the required detonation intensity for transitioning to the injected fresh mixture region. As shown in Fig.5,the direct initiation process of case 7 and case 8 are presented to make a comparison. In case 7, the heat release in the reaction zone is not enough to support the post-shock when passing through the high-compressive injected region, which results in the decoupling of the post-shock and reaction zone. The intensity of post-shock decreases continuously and could not block the injection of fresh mixture.Hence the ignition process failed.In case 8, the heat release of the reactant is sufficient to support the post-shock and block the injection of the mixture. The new detonation front has transitioned to the injected mixture and become stable due to the coupling of the post-shock and reaction zone.This also explains why experimental studies needed to add some highactivity fuel or oxidant to enhance the detonability of this mixture.The high-activity fuel or oxidant could provide help in the transition process and enhance the Deflagration to Detonation process.

3.2. Instabilities analysis

Fig. 7. Time evolutions of Case 2(total pressure 1 MPa).

The instabilities in rotating detonation waves are observed and significantly affected by the total pressures in the simulations.From the results, it can be found that all the successful cases could maintain a single wave mode. Fig. 6 shows the time history of the pressure and temperature at the probe near the head end. In each figure, the red line, blue line, and black line represent the equivalence ratio of 1, 1.4, and 0.9 respectively. The rotating detonation waves in these cases all propagate over ten cycles. The pressure peak is coupled with the temperature peak,which reveals the selfsustained propagation of rotating detonation waves.Increasing the inlet total pressure will obtain a higher von Neumann peak pressure. However, the changing in equivalence ratio shows a small impact on the peak pressure. The cases with the equivalence ratio of 0.9 which means RDE operates at lean-fuel condition will obtain a weaker rotating detonation wave and have difficulties in initiation of detonation waves under high inlet total pressures.It can be seen that some fluctuations exist in the pressure curve in the enlarged figure of cases 1-3.These fluctuations become smaller in case 4-6 and almost invisible in case 8-9. This phenomenon indicates that increasing total pressures could decrease the instabilities in the rotating detonation waves. It could be found that the equivalence ratios have a negligible effect on instabilities because the fluctuations have the same amplitude under same inlet total pressures.

Fig. 8. Time evolutions of Case 5(total pressure 1.5 MPa).

Fig. 9. Time evolutions of Case 8(total pressure 2 MPa).

Fig. 7 shows the time evolutions of pressure and temperature contours of case 2 when the rotating detonation wave has propagated over ten cycles.It could be seen that the rotating detonation wave is unstable.The instabilities are induced by the comparatively high pressure and reflected shock wave after the detonation wave[38].The fresh triangular mixture layer is irregular and interrupted by a burned gas bump.From Fig.7(a1)to Fig.7(a2),the detonation height drops quickly because of the burned gas bump.The incident angle between the deflagration surface and the detonation wave decreases.Therefore,the reflecting shock and the slip line all move forward, which leads to the clockwise spinning of the detonation wave structure. As a result, the pressure of the inlet boundary downstream of the detonation wave decreases, and the fresh mixture begins to inject in a long distance approximately from x=0.01 m to x=0.1 m.This simultaneous injection in such a longdistance brings the problem that the fresh mixture will lose a regular triangular shape. As shown in Fig. 7(a3), the detonation wave has passed the burned gas bump and left a high-temperature region at the upper end of the detonation wave. Meanwhile, the detonation front height begins to rise.The changing of the incident angle will force the reflecting shock and slip line to move backward,which leads to the counterclockwise spinning of the detonation wave structure. As a result, a reflecting shock appears at the inlet region downstream of the detonation wave and blocks the injection of the fresh mixture. This block effect interacts with the release effect caused by the drop of the detonation front height,resulting in the appearance of a new burned gas bump as shown in Fig. 7(a5).These features influence each other and cause the instabilities of the rotating detonation waves. Due to the self-adjusting mechanism predicted by numerical simulations[8]and experiments[55],the instabilities will not quench the rotating detonation waves but periodically exist with the detonation waves under certain circumstances.

As can be seen in Fig. 8(a1-a5), increasing the total pressureslightly weakens the instabilities of the detonation wave.A similar phenomenon has been observed in the numerical research of Zhao et al.[39].The contact surface is still distorted and the fresh mixture layer shows an irregular shape. However, the burned gas bump becomes higher than case 2, which reduces the influence of the burned gas bump by adjusting the incident angle when the detonation wave passes it. As shown in Fig. 9(a1-a5), the rotating detonation wave becomes more stable with the increase of total pressure.The contact surface is smooth and the fresh mixture layer maintains a triangular shape throughout all the times.This proves the improvement of increasing the total pressure on the instabilities of rotating detonation waves. The reason account for this is considered as below. First, the intensity of the detonation wave gets higher due to the higher reactant density in the fresh mixture layer.A high-intensity detonation wave is more stable and will be less affected by the downstream perturbations.The second,the high compression in the fresh mixture layer with high pressure so that the fuel mass rates are not greatly affected by the downstream perturbations[39].The third,the burned gas bump tends to be smoothed by the high injection total pressure,which decreases the oscillation of the detonation front height.

Table 3 Propulsive performance.

Fig.10. Thrust and fuel mass rate.

Fig.11. Specific impulse and detonation velocity.

3.3. Propulsive performance analysis

The propulsive performance including the thrust and specific impulse is investigated by the cases with self-sustained rotating detonation waves. The thrust at the chamber outlet is obtained by

where F is integrated along the outlet of the chamber and p∞is the ambient air pressure.The fuel-based specific impulse is calculated by

where g is the gravity acceleration. The fuel mass flow rate on the inlet is defined as

where ρfis the fuel density. Note that these parameters are timedependent quantities. The calculated average propulsive parameters calculated are presented in Table 3. The standard deviation is also presented to show the fluctuations induced by the instabilities.As can be seen,the deviation of specific impulse decreases with the increasing of the total pressure, which predicts the instabilities of rotating detonation waves are weakened.

As shown in Fig.10,the thrust almost increases linearly with the increase of total pressure while the equivalence ratio shows a small impact on the thrust.Even though more fuel is injected,no obvious increase is observed in the rich fuel cases with an equivalence ratio of 1.4.This will greatly decrease the specific impulse of the rich fuel cases. The specific impulse is around 1800-1900 s when the equivalence ratio is 1.4. It could be found that the total pressure shows a small impact on the specific impulse although it will influence the ignition and instabilities of rotating detonation engines.The detonation velocity is highest in the cases with the equivalence ratio of 1. Cases with the equivalence ratio of 0.9 and 1.4 all obtained a lower detonation velocity while the cases with the equivalence ratio of 0.9 obtained the lowest detonation velocity(see Fig.11).

4. Conclusion

The present paper numerically investigated the effects of the total pressures and equivalence ratios on the kerosene/air rotating detonation engines by a paralleling CE/SE method. A hybrid MPI + OpenMP parallel model was applied to accelerate the simulations on the HPC systems. The effects of the total pressures(1 MPa,1.5 MPa, 2 MPa) and equivalence ratios (0.9,1,1.4) on the operation characteristics of RDE including the ignition,instabilities,and propulsive performance are analyzed. The results are as follows:

(1) The instabilities in rotating detonation waves are observed in cases with low total pressure(1 MPa)and weakened with the increase of total pressures. In low total pressure cases, the detonation front height oscillates when passing the burned gas bump. The drop and rise of the detonation front height will result in a release effect and block effect on the injection of the fresh mixture.In this way,a new burned gas bump will be formed.Increasing the total pressure helps to smooth the burned gas bump and reduce the influence from the downstream. This indicates a feasible way to reduce the instabilities is to appropriately increase the total pressure.The thrust increases almost linearly with the increase of total pressure but it obtains a similar specific impulse with different total pressures.

(2) The equivalence ratio shows a negligible influence on the instabilities,but it affects the ignition process and detonation velocity. Due to the low-activity nature of the kerosene/air mixture, the rotating detonation waves are harder to be initiated directly under lean fuel conditions. The highest detonation velocity is obtained with the equivalence ratio of 1, followed by the cases with an equivalence ratio of 1.4.Cases with an equivalence ratio of 0.9 obtain the lowest detonation velocity. The propulsive performance results show little difference in the thrust among cases with the equivalence ratios of 0.9,1, and 1.4.

(3) The hybrid MPI + OpenMP strategy could maintain a high parallel speedup and efficiency with the increase of CPU numbers. The different combinations of MPI processes and OpenMP threads show a significant influence on the parallel performance. In our test, the best parallel performance is obtained with the combination of 4 processes×6 threads in each node. This advanced hybrid MPI + OpenMP parallel model could provide further help for the more complex simulations such as two-phase kerosene/air RDE.

The current work systematically explores the effects of the total pressures and equivalence ratios on the operation characteristics of kerosene/air rotating detonation engines, which is very useful for the understanding of the air-breathing kerosene/air RDE.However,we do not take the liquid phase and some practical factors into consideration. In our future work, numerical simulations on the two-phase kerosene/air RDE with more practical factors will be conducted to further investigate the rotating detonation phenomenon.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

The authors would like to acknowledge the National Natural Science Foundation of China (Grant Nos.11802137, 11702143), the Postgraduate Research & Practice Innovation Program of Jiangsu Province(Grant No.KYCX19_0292),the Natural Science Foundation for Young Scientists of Jiangsu Province of China (Grant No.BK20190468), the Fundamental Research Funds for the Central Universities(GrantNos.30918011343,30919011259,309190112A1). The authors gratefully acknowledge financial support from the China Scholarship Council. The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore(https://www.nscc.sg).