Fully resolved numerical simulations with finite rate chemical reactions and detailed molecular diffusion have been conducted for a series of laminar premixed hydrogen–air flames under atmospheric condition. The objective of the work is to study the influence of unsteadiness of the flow on the local and global flame dynamics. Two equivalence ratios with and 4 are considered, leading to a negative and a positive Markstein number Ma0 at steady-state condition. The flames are excited with oscillating inflows at pre-defined frequencies f to assess the effect of unsteady flame stretch on flame dynamics. The Damköhler number, defined by the ratio of the inverse frequency of the oscillations and flame transit time, is used to characterize the interactions between the flow and the chemical reactions based on their time scales. For both lean and rich flame conditions, the local flame speed Sl is less sensitive to the flame stretch in an unsteady flow, which results in a reduced magnitude of the Markstein number . In addition, is smallest when the time scale of the flow approaches the intrinsic time scale of the flame ( ). The global consumption speed St, computed from integration of the fuel burning rate over the whole computational domain, yields a phase delay and a damped oscillation with respect to the unsteady inflow. The phase delay increases with f or decreasing Da, whereas the reverse trend has been found for the oscillation amplitude of St. The flame is not able to follow the unsteady flow or adjust its flame surface at high excitation frequencies with Da <1, and vice versa in the low frequency range with . An efficiency factor E has been introduced to model the damped response of the flame due to flow unsteadiness, which reproduces the asymptotic behavior of at and at . The simulation results reveal that the fluctuation time scale plays a significant role in elucidating the effect of flame–flow interaction, which should be considered for turbulent combustion modeling.
I. INTRODUCTION
Turbulent premixed combustion is controlled by flame–turbulence interactions leading to changes in local and global flame dynamics and consumption rates. The flame surface is stretched by tangential straining and/or curvature caused by the non-uniform, unsteady flow, altering the internal flame structure as well as the local and global burning rate of fuel. At high pressures, the flame becomes thinner, and highly turbulent fluctuations may lead to flame wrinkles with extremely large curvatures and flame stretch. Accordingly, the classic flamelet approach, which employs stationary laminar unstretched flames for representing the turbulent flame structure, is not suited in this case.
The flame dynamics in terms of the correlation between laminar and turbulent flame speed with flow stretch has been extensively studied in the last few decades.1–10 Most of these works were focused on analyzing single snapshots (e.g., for spherically expanding flames) or time-mean flame structures (e.g., for Bunsen-type flames), which assign the influence of flame stretch to integral properties such as the turbulent length scale Lt or the root mean square (rms) values of the velocity fluctuations . Hence, the time history of the unsteady flow and its effect on the flame dynamics have been neglected. However, as shown in Ref. 11 by detailed numerical simulation of oscillatory plane-jet H2/air flames, the flame displacement is attenuated strongly in an unsteady flow, particularly in the high-frequency range of fluctuations. The flame requires a relaxation time to adapt to the imbalance of the prevailing diffusion–reaction processes caused by the unsteady flow. Similar behavior has been reported in Refs. 8, 9, and 12, applying direct numerical simulation (DNS) of 1D counter-flow oscillating, premixed H2/air flames, and 3D turbulent methane/air flames.1 It has been shown in the latter case that motion of the flame surface was weakened at fluctuation frequencies above the inverse of the flame transit time that is defined by the ratio of the unstretched laminar flame thickness and unstretched laminar flame speed . Therefore, the time scale of the flow represents an essential parameter for studying the effect of flame–turbulence interactions, which should be considered in addition to Lt and . For studying the effect of the flow time scale on flame–turbulence interactions, it is regarded as the oscillation period of the velocity fluctuation and is normalized by the chemical time scale to obtain the Damköhler number .
Although implicitly including an integral turbulent time scale via , present models are not able to predict the observed effect of flow unsteadiness on flame dynamics due to incomplete consideration of the underlying flow history. Instead, most of the available combustion models assume the fast chemistry limit with , i.e., the instantaneous adaption of the chemical reactions to the unsteady flow. Therefore, the objective of the current work is to explore the effect of unsteady stretch on the local and global flame dynamics, which could be incorporated into general combustion modeling concepts. For that purpose, detailed numerical simulations have been conducted for a variety of premixed Bunsen-type, hydrogen/air flames, being excited with harmonically fluctuating inflow velocities at prescribed frequencies. The oscillating inflow is employed to emulate the unsteady nature of turbulent flow. Compared to the previous 1D studies conducted for 1D purely strained counter-flow flames,9 this work assesses quantitatively the local and global flame dynamics by means of the time-mean Markstein number and the global consumption speed. Moreover, the study is conducted for a more general 2D Bunsen-type flame setup, which covers both the tangential straining and curvature effects. In this way, correlations of the local and global flame speed with unsteadiness of the flow can be quantitatively assessed. The results can be used to derive first-order estimations for modeling the time history effect of the flow on premixed flame propagation.
The work is organized as follows: Sec. II presents basic parameters used in this work for characterizing flame dynamics in terms of flame speed and stretch. The numerical setups used for the simulations are shown in Sec. III. The simulation results are discussed in Sec. IV with respect to the local and global flame dynamics. Finally, the main conclusions are summarized in Sec. V.
II. THEORETICAL BASICS
III. SIMULATION SETUP
A. Computational grid
Flame parameters from solutions of 1D H2/air flames.
. | (m/s) . | (mm) . | τc (ms) . | fc (Hz) . | . |
---|---|---|---|---|---|
0.5 | 0.54 | 0.42 | 0.78 | 1287 | 0.0018 |
4 | 1.72 | 0.44 | 0.26 | 3912 | 0.081 |
. | (m/s) . | (mm) . | τc (ms) . | fc (Hz) . | . |
---|---|---|---|---|---|
0.5 | 0.54 | 0.42 | 0.78 | 1287 | 0.0018 |
4 | 1.72 | 0.44 | 0.26 | 3912 | 0.081 |
As illustrated in Fig. 1, the numerical simulations have been conducted for a Bunsen-type flame. The computational domain is built from a wedge with an opening angle of 3°, assuming axial symmetry for a laminar round-jet flame. Only one computational cell is used in the circumferential direction, so that the simulation domain is essentially two dimensional. Since the simulations are performed in the framework of OpenFOAM, a wedge-shaped domain is required for OpenFOAM to treat the governing equations in cylindrical coordinates. y and x indicate the streamwise and radial coordinates. The dimensions of the boundaries have been set as multiples of , as indicated in Fig. 1. The domain length is set to be three times the flame height estimated from and , assuming a perfect conical flame. An equidistant grid is used for the simulations, which resolves the flame thickness with 20 cells with a grid resolution of m.
B. Boundary conditions
C. Numerical setups
The simulations have been carried out with an in-house solver developed in the framework of the open-source program OpenFOAM,14–16 which employs detailed calculations of the chemical reaction rates. The reaction mechanism developed by Li et al.17 is used for hydrogen/air combustion. The chemical source terms are treated via an operator splitting approach, where the governing equations for batch reactors are integrated in each computational cell over the computational fluid dynamics (CFD) simulation time step. Using the species concentrations in each cell at the beginning and at the end of the integration, a chemical source term that is averaged over the CFD time step can be computed and used as a source term in the governing equations. A more detailed description of this procedure can be found in Refs. 14–16 and 18.
The mixture-averaged transport model is used to consider differential diffusion, with transport coefficients of each species computed from the Chapman–Enskog solution of the kinetic gas theory. Soret and Dufour effects are not considered in this work. The latter one can generally be neglected as the contribution to the energy transport is negligible, even for hydrogen flames. As this paper focuses on very lean and rich hydrogen flames, the contribution of the Soret effect is less pronounced than at near-stoichiometric conditions. For example, the difference in the laminar flame speed with and without the Soret effect at is about 11%, while the difference at is 3% and at is 6%. The solver has been validated in previous works.15,18–29 The time steps are 0.5 s for the flame and 0.25 s for the flame, leading to a maximum convective CFL (Courant–Friedrichs–Lewy) number of 0.2 and a maximum Fourier number of 0.35.
A fully implicit scheme of second order accuracy for the time integration and a fourth order, unbounded interpolation scheme for discretization of the convective and diffusive terms have been used. The numerical solver employed in this work uses a pressure–velocity coupling via the PIMPLE (SIMPLE-PISO) algorithm. While this scheme is fully compressible in the sense of fully coupling density and pressure, it behaves more like a low Mach-number approach with regard to acoustics. Due to the implicit treatment and PIMPLE approach, the simulation time step is based on the convective CFL number, not the acoustic CFL number. Therefore, acoustic waves are generally not resolved since they are not the focus of this work. Instead, the non-reflective boundary at the outlet is used to prevent reflection in the case of nonphysical numerical oscillations. Additionally, a sponge region near the outlet avoids reflection of pressure waves back into the domain.
As the simulation is fully compressible, the oscillating mass flow at the inlet results in a fluctuation of the pressure inside the domain due to momentum conservation. The pressure fluctuation is particularly large in the case of high excitation frequencies. However, as the amplitude of the oscillating mass flow and the length of the computational domain are small in this work, the maximum pressure variations observed in the computational domain are on the order of 100 Pa. Therefore, the effect of pressure on chemical reaction rates is negligible.
IV. RESULTS AND DISCUSSION
A. Overall flame structure
Figure 2 presents contours of the calculated heat release rate and temperature T from the steady-state solutions (f = 0 Hz) for at the top and at the bottom. The solid lines in the temperature fields on the right denote the flame fronts given by the iso-surfaces of (see Table I). For the rich flame with has a maximum at the flame tip and a minimum at the flame base. The behavior is, however, reversed for the lean flame with , where the flame is extinguished at the tip, and is largest at the flame base due to the high diffusivity of hydrogen. The flame is negatively curved at the flame tip, which results in a defocused diffusion of hydrogen from the fresh to the burnt side. Therefore, the local equivalence ratio is smallest at the flame tip, which leads to the highest for the rich flame and the smallest for the lean flame, as shown in Fig. 3 by the contour-plots of , on the left for the and on the right for the flame. Due to the high diffusivity of hydrogen, is even smaller than that at the lower ignition limit for , so that, the flame extinguishes at the flame tip. On the contrary, the flame yields the highest positive curvature at the flame base, which results in the largest . The correlation of heat release or burning rate with regard to flame curvature results in a negative Markstein number for the current lean flame (Sl increases with stretch) and a positive Markstein number for the considered rich flames (Sl decreases with stretch).
Contours of heat release rate (left) and temperature (right) of steady-state H2/air flames for (top) and (bottom).
Contours of heat release rate (left) and temperature (right) of steady-state H2/air flames for (top) and (bottom).
Contours of calculated local equivalence ratio of steady-state H2/air flames for (left) and (right).
Contours of calculated local equivalence ratio of steady-state H2/air flames for (left) and (right).
Figure 4 provides an overview of instantaneous contours of calculated for the unsteady lean-premixed hydrogen–air flame at different f, where the flame shape deviates from that of the stationary flame shown in Fig. 2. From left to right, four snapshots of are used to represent one whole oscillation period for each Da, having a time interval of with . For the low frequency case with f = 125 Hz (see the first row with in Fig. 4), the flame is able to adapt promptly to the oscillating inflow and moves periodically as a whole with a large displacement from its time-mean position. As f approaches the natural frequency of the flame (see the second and third rows with Da = 2.5 and 1 in Fig. 4), the flame enters a pulsating mode and becomes wrinkled. With further increased f (see the lowest row with Da = 0.25 in Fig. 4), the flame turns into its steady-state solution without discernible motion. The same behavior with respect to flame response to unsteady flow has been found for the case with . Therefore, the overall flame response to an unsteady flow, in terms of changing its flame surface area through stretch, is attenuated with the increased fluctuation frequency of the flow (see also Fig. 7). This effect will be further elucidated in Sec. IV C.
Snapshots of calculated heat release rate of lean hydrogen/air flames at different excitation frequencies for .
Snapshots of calculated heat release rate of lean hydrogen/air flames at different excitation frequencies for .
B. Local flame dynamics and Markstein number
The flame particle (FP) method proposed in Refs. 13, 30, and 31 has been applied in this work to assess the effect of unsteady flow on the local flame dynamics. The method introduces a number of virtual, massless particles along the flame surface and tracks them during flame propagation. In this way, the time evolution of each individual flame surface element along with variations of flow stretch conditions and burning rates can be recorded, allowing a detailed evaluation of the time history effect of the flow on the flame's local dynamics. Figure 5 plots against Ka for the (top) and (bottom) hydrogen/air flames, which are calculated from the FPs using Eqs. (1) and (5). For the “on-the-fly” computation of Sl with Eq. (5), the line integration along the flame normal direction is performed locally over a distance of on both sides of the flame surface. The black points indicate the data obtained from the steady-state solutions, where increases with Ka for and vice versa for . Therefore, the Markstein number for the stationary flames is negative for the lean flame and positive for the rich flame considered in this work. The flame is stabilized at the burner orifice and yields a radius of approximately 3 mm (8 ). Therefore, cellular structures attributed to thermodiffusive instabilities cannot grow for the steady-state lean flame case with . For the oscillating lean flames, the corrugations on the flame surface are enhanced due to the thermodiffusive effect.
Correlations of local consumption speed with normalized stretch rate calculated from detailed numerical simulations of oscillating hydrogen–air flames at (top) and Φ = 4 (bottom).
Correlations of local consumption speed with normalized stretch rate calculated from detailed numerical simulations of oscillating hydrogen–air flames at (top) and Φ = 4 (bottom).
For the unsteady solutions plotted in Fig. 5, the FPs are seeded at the base of the flame at a time distance of with . These FPs, therefore, follow different spatial trajectories in time and experience different stretch conditions. This leads to a scattering for the distributions of vs Ka for the oscillating compared with the steady-state flames, as shown in Fig. 5. The diffusion–reaction balance prevailing the flame's local dynamics is disturbed by the unsteady velocity fluctuations. In this case, the flame requires a certain relaxation time for adjusting its internal structure and dynamics according to the transiently changing stretch, so that a previous state of is recorded by the FPs for the current Ka.
As shown in Fig. 5, the scattering of is particularly large when Da approaches unity ( ). In contrast, resembles the steady-state solutions for Da <1 and . The results indicate that the largest deviations of the local flame dynamics compared with that of the steady-state solution are due to the unsteady flow condition with ( ) or when the turnaround frequency of the flow approaches the flame's natural frequency. The internal structure of the flame sustained by the intrinsic reaction–diffusion balance yields the smallest resistance against external flow disturbances in this case, so that the flame is excited and falls into resonance with the unsteady flow fluctuations. Accordingly, the internal flame structure and the flame dynamics are altered considerably at , which leads to enhanced wrinkling of the flame (see also Fig. 4 for the case with Da = 1) and strong scattering of with respect to Ka.
Figure 6 plots the Markstein number calculated from linear fitting of all data pairs of vs Ka assigned to each excitation frequency shown in Fig. 5. As the linear Markstein correlation is generally valid only for small stretch values with Ka close to 0, the fitting procedure has been performed for the range of . In this way, corresponds to an unsteady or “turbulent” Ma averaged over different states of the flame during one oscillation period. The steady-state Markstein numbers Ma0 used in the ordinate axis of Fig. 6 for normalizing are for and for . is less than unity in general under unsteady flow conditions, which indicates a decreased sensitivity of the local flame dynamics to stretch. The smallest values of are found for Da close to unity, where the oscillation frequency matches with the intrinsic frequency of the flame. The flame is excited and wrinkled in this case, as shown in the third row of Fig. 4. Consequently, vs Ka reveals a scattering in Fig. 5. The correlation between and Ka is partly positive in the range of for due to formation of the enclosed envelope structures for and Ka, so that a significant decrease in at Da = 1 can be detected in Fig. 6 for this case.
Time-averaged Markstein numbers calculated from linear-fitting of data pairs of vs Ka shown in Fig. 5.
Time-averaged Markstein numbers calculated from linear-fitting of data pairs of vs Ka shown in Fig. 5.
C. Global flame response to unsteady stretch
Figure 7 depicts the temporal evolution of the global consumption speed St calculated from integration of the fuel reaction rate over the whole domain, see Eq. (6). In this way, St corresponds to the bulk flow rate of fuel prescribed at the inflow, assuming a balance of incoming and consumed fuel mass. The time axis is normalized with the oscillation period in order to compare St calculated from different excitation frequencies or Da. The dashed lines with a sinusoidal function in Fig. 7 show the oscillation curve of the bulk inflow velocity . St exhibits a phase delay with respect to the inflow, as indicated in the lower part of Fig. 7, which strictly increases with the increasing f or decreasing Da. In addition, St yields a damped oscillation amplitude compared with that of the inflow. In the cases with , the flame cannot respond to the highly unsteady flow and approaches its stationary solution, so that is decreased significantly and approaches 0. Figure 7 reveals that the ability of the flame to follow the flow fluctuation is strongly attenuated under the conditions of . In this case, the oscillation amplitude of St is significantly weakened, but the time periods remain as given by .
Time developments of global consumption speeds calculated from simulations of oscillating hydrogen/air flames at different forcing frequencies.
Time developments of global consumption speeds calculated from simulations of oscillating hydrogen/air flames at different forcing frequencies.
The same behavior is illustrated in the phase diagram in Fig. 8, where St (solid curves in Fig. 7) is plotted against (dashed curves in Fig. 7). The delayed response of global flame dynamics in terms of St to the unsteady inflow results in a hysteresis loop with a tilted elliptical shape. At the same time, the oscillation amplitude of St is reduced with the decreasing Da or increasing f, which is close to 0 for (quasi-horizontal curves in Fig. 8). The results reveal that the overall flame dynamics are sensitive to flow stretch only at low excitation frequencies, whereas its response to the flow is strongly attenuated in the case of high frequencies.
Global consumption speed St plotted against the inflow bulk velocity at different Da.
Global consumption speed St plotted against the inflow bulk velocity at different Da.
In summary, the flame wrinkling or enlargement of flame surface caused by turbulent flow is mainly attributed to unsteady fluctuations with Da >1. The flame regime “ideally stirred reactor” proposed in the classic Borghi–Peters regime diagram for Da <1 and , which is characterized by a considerably thickened flame brush and homogeneously mixed chemical state, is not visible.
D. Correlation of global consumption speed with total flame surface area
Figure 9 depicts as a function of the total flame surface area , where A0 is the area of the inlet boundary and the black dots represent the steady-state solutions. The solid line indicates the flamelet relation with , assuming that the local consumption speed Sl is equal to the unstretched laminar flame speed . Although not shown here, At oscillates in a similar way as St in time (as shown in Fig. 7) and no distinct phase delay between At and St has been observed. Therefore, and yield a quasi-linear correlation in Fig. 9 but deviate from equality. The result reveals that the increase in St is due to the increase in At and almost not affected by the unsteady nature of the flow.
Global consumption speed plotted against the total flame surface area of oscillatory H2/air premixed flames at different Da.
Global consumption speed plotted against the total flame surface area of oscillatory H2/air premixed flames at different Da.
The current round-jet flame experiences an overall high negative stretch due to curvature in the circumferential direction, which is compensated by the positive stretch component due to tangential straining, see Eq. (1). For the flame, over in Fig. 9 lies closely together with the reference curve (solid line) with , indicating that the global stretch Kat considering both curvature and straining of the whole flame surface is close to 0. In the case with , the bulk flow velocity at the inlet ( , see Sec. III and Table I) is much higher than that of the flame. Therefore, the rich flame yields overall a positive Kat due to the higher positive tangential strain compared with the negative stretch caused by flame curvature. This leads to I <1 according to Ma >0 for the case with Φ = 4. As a result, the correlation curves of vs in Fig. 9 lie below the reference line indicated with . Note that as the lean flame with is locally extinguished at its tip (see Figs. 2 and 4), the flame surface area evaluated from the iso-surface of is overestimated. As a consequence, the distributions of vs in Fig. 9 for the lean flame should be shifted toward smaller , which results in an increased I. As a conclusion, the averaged flame speed over the flame surface area is larger than in the case of Ma <0 and smaller than for flames with Ma >0, which is valid for both stationary and transient flames.
E. Efficiency of flame response to unsteady stretch
Efficiency factor considering the effect of flame response to unsteady stretch in dependence of Da.
Efficiency factor considering the effect of flame response to unsteady stretch in dependence of Da.
According to the results shown in Fig. 10, the flame can follow the flow fluctuations and change its flame surface area only under the conditions of Da >10 with ; the ability of the flame to respond to velocity fluctuations is strongly weakened for Da <1 with . This behavior for the effectiveness of global flame response to unsteady flow is illustrated in Fig. 11. The result suggests an extension of the regime diagram for flame–turbulence interaction with the asymptotic transition from a strongly responsive flame at Da >10 to a weakly responsive flame at Da <1 by a given turbulent Reynolds number. Its influence is particularly important in highly unsteady turbulent flows, where the energy cascade covers a wide range of time scales or Da. In this case, the flamelet correlation assuming a proportional increase in St with ( ) is strongly violated due to Da <1. In addition, even an excessively large stretch rate will not lead to flame quenching at Da <1 because the flame is not able to sense the highly unsteady flow stretch.
Effectiveness of flame–turbulence interactions considering the influence of flow unsteadiness in terms of Da.
Effectiveness of flame–turbulence interactions considering the influence of flow unsteadiness in terms of Da.
F. Discussion
The current work has reproduced the existing findings from the previous studies by Im and Chen9 applied to 1D oscillating counter-flow flames, where the flame speed has been shown to be less sensitive to stretch under unsteady flow conditions. In the current work, a round-jet flame setup is used, which leads to a substantially extended stretch range, considering both effects caused by tangential strain and curvature, and covering both negative and positive stretch. Furthermore, the global flame response in terms of the overall flame surface area is used here to evaluate the global flame dynamics, whereas the concept of flame surface is less meaningful for the planar flame setup. Therefore, the current work extends considerably the available knowledge based on 1D studies.
The obtained results can be incorporated into combustion modeling concepts in order to account for the time history effect on the response of global flame dynamics. For that purpose, a fitting function may be derived to model the correlation of E with Da shown in Fig. 10, which can be used to evaluate a cell-averaged burning rate based on St, for instance, in the framework of a RANS (Reynolds-Averaged Navier–Stokes) or LES (Large Eddy Simulation) concept.
Note that Da yields a non-uniform distribution along the flame surface in a general turbulent flame: on one hand, the broadband fluctuations of the flow cause variation of the local turbulent time scale, and on the other hand, the flame thickness and the flame speed can be changed locally by the turbulent flow, which results in a different chemical time scale. Furthermore, the flame–turbulence interactions take place over a wide range of time scales in a real turbulent flame, whereas a fixed forcing frequency is applied to interact with the flame in the current study.
V. CONCLUSIONS
Detailed numerical simulations have been performed for a series of oscillating premixed H2/air flames with pre-defined frequencies f. The objective is to assess systematically the impact of flow unsteadiness on the local and global flame dynamics. The flame requires a relaxation time to respond to the oscillating flow so that flame speed Sl is less responsive to stretch under unsteady conditions. The statement is valid for both lean and rich premixed flames with oblique flow and strained and curved flame structure. As a result, the magnitude of the Markstein number decreases in an unsteady flow. The flame is excited with an enhanced wrinkling, and is smallest when f approaches the natural frequency of the flame, i.e., .
The global consumption speed St, which correlates almost linearly with the total flame surface area, exhibits a phase delay of and a reduced amplitude of with respect to the oscillating inflow. increases with f, whereas the reversed trend is found for . The damped oscillation of St due to flow unsteadiness has been modeled with an efficiency factor E in dependence of Da, which reveals the asymptotic limits with for and for .
In conclusion, the time history effect or unsteadiness of the flow plays an essential role for the correct description of local and global flame dynamics regarding flame–flow interaction. The generally used flamelet assumption is valid only for large Da, having E >0.8 for Da >10 in the current work. This work has focused on premixed hydrogen/air flames excited at fixed single frequencies. Further studies for additional fuel/air mixtures and considering the interactions between different flow time scales are required to justify the general validity of the obtained results.
ACKNOWLEDGMENTS
The authors acknowledge the financial support of the Helmholtz Association of German Research Centers (HGF), within the research field MTET (Materials and Technologies for the Energy Transition), subtopic “Anthropogenic Carbon Cycle” (38.05.01). This work utilizes computing resources from the supercomputers Hawk at HLRS Stuttgart and HoreKa at SCC/KIT. The work leading to this publication was supported by the PRIME programme of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Feichi Zhang: Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review and editing (lead). Thorsten Zirwes: Investigation (lead); Methodology (lead); Validation (lead); Writing – review and editing (supporting). Yiqing Wang: Investigation (equal); Methodology (equal); Software (equal). Zheng Chen: Conceptualization (equal); Investigation (equal); Methodology (equal). Henning Bockhorn: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – review and editing (supporting). Dimosthenis Trimis: Funding acquisition (equal); Project administration (equal); Supervision (equal). Dieter Stapf: Funding acquisition (equal); Project administration (equal); Supervision (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.