As a common water conveyance structure, aqueduct is widely used in cross-regional water conveyance and diversion projects. However, large double-trough aqueduct structures are prone to damage under earthquake, resulting in water leakage and even interruption, endangering the safety of people’s lives and property in surrounding areas. In this paper, the nonlinear reinforced-concrete material program is compiled based on the FORTRAN language. Taking the Shuangji River aqueduct project as an example, the numerical model of the fiber beam element of the double-trough aqueduct structure is established, and the nonlinear dynamic time-history response analysis of different ground motion inputs is carried out. Based on the midspan displacement and acceleration time history, pier bottom bending moment and shear force time history, the nonlinear response law of large aqueduct structure is studied. A criterion for determining the overall collapse of aqueduct structure based on effective energy is proposed. The correctness and feasibility of the nonlinear numerical model of the aqueduct structure and the dynamic catastrophe analysis method are verified by focusing on the overall seismic response of the aqueduct structure and the material yield status of key parts, and comparing with the collapse criterion based on the limit value of pier top displacement angle. The conclusions will provide a basis for seismic analysis, optimization design, and operation and maintenance safety assessment of large aqueduct structures and have important theoretical significance and engineering application value.

1. Introduction

Large aqueduct structures, as important structures for transporting water across complex terrain, have been widely used in the South-to-North Water Transfer Project. Many large reinforced-concrete aqueduct structures have been or will be built in high-seismic-intensity areas, and the seismic fortification intensity of some areas is as high as 9°. Under the action of strong earthquakes, large reinforced-concrete aqueduct structures are prone to damage and destruction. In addition, various problems occurred, such as the short-time concentrated discharge of water from the water transmission line, difficulty in postearthquake repair, and local secondary disasters [1, 2]. Therefore, the study of the nonlinear responses of large double-channel aqueduct structures under earthquakes is urgently required.

At present, the seismic responses of reinforced-concrete aqueduct structures are mainly studied through experiments or numerical simulations using finite element software. Li et al. [3] used the finite element software ANSYS, by establishing a numerical model considering the interaction between fluid and structure, and the time history and response spectrum of the seismic effect of aqueduct structure are calculated. Several authors studied the effects of soil-aqueduct-water coupling via numerical analysis of the transverse channels of aqueduct structures and performed pseudo-dynamic tests of structure models. The results showed that the water conditions can play a role in shock absorption [46]. Wang et al. [7, 8] established a dynamic analysis model of beam segment elements considering the impact effect at the expansion joint of an aqueduct, and they discussed the impact of the collision at the expansion joint of adjacent aqueducts on the seismic responses of aqueduct structures under the seismic action of a longitudinal channel. Wang et al. [9] studied the dynamic interactions between an aqueduct and a water body through scale tests using a shaking table model. The above studies explored the seismic performance of the aqueduct structure and its influencing factors from different aspects and provided an important reference for the seismic design and evaluation of the aqueduct structure. However, most studies in numerical simulations used linear elastic models, lacked consideration of the nonlinearity of reinforced-concrete materials, and there is still a big gap between the theoretical analysis results obtained and the actual situation, and solid elements common numerical simulation calculation time and high cost. This paper used ABAQUS to compile a reinforced-concrete material subprogram suitable for aqueduct structures and established a numerical model of a large-scale double-trough aqueduct structure with refined fiber beam elements, and the nonlinear response analysis method of aqueduct structure under earthquake is studied.

In recent years, some researchers have also carried out fruitful research on the seismic reliability analysis of aqueduct structures. Zhang et al. [10, 11] studied the reliability of the water conveyance function of the aqueduct structure through the physical response of the structure under the action of an earthquake. Ma et al. [12] established the reliability analysis model of the aqueduct system based on the main failure mode and the comprehensive correlation coefficient method and discussed the system reliability of the multiside walls of the piles and beams of the Minghe Aqueduct. However, the research on the nonlinear analysis of aqueduct structure and the judgment criterion of the overall structure failure under the action of the strong earthquake is relatively insufficient. Key problems such as strain localization and size effect in seismic damage analysis of traditional aqueduct structures have not been effectively resolved, which hinders the development of current seismic research on aqueduct structures. In fact, the nonlinear response analysis of the aqueduct structure often refers to the research results of the seismic response of the bridge structure. Parool and Rai [13] carried out a study on the seismic performance of a multispan steel-supported simple-supported girder bridge. The seismic response analysis of a four-span bridge model showed that the failure of the main-span steel support was the main cause of bridge failure; Li et al. [14] studied the dynamic response of continuous girder bridge structures under strong earthquakes, designed and produced a two-span continuous girder bridge scale model, considered different types of supports, and selected four types of seismic waves reflecting the soil characteristics of different sites and the seismic simulation shaking table array test of the scaled model bridge under different ground motion input test. At present, since a unified collapse criterion has not yet been formed, a simplified collapse criterion is often used from the perspective of macroscopic physical quantities, and it is difficult to determine the critical limit state of the overall structure collapse by numerical or experimental analysis for a certain component. Moreover, due to its “top-heavy” structure, the aqueduct structure has many special features compared with ordinary bridge structures in terms of working mode. The collapse process under the action of catastrophic earthquakes is often extremely complicated. Therefore, this paper based on the established refined numerical model considering material nonlinearity deeply discusses the failure mechanism of large-scale double-trough aqueduct structure under catastrophic earthquake. From the perspective of effective characteristic energy and effective input energy, a material-unit-structure multiscale aqueduct structure collapse criterion is proposed, and a refined analysis method for the whole process of large-scale aqueduct structure dynamic catastrophe is established.

2. Fiber Beam Theory

Three-dimensional solid element models and fiber beam element models are mainly used in the numerical analysis of reinforced-concrete beam and column members. In recent years, it has been found that fiber beam elements can reflect the typical nonlinear mechanical behaviors of reinforced-concrete structures with small numbers of elements [15]. The fiber beam element model is shown in Figure 1; that is, the constructed cross section is divided into several fiber sections along the length of the rod. Based on the assumption of a flat section, the relationship between the force and deformation of the entire section of the member can be calculated according to the uniaxial stress-strain relationship (see Section 4 for the cyclic stress-strain relations) of the material corresponding to each fiber [16].

The model is based on the following assumptions: (1) the deformation of the fiber section meets the assumption of the plane section, and the strain of all fibers in the section can be obtained only by determining the strain and curvature in the middle of the section. (2) The slip between reinforcement fiber and concrete fiber is ignored, and the influence of section shear and torsional deformation is ignored. (3) Each fiber material satisfies the uniaxial stress-strain constitutive relationship.

In this model, the element is divided into several integral segments, and the cross-sectional force at the integral point [17] is obtained through a linear interpolation method. The linear interpolation function is shown as follows:

The flexibility matrix of the cross section is integrated along the axial direction of the member to obtain the element flexibility matrix, as follows:

The unbalanced force vector of the cross section is transformed into the residual deformation, and the deformation increment of the element in the next iteration is obtained by integrating along the axial direction of the member, as follows:

The section stiffness matrix is shown as follows:

The stiffness matrix of the thin-walled fiber beam elements of the aqueduct thin-walled structure is fully considered in the model. The coupled vibrations of the lateral bending and the torsion and warpage deformation of the aqueduct thin-walled structure are fully considered. The deformation of the aqueduct section satisfies the assumption of a flat section. The reinforcing steel bar and concrete fibers are in a uniaxial stress state. The elastic and plastic tensile constitutive model of concrete and the steel hardening constitutive model under repeated loading are adopted. Because the buckling deformation effect on the overall deformation of the aqueduct thin-wall structure is larger, the deformation at any point within the unit is obtained by introducing the warping function of the displacement, as follows:where , , represent the linear displacement function in three directions; , , are the angular displacement function; , , are the angular velocity.

When the aqueduct structure begins to undergo elastoplastic deformation and large-displacement deformation, the influence of the material nonlinearity and geometric nonlinearity on the stiffness matrix is considered. The strain matrix is as follows:where and are the total strain matrix and linear microdisplacement strain matrix, and is a geometrically nonlinear strain matrix. The stiffness matrix of the fiber beam segment element model of the aqueduct thin-walled structure considering geometric nonlinear deformation is as follows:

3. Determination Criterion of Structural Dynamic Collapse Based on Energy

Many studies have shown that the failure modes and collapse mechanisms of structures in the collapse process are extremely complex [18], and most of them are caused by the comprehensive action of many factors. It is difficult to accurately evaluate the nonlinear collapse behaviors of structures under catastrophic earthquake loads by using a single macrostructural response parameter.

The response of the structure under an earthquake is actually a process of energy input, transformation, and consumption. For a general multidegree-of-freedom dynamic system, the equation of motion is as follows:where is the mass matrix of the system, is the damping matrix of the system, is the restoring force vector, is the external excitation vector, , , and , respectively, are the acceleration, velocity, and displacement vectors of the nodes.

The energy balance equation is obtained by multiplying both sides at the same time and integrating 0 to t time:

Among them, represents the kinetic energy of the system at time t; represents cumulative plastic energy consumption during the period; represents the deformation energy of the system; represents the total energy input from outside during the period in the nonlinear analysis of structures; can be obtained by integrating the strain energy of the unit body in the region of the structure:

Among them, is stress tensor, is strain tensor, and represent all domains to be solved.

The total strain at time t is decomposed into elastic part and plastic part , and thenwhere represents the elastic deformation energy of the system, and represents the plastic deformation energy of the system.

A dynamic critical collapse state criterion based on the effective intrinsic energy and the effective input energy of the structure was proposed in the literature [19, 20].

The effective intrinsic energy is defined as follows:

The effective input energy is defined as follows:where is the kinetic energy of the system, and is the elastic deformation energy of the system. The criterion for the overall collapse of the structural system is as follows:

Therefore, the function that determines the critical state of the structural collapse can be defined as follows:when for all , the structure maintains dynamic stability. The structure collapses when .

4. Material Constitutive Relationship

4.1. Rebar Reinforcement

The reinforcement constitutive model proposed in the literature [21, 22] was adopted in this paper (as shown in Figure2). The steel skeleton line is composed of two sections: OA (elastic section) and AB (strengthening section). In the OA section, the elastic modulus is , and the yield of the steel bar is ignored. When entering the AB section, the stiffness coefficient is 0.01. Unloading begins when the maximum strain point is reached.

4.2. Concrete Constitutive Relationship

To better simulate the material failure behavior of a large aqueduct structure under strong earthquakes, a concrete constitutive model [23, 24] considering the tensile strength was selected in this paper (as shown in Figure 3). The initial modulus of elasticity is the slope of the tangent line of the skeleton line at the origin (). The tensile skeleton line is a double line with a softening section, and the elastic modulus of the rising section is . After reaching the axial tensile strength , the line enters a descending section with a stiffness of . The compression skeleton line is determined by the following formulas:where is the axial compressive strength of concrete, and is the concrete strength corresponding to the ultimate compressive strain .

The degradation of the unloading and reloading stiffness is controlled by a virtual point R, whose position is determined by the damage parameter . The corresponding strain at point R is calculated according to (18), and the corresponding stress is . The value range of is shown in (19).

After R point is determined, the loading and unloading rules within the compression skeleton line are as follows: the intersection point of the line between the unloading point and the R point with the horizontal axis is the residual strain , and the slope of the line is the damage stiffness, denoted as , where is the compression damage coefficient. When unloading, it is first unloaded according to the initial stiffness , and when the unloading reaches the residual strain point and the stiffness is , the unloading stiffness was changed to . During reloading midway through unloading, the initial stiffness was added load until the connection between the unloading point and R point is reached, then load according to the damage stiffness (as shown in Figure 2).

According to the paper [23], the constraint index is adopted to consider the constraint effect of stirrup on concrete in reinforced-concrete beam-column members, and the concrete stress-strain curve equation in the code is modified. The formula for calculating is as follows:where is the yield strength of stirrup; is the volume ratio of stirrup in concrete beam-column members.

In this case, the compressive strength and the peak strain are as follows:

5. Engineering Examples

In the following examples, the upper structure was a prestressed concrete rectangular groove structure, and the body was C50 concrete. The elastic modulus was 3.45 × 1010 Pa, the bottom width of the cross section was 5.00 m, the depth of the groove was 4.25 m, the total length was 197.80 m, and the dimension of the cross section through which water passed was 4.00 m × 3.40 m × 2 (the number of wide and high slots). The lower part was the foundation of the bored pile supported by a single row frame. The concrete density was 2500 kg/m3, the elastic modulus of the material was designed as 2.55 × 1010 Pa, and Poisson’s ratio was 0.3. The support for a single row of the three-column structure had a height of 10.40 m and a circular cross section with a diameter of 800 mm. The cover beam and the connection beam were arranged between the two columns, both of which were rectangular sections. The dimensions of the upper beam section were 1100 mm × 1000 mm. The elevation structure diagram and cross section diagram of the aqueduct are shown in Figures 4 and 5, respectively.

The fiber beam element in ABAQUS was used to establish a numerical model of the above-mentioned double-channel aqueduct structure, as shown in Figure 6. The cross sections of the elements constructed in the finite element model were consistent with an actual engineering design. The B31 [25] beam element was selected for the erection and cover beams, as shown in Figures 7 and 8 [26]; the open thin-walled beam element B31OS [25] was selected for the aqueduct body, as shown in Figure 9; and rebar was added to the concrete through common node technology based on the concept of equivalent area and equivalent position [25]. The rubber bearing was replaced by a spring element, and the horizontal stiffness was set as K = 1.53 × 106 N·m. In the numerical model, there were 121 nodes and 110 units in total. Using the Westergaard model, the hydrodynamic pressure in the trough is applied to the aqueduct wall in the additional water mass method.

To explore the nonlinear dynamic responses of large aqueduct structures under strong earthquakes, seismic waves of four different site types were input [27]: Northridge wave, El Centro wave, Taft wave, and Tianjin wave. They were input along the transverse channel with the maximum amplitude modulated to 0.6 g. The damping was Rayleigh damping with a damping ratio ξ = 0.05.

6. Result Analysis

6.1. Frequency and Mode Shape of Natural Vibrations of Aqueduct

It can be seen from Table 1 that the first, second, and fourth orders of the finite element model of the double-slot aqueduct structure are lateral vibration; the third, fifth and seventh orders are lengthways vibration; the sixth and tenth orders are vertical vibration; and the eighth and ninth orders are pier bending vibration. This showed that the lateral stiffness of the double-channel aqueduct structure was relatively weak, and the aqueduct shelf was prone to bending failure under the action of a strong earthquake. Therefore, the dynamic nonlinear analysis of the double-channel aqueduct structure under the action of a lateral earthquake will be carried out.

6.2. Stress and Strain of the Fiber at the Column Bottom Section of the Aqueduct

Since the bottom of the column of a large aqueduct structure could easily enter the elastoplastic deformation regime, the section with the maximum plastic strain of the concrete in the column member was extracted for analysis. Figures 10 and 11 show the stress-strain curves of the concrete fiber at the bottom of the shelf. The stress-strain curves show that the compressive strength of the concrete element was much greater than the tensile strength. When the initial seismic excitation was small, the structure was in the linear elastic stage, and the concrete unit entered the plastic stage with the increase in the seismic wave peak value. Under the action of the El Centro wave, the positive tensile strength of the material was nonlinear, but at this time, the seismic excitation was small, and there was no residual deformation after unloading. Under the action of the El Centro wave, the concrete unit was subject to tensile yield. At this time, the tensile stress of the section was reduced to zero, but the strain continued to increase until the concrete crack exited the work. When the structure was subjected to a reverse force, the crack surface was closed, and the crack surface effect was generated. The tensile strength of concrete could be recovered when the structure was again under positive tension. The hysteresis curve is disordered because of the high axial compression ratio of aqueduct structure, and the irregular elastic-plastic deformation of reinforced-concrete members occurs when the local seismic excitation increases.

6.3. Nonlinear Response Analysis of Aqueduct Structure

To compare the seismic response differences between the nonlinear state and the linear elastic state of the aqueduct, two calculation models were used: ① the elastic model, in which the structural materials always retained linear elasticity during the earthquake, and ② the elastoplastic model, in which the material exhibited nonlinear characteristics during the earthquake.

Figures 1215 show the time-history curves of the aqueduct structure under elastic and elastoplastic conditions due to lateral ground motion. The displacement, velocity, and acceleration under elastoplastic conditions were greater than those predicted by the elastic model. The peak acceleration response increased by 75%, the peak velocity response increased by 50%, and the peak displacement response increased by 68%. The reason was mainly that after entering the elastoplastic phase, the stiffness matrix changed, and the influence of the concrete after reaching the yield stiffness decreased. The amplitude of the seismic wave was small in the initial stage of the earthquake, the aqueduct structure was in the elastic stage, and the seismic responses under the elastic and elastoplastic conditions were basically consistent. With the increase in the seismic amplitude, the tensile yield of the structure entered the nonlinear phase. The acceleration and displacement responses of the structure were significantly greater in the nonlinear phase than those under linear elastic conditions.

Table 2 compares the values of the shear force, bending moment, acceleration, velocity, and displacement of the cross section at different positions under the action of the El Centro wave under elastic and elastoplastic conditions. The results showed that the peak responses of the shear force and bending moment under elastic conditions were larger than those under elastoplastic conditions. The reason was that the elements in the total stiffness matrix of the structure decreased significantly after entering the nonlinear state, and plastic damage occurred due to the large plastic deformation of the members, which led to a corresponding increase in the displacement and acceleration. In addition, the plastic damage increased the damping of the structure and the plastic energy dissipation, which reduced the decrease in the internal force of the structure caused by the earthquake to a certain extent. By comparison, it was found that the peak response of the midspan section was the largest.

In order to compare and study the influence of seismic waves of different site types on the nonlinear seismic response of large double-slot aqueduct structure, four seismic waves suitable for different site types are selected (as shown in Figure 16). Northridge wave of weak site type, Taft wave of medium-soft site type, EL Centro wave of medium-hard site type, Tianjin wave of hard site type, input along the transverse trough, the maximum amplitude are adjusted to 0.6 g. At the same time, using Rayleigh damping, damping ratio ξ = 0.05.

Figures 1720 show that the seismic responses of the aqueduct structures were different under the excitation of different seismic waves at different sites. The peak displacement responses of the Tianjin, Northridge, Taft, and El Centro waves were 0.23, 0.145, 0.078, and 0.11 m. The displacement responses of the aqueduct structure at the soft site were larger, and the maximum difference between the displacement responses at this and other sites was 0.149 m. Therefore, the influence of the site type on the nonlinear response of the aqueduct structure should be considered in the seismic design. Figure 21 compared the displacement response of the elastic model and the elastic-plastic model of the aqueduct structure under the ground motion excitation of four different site types. It can be clearly seen that the displacement response of the elastic-plastic model of the aqueduct structure under different seismic waves is greater than that of the elastic model. The aqueduct structure can dissipate the seismic energy in the form of plastic energy dissipation under the earthquake action to reduce the damage and failure of the structure under the earthquake.

7. Determination of Dynamic Collapse of Aqueduct Structure Based on Effective Energy

Based on the numerical model of the fiber beam element of the aqueduct structure, four Northridge waves with different amplitude sizes were input to simulate the collapse of the structure. The seismic amplitude of each working condition is shown in Table 3.

An aqueduct structure will enter the strong nonlinear stage under a disastrous earthquake load. At this time, the collapse of the structure cannot be accurately predicted by physical responses, such as the displacement, velocity, and acceleration. Figure 22 shows the characteristics of the effective energy of the aqueduct structure and the effective energy input for comparison. The effective energy values of the structure under conditions 1 and 2 were always less than the input energy. Under condition 3, the effective energy exceeded the input energy for the first time at 5.7 s, whereas this occurred at 5.3 s under condition 4. According to the effective energy collapse criterion, it was determined that the structure under conditions 1 and 2 maintained dynamic stability, while the aqueduct structure under conditions 3 and 4 collapsed.

Figure 23 shows the time-history curves of the energy of the structural system under conditions 1 (dynamic stability) and 3 (dynamic instability). Among them, ALLKE represents the kinetic energy in the system, ALLPD represents the plastic energy in the system, ALLSE represents the elastic energy in the system, and ALLWK represents the energy exerted by the external force on the entire model system. The structure was basically in the linear elastic stage in the initial stage of the earthquake, and the energy in the system also reached the maximum value at the moment of the earthquake’s peak action. Under condition 3, when the aqueduct structure was on the verge of collapse, the energy of the system suddenly changed. Under condition 1, the plastic energy dissipation of the structure continued to increase, with a maximum value of 192 kJ. Under condition 3, the plastic energy dissipation of the structure rapidly increased to 413 kJ within 4.6 s. At this time, the structure experienced dynamic instability, while the plastic energy dissipation basically remained unchanged. The plastic energy dissipation played an important role before the structural dynamic instability occurred.

In order to verify the reliability of the collapse criterion of the aqueduct structure based on effective energy, Figure 24 shows the time-history curves of the displacement angle at the pier top of the aqueduct under four working conditions. It can be seen from Figure 24 that the displacement of the pier top of the aqueduct structure does not appear drift phenomenon in condition 1 and condition 2, which reflects that the aqueduct structure does not collapse under these two conditions. However, it is found from the figure that the time history of the displacement angle at the pier top of the aqueduct in Condition 3 and Condition 4 shows a large drift, indicating that the aqueduct structure is unstable and collapsed due to structural stiffness degradation and material softening. Referring to the previous research results, the initial collapse of the aqueduct structure is roughly defined with the displacement angle of the pier top exceeding 5% for the first time. It can be seen from Figure 24 that the initial collapse time of the aqueduct structure in Condition 3 is about 5.8 s, and that in Condition 4 is about 5.3 s. By comparison, it is found that the collapse criteria of the aqueduct structure based on the effective energy criterion and the collapse criterion of the aqueduct structure based on the displacement angle of pier top can accurately determine whether the aqueduct structure collapses or not. When the limit value of pier top displacement angle is defined as 5%, the prediction of the initial collapse time of the two methods is basically consistent, which verifies the reliability of the collapse determination method of aqueduct structure based on an effective energy criterion. The above results show that the collapse criterion based on effective energy can accurately predict the dynamic stability of the aqueduct structure under an earthquake, whether it is under collapse or not.

The material stress-strain curve can further explain the collapse behavior of the aqueduct structure from the microscopic point of view of material yield. Taking condition 3 as an example, the stress-strain curves of concrete elements at the bottom and top of the left pier and middle pier are extracted, respectively, as shown in Figure 25. It can be seen that the increase in plastic energy dissipation caused by plastic damage of a large number of concrete members in condition 3 is basically consistent with the change in plastic energy dissipation in Figure 25. Compared with the stress-strain curves of the four parts, the pier in the aqueduct structure does not yield in the early stage of the earthquake, and the top of the pier is mainly a compressive failure, and the bottom of the pier is more prone to tensile failure in the early stage of earthquake. It can be seen that the concrete elements at the bottom of the left pier, the bottom of the middle pier, and the top of the middle pier are subjected to tensile yield under the action of the disastrous earthquake, and the compression does not reach the yield state. The concrete element at the top of the left pier reaches the yield state under tension and compression, and there is a certain degree of yield failure in the four parts.

8. Conclusion

In this paper, the nonlinear dynamic analysis model of double-slot aqueduct structure based on fiber beam element is established by taking the nonlinear development process analysis of large aqueduct structure as the breakthrough point. Taking Shuangji River aqueduct in the middle route of South-to-North Water Diversion Project as an example, the dynamic nonlinear whole-process response analysis under different working conditions is carried out, and the collapse criterion of the aqueduct structure based on effective energy is proposed. The main work and research conclusions of this paper are summarized as follows:(1)The seismic response time history calculation results of the linear elastic and nonlinear model of the aqueduct structure are compared and analyzed. The nonlinear response value of acceleration, velocity, and displacement of the aqueduct structure under seismic action is greater than that of linear elastic dynamic response value, and the internal force response value is less than that of the linear elastic response value.(2)According to the four site types of hard site, medium-hard site, medium-soft site and soft site, four representative seismic waves are selected to study the nonlinear response of the double-slot aqueduct structure under different sites and different amplitude-modulated ground motions. The nonlinear response value of the aqueduct structure increases with the weakening of the site conditions.(3)The criterion of overall instability and collapse of aqueduct structure based on effective energy is proposed, and compared with the maximum displacement angle criterion. The results show that when the effective characteristic energy formed by the aqueduct structure under earthquake action is always less than the effective input energy, the aqueduct structure will maintain dynamic stability; when the effective characteristic energy in the aqueduct structure system exceeds the effective input energy for the first time, the aqueduct structure will collapse at this time.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Jianguo Xu and Jinpeng Zhang contributed to conceptualization; Jinpeng Zhang and Chong Wu contributed to data curation; Chong Wu and Jianguo Xu contributed to formal analysis; Yupeng Geng and Jinpeng Zhang contributed to funding acquisition; Chong Wu contributed to methodology; Chong Wu contributed to original draft preparation; Chunyu Zhang and Jianguo Xu contributed to review and editing; Yupeng Geng contributed to project management.


This research was supported by the National Natural Science Foundation of China (No. 52079128) and the Science and Technology Project of Henan Province (Nos. 182102210013 and 212102310289). The authors thank LetPub (https://www.letpub.com/) for its linguistic assistance during the preparation of this manuscript.