Abstract

To study the shear force distribution laws of a box girder with a single-box multichamber (SB-MC) configuration for different supporting conditions, numbers of webs, stiffness of end diaphragm, and web thickness values, a box girder with SB-MC was numerically simulated using three-dimensional finite element model. According to the comparison results of web shear force, the concept of η, a shear-increased coefficient for webs, was introduced. The results show that supporting conditions and chambers have a significant impact on the shear-increased coefficient η, and end diaphragm must be set up in the 3D finite element model when calculating η. Nonlinear analysis shows that in the elastic phase, the shear-increased coefficient η basically does not change, but in the cracking stage, the coefficient η of each web changes with the degree of web cracking, and side-webs (S-Webs) reach the ultimate load first. The variation of the web thickness hardly affects the distribution of the shear force, so the method to adjust the web thickness of S-Web was proposed according to the result of shear-increased coefficient η to improve the shear resistance of the box girder.

1. Introduction

Box girder is an important type of bridge structure. For box girder, many scholars have conducted research on shear lag, effective width [15], and live load distribution factors [68]. There are multiple webs in the cross section of a single-box and multichamber (SB-MC) box girder, and their spatial mechanical behavior is complicated. The laws of their shear distribution have important significance and value to effectively ensure the excellent work of the web.

To study the shear resistance, failure modes, and internal force redistribution of box-shaped bridges, the finite element model of the box girder was established and numerically analyzed [9]. This study revealed the effects of factors such as the shear-span ratio, concrete strength, stirrup ratio, and web thickness on the failure modes and shear capacity of the box girder.

The shear resistance of the box girder was experimentally analyzed [10]. The results show that it is unsafe to use the equivalent I-beam to calculate the shear capacity of the box girder. Then, methods based on the modified pressure field theory were compared with the experimental results, including the double-section method, single-section method, simplified analysis method, and finite element method. The conclusion is that the analysis methods based on the modified pressure field theory could better predict the shear capacity of box beams.

A comparative analysis of the shear design codes of various countries was carried out, and the influence of various influencing factors on the calculation of shear capacity of different national codes was analyzed [11].

There were few results on the shear distribution in the webs of the box girder with SB-MC, so the regular processing method [1214] is each web equally shares the shear force in the shear calculations of the box girder. However, a preliminary study showed that the shear force is not evenly distributed in webs because of factors such as the supporting conditions, number of webs, stiffness of end diaphragm, and web thickness.

Since the finite element software performed well in civil engineering simulations, many practical projects were calculated and analyzed using this method, and the results proved to be positive [1517]. This study uses the MIDAS/FEA software to calculate and analyze the models of the SB-MC box girder under various conditions. Based on the summary of actual sharing of shear forces of each web, the laws of shear distribution of the box girder are obtained, and related design suggestions are proposed.

2. Research on the Shear Distribution of the Box Girder

2.1. Extraction Method of the Shear Force in Webs with the Three-Dimensional Finite Element Model

The stress analysis of a box girder is relatively complicated. With the conventional theory for the box girder, only the analysis results under specific conditions of a typical section can be obtained. It is difficult to obtain the effects of factors such as the change in internal force of the full bridge section and the difficulty in reflecting different supporting conditions. The main goal of this party is to use the software Midas/FEA to build a linear elastic three-dimensional finite element models of the single-box multichamber (SB-MC) box girder under the action of self-weight of the bridge and second-stage dead load in different cases, such as different supporting conditions, number of webs, stiffness of end diaphragm, and web thickness (Figure 1).

Using the function sum of internal forces in specific location in the software to extract the shear forces in each web is varied with the support conditions, numbers of webs, stiffness of end diaphragm, and web thicknesses. Then, the shear distribution of the web could be obtained.

The three-dimensional finite element model can well reflect the mechanical behavior of the structure, but in general, the postprocessing of the three-dimensional finite element analysis only provides the displacement, stress, and strain. The internal forces such as bending moment, axial force, and shear force are difficult to obtain. To study the distribution laws of the shear force in the SB-MC box girder, it is necessary to extract the shear values of each web of the box girder model in a specific section using the function sum of internal forces in specific location.

Taking the box girder with a single box and three chambers as an example, the sum of the internal forces in the specified section of the box girder was extracted. The main process is as follows: first, in the postprocessing after the model analysis, the objective section of the box girder was obtained by a flat surface, which is established by three points and intersects the box girder (Figure 2(a)).

Then, the objective section was selected by the function sum of internal forces in specific location. Finally, the sum of the internal forces in the selected section can be outputted by clicking Text. The exported internal forces mainly contain Fx (axial force in the X direction), Fy (shear force in the Y direction), Fz (shear force in the Z direction), Mx (bending moment around the X-axis), My (bending moment around the Y-axis), and Mz (bending moment around the Z-axis), as shown in Figure 3.

The research on shear distribution of webs need to extract the shear force of each web and section. Before the sum of the internal forces Vi of web, only the area of the objective web that requires the extraction of shear force was activated. Then, in the sum of internal forces, the result is the internal forces of the objective web, i.e., the shear force of the objective web can be accessed, as shown in Figure 2(b). If the sum of the internal forces Vt of the section needs to be extracted, the entire section of the girder needs to be activated (Figure 2(a)).

2.2. Models to Analyze the Shear Distribution of Webs

This study mainly extracts the shear values shared by all webs of the SB-MC box girder under different conditions of support cases, numbers of webs, and web thicknesses; then, we summarize the shear distribution laws of the web.

The basic parameters of the box girder are span: 25 m; beam height: 1.60 m; roof thickness: 20 cm; and floor thickness: 20 cm. The middle segment of the girder is 10 m long, the segment with widening webs is 2.5 m long, the segment with widened webs is 5.0 (4.6, 4.4, 4.0, 3.8, and 3.5) m long, and end diaphragm is 0.0 (0.4, 0.8, 1.0, 1.2, and 1.5) m long. The variables studied are the support conditions, web thickness, the stiffness of end diaphragm, and number of webs. The analysis model corresponding to each variable is described as follows.

2.2.1. Models with Different Support Cases

The cross section of the box girder has a single box and five chambers, whose roof width is 20.4 m, floor width is 15.4 m, the thickness of end diaphragm is 1.2 m, and center distance between two adjacent chambers is 300 cm. The thickness of the web linearly changes from 40 cm to 60 cm. Models with different support conditions are introduced as follows:(1)Model SPT-2 used two bearings, whose center distance is 14.2 m, as shown in Figure 4(a)(2)Model SPT-3 used three bearings, whose center distance is 7.1 m, as shown in Figure 4(b)(3)Model SPT-4 used four bearings, whose center distance is 4.7 m, as shown in Figure 4(c)

2.2.2. Models with Different Numbers of Chambers

The cross section of the box girder has a single box and five chambers, whose roof width is 20.4 m, and the floor width is 15.4 m. The thickness of the web linearly changes from 40 cm to 60 cm, and double bearings were used for this model with a center distance of 14.2 m.

Four models were comparatively analyzed: the box girder with a single box and two chambers (three webs with a web spacing of 740 cm, which corresponds to model CHB-2), the box girder with a single box and three chambers (four webs with a web spacing of 483 cm, which corresponds to model CHB-3), the box girder with a single box and four chambers (five webs with a web spacing of 385 cm, which corresponds to model CHB-4), and the box girder with a single box and five chambers (six webs with a web spacing of 300 cm, which corresponds to model CHB-5).

2.2.3. Models with Different Thicknesses of the Web

Based on model CHB-4, a comparative analysis of four models with the web thickness as the variable was performed using a three-bearing form with a support distance of 7.1 m. The four models are introduced as follows: model WebThk-35(55), with various web thicknesses of 35–55 cm; model WebThk-40(60), with web thickness of 40–60 cm; model WebThk-45(65) with web thickness of 45–65 cm; and model WebThk-50(70) with web thickness of 50–70 cm.

2.2.4. Models with Different Stiffnesses of End Diaphragm

In order to obtain the influence of stiffness of end diaphragm on the shear force distribution, six analytical models, model DPH-0, model DPH-0.4, model DPH-0.8, model DPH-1.0, model DPH-1.2, and model DPH-1.5, are proposed, corresponding to thicknesses of 0.0 m, 0.4 m, 0.8 m, 1.0 m, 1.2 m, and 1.5 m, respectively.

2.3. Results and Discussion to Analyze the Shear Distribution of Webs

The names of each web are stipulated with a single-box five-chamber configuration as an example (Figure 5). The two outermost webs are named side-webs (S-Webs), the webs close to S-Webs are called secondary side-webs (SS-Webs), and middle-webs (M-Webs) are at the center of the cross section.

Twenty groups of cross sections of the box girder were taken at intervals of (1.3 + 18 × 1.0 + 1.3) m from the midpoint of the bridge midspan to both sides of the bearing along the portrait direction of the bridge to analyze the shear distribution of each web in the portrait direction.

2.3.1. Shear Distribution of Webs with Different Support Cases

Model SPT-2, model SPT-3, and model SPT-4 are the models of the box girder with different bearing spacings, whose form of cross section is single box and five chamber. The roof width of the box girder is 20.4 m, the floor width is 15.4 m, and the center distance of chambers is 300 cm.

According to the conventional design method, the shearing force shared by each web is averaged using formula (1). However, a preliminary study showed that the shear force is not evenly distributed in each web; this paper defined the shear-increased coefficient η of each web, which characterizes the ratio of the actual shear force value of each web to the mean value of the webs, as shown in formula (2):where is the nominal shear value equally distributed to each web (MPa); is the shear value of the full section of the girder (MPa); is the number of webs; and is the shear value of single web (MPa).

From the 20 groups of cross sections in Figure 5, the shear forces shared by the S-Webs, SS-Webs, and M-Webs of the box girder in model SPT-2, model SPT-3, and model SPT-4 were extracted. The variation of the shear values and shear-increased coefficient η of each web from model SPT-2 to model SPT-4 along the portrait directions of the girder is shown in Figures 6(a)6(f).

The average value curve shows η = 1.0. As observed in Figure 6,(1)The shear forces shared by each web of the box girder are not evenly distributed, and the shear-increased coefficient η of both SS-Webs and M-Webs is less than 1.0. However, the shear-increased coefficient η of S-Webs is more than 1.0, which is unsafe in bridge design.(2)In model SPT-2, the shear-increased coefficient η of S-Webs is 1.58 in the midspan and 1.64 in the sections at the bearings, and the increase of these values is remarkable compared with the mean of 1.0. In addition, the shear-increased coefficient η of S-Webs in model SPT-3 and model SPT-4 is 1.31 and 1.27 in the section at the bearings, respectively, which cannot be ignored in the design.

The shear force shared by each web with double bearings is most different because the web near the bearings directly bears the reaction force from the bearings. According to Saint-Venant’s principle, the shear-increased coefficient η of each web gradually decreases from the bearings to the midpoint of the bridge.

2.3.2. Shear Distribution of Webs with Different Numbers of Chambers Cases

Model CHB-2, model CHB-3, model CHB-4, and model CHB-5 are the models of the box girder with different numbers of chambers, which used double bearings with a center distance of 14.2 m. The roof width of the box girder is 20.4 m, the floor width is 15.4 m, and its center distance of the chambers is 300 cm. Models CHB-2 to CHB-5 are described in detail in Section 2.2.2.

It can be seen from the calculation results that all shear-increased coefficients η of S-Webs in models CHB-2 to CHB-5 are greater than 1.0, and their values for SS-Webs and M-Webs are less than 1.0. As shown in Figure 7, with the increase in number of webs, the shear-increased coefficient η of S-Webs increases. In contrast, the shear-increased coefficient η of SS-Webs and M-Webs decreases with the increase in number of webs.

2.3.3. Shear Distribution with Different Thicknesses of the Web Cases

Model WebThk-35(55), model WebThk-40(60), model WebThk-45(65), and model WebThk-50(70) are the models of the box girder with different web thicknesses, whose cross-sectional forms are a single box with four chambers, three bearings, and a center distance of 7.1 m. The roof width of the box girder is 20.4 m, and the floor width is 15.4 m. Models WebThk-35(55) to WebThk-50(70) are described in Section 2.2.3. The results of different thicknesses of the web are shown in Figure 8.

The shear distribution results are as follows: the largest shear value is assigned to S-Webs, and the shear-increased coefficient η of S-Webs in models WebThk-35(55) to WebThk-50(70) is greater than 1.0. As shown in Figure 8, the shear-increased coefficient η values differ by less than 1%. Therefore, the effect of the web thickness on the shear distribution of the web is unremarkable.

2.3.4. Shear Distribution with Different Stiffnesses of End Diaphragm Cases

Both ends of the box beam are generally provided with end diaphragms of a certain thickness.

A certain number of supports are arranged below the end of the box beam, and the reaction force at the support has a longitudinal and transverse effect on the main beam. End diaphragms at the end of the box beam can play a transverse role against the reaction force of support, and at the same time, due to the action of end diaphragm, the reaction force of support can be transmitted to the webs relatively uniform. Therefore, whether end diaphragms are set or not, the stiffness of end diaphragm to the shear-increase coefficient η is worthy of further study. Six models of end diaphragms were studied in the effects of stiffness. Calculations prove that(1)When end diaphragms are not set, the shear-increased coefficient η of the S-Webs is the largest, reaching 1.37 at the internal forces extraction section 10 and 1.24 at the middle of the span (corresponding to internal forces extraction section 0).(2)After end diaphragms are set, the shear-increased coefficient η of the S-Webs is rapidly decreased. For example, for section 10, when the thickness of end diaphragms is 0.4 m, the shear-increased coefficient η of the S-Webs is 1.12, which is reduced by 18.2% with respect to the case where end diaphragms are not provided.(3)Within the range of thicknesses commonly used for end diaphragms, that is, in the range of 1.0∼1.5 m, the shear-increased coefficient η of the S-Webs is not obvious, and the difference is less than 1% (Figure 9).

3. Analysis of the Shear Failure of the SB-MC Box Girder

For the uneven distribution of shear force in the webs of an SB-MC box girder, a three-dimensional nonlinear model was established. Through the nonlinear analysis from loading to failure, we attempt to find the effect of the uneven distribution shear on the damage process of the webs and the redistribution of shear values in the webs.

3.1. Finite Element Model for the Nonlinear Analysis of the SB-MC Box Girder
3.1.1. Constitutive Tensile Relationship of Concrete Material

The nonlinear properties of the concrete material are considered in the Midas/FEA analysis to determine the causes of cracking. The total strain crack model in the concrete material parameters is a type of discrete crack model, which finds the crack behavior according to a fixed crack model or a rotating crack model. Once determined, the cracking direction does not change in any case for the fixed crack model; although less widespread, it was successfully applied in some problems of reinforced concrete structures [1820]. In comparison, the rotating crack model requires no knowledge of past cracking states, so its calculation process is much simpler and convergent [2123]. In this study, the rotating crack model was selected.

Based on fracture energy theory, a parabolic hardening-softening model is used in the total stress cracking model, which depends on three parameters of compressive strength fc, compressive fracture energy Gc, and characteristic element length h, as shown in Figure 10.

The strain corresponding to 1/3 of the maximum compressive strength fc is shown in equation (3), and the strain corresponding to the maximum compressive strength fc is shown in equation (4):

As can be seen from the above formula, the values are independent of cell size and compressive fracture energy. Thus, the ultimate strain at the softening end position of the compression zone is shown in the following equation:

A linear softening model is used in the total stress cracking model, which exhibits softening when the tensile strength of the base material is exceeded, as shown in Figure 11.

The major constitutive equations for the linearized softening tensile model are shown in the following equations:where is the tensile fracture energy of concrete, is the elastic modulus of concrete, is the crack width, is the tensile stress in the concrete, and is the minimum ultimate crack strain.

3.1.2. Constitutive Relationship of the Steel Bar

The von Mises model is widely used in the analysis of metallic materials. The von Mises model in MIDAS/FEA requires the definition of the hardening curve of the material yield stress fy,r. Tri-fold line stress-strain model of the hot-rolled steel bar was used to generate the constitutive model of the steel bar, as shown in Figure 12.

3.1.3. Bonding Model of Reinforced Concrete

Because the slippage between concrete and steel in reinforced concrete structures is impressed by cracks, a bond-slip model based on the deformation theory of plasticity can be used to simulate this effect. Assuming that the tangential direction is nonlinear and the normal direction of materials is linear elasticity, the bond-slip model can be expressed as follows:

By determining the right side of equations (7) and (8) for the relative displacement, the tangent stiffness is as follows:

3.1.4. A Model for Calculation

The object of nonlinear analysis is a box girder with a single box and four chambers. The basic parameters of the box girder are introduced as follows. Its span is 25 m, the beam height is 1.60 m, the roof and floor are 20 cm thick, and the thickness of the web from its midspan to the side-span is 40–60 cm in a linear gradient relationship. The middle segment of the girder is 10 m long, the segment with widening webs is 2.5 m long, the segment with widened webs is 3.8 m long, and end diaphragm is 1.2 m long. The distance between the chambers is 300 cm, the width of the roof is 20.4 m, and the width of the floor is 15.4 m. The model used double bearings with a bearing spacing of 14.2 m.

The stirrups are 12 mm in diameter, and the diameter of longitudinal steel bars of the bottom plate, the top plate, and the other is 32 mm, 16 mm, and 12 mm, respectively. Reinforcement configuration of the stirrups and longitudinal steel bars is shown in Figure 13.

The shear-span ratio of the model is 2.0. To avoid the stress concentration at the areas for loading and supporting, elastic blocks were used in the corresponding position, as shown in Figure 13. The sum of the applied vertical loads on a single loaded block is 70500 kN. Because the roof of the model is 20.4 m wide and the elastic block for loading is 0.475 m wide, the surface pressure applied to the surface of the elastic blocks is approximately 7275.5 kN/m2. The displacement criterion was used as the convergence condition and a modified Newton–Raphson iteration method was used to analyze the model in nonlinearity.

3.2. Results and Discussion of the Nonlinear Analysis
3.2.1. Curves of Load-Displacement

The vertical displacements of each web in the underside of the beam at the midspan were first extracted, as shown in Figure 14. Through the function multistep isosurface in the postprocessing of MIDAS/FEA, the crack development process of the box girder with a single box and four chambers was extracted according to the load step. The failure mode and cracking development of the beam are shown in Figures 15 and 16.

Figure 14 shows the following: (1) the entire loading process went through three phases (Stage I, Stage II, and Stage III) and three feature points (a, b, and c). In Stage I, corresponding to section 0‐a in Figure 14 for the elastic stage, the displacement increases linearly with the load; after the cracking characteristic point a (load factor is 0.22), the girder enters the cracking stage, corresponding to Stage II, as seen in Figure 16, the S-Webs firstly crack (Figure 16(a)), and then the SS-Webs and M-Webs sequentially crack (Figures 16(b) and 16(c)). As the load increases, the degree of cracking of the girder increases (Figures 15 and 16(d)). After reaching point b, the girder enters Stage III (the stirrups in shear span (corresponding to sections 8, 9, and 10, as seen Figure 5) and longitudinal rebars in the midsole plate yield) and further reaches the shear capacity.

3.2.2. Shear Force Redistribution during Loading

MIDAS/can extract not only the internal force of a solid element in the elastic phase but also the cracking phase. Under each load step, the shear force of the S-Webs, SS-Webs, and M-Webs of section no. 9 (Figure 5) was extracted, and the shear-increased coefficient η was calculated. Figure 17(a) shows the variation of the shear-increased coefficient η with loading during the whole process. Figure 17(b) clearly shows the change of the shear-increased coefficient η when the loading force was from 10000 to 30000 kN before and after cracking. The characteristics of the change in the shear-increased coefficient η are summarized as follows:(1)In the elastic phase I (corresponding to oa), the shear-increased coefficient η remains basically unchanged and is 1.25, 0.84, and 0.81 to S-Web, SS-Web, and M-Web, respectively.(2)When loaded to point a, S-Web first cracks; its shear-increased coefficient η reduces to 1.19, but the shear-increased coefficient η of SS-Web and M-Web increase to 0.9 and 0.82, respectively.(3)With the further loading, the girder continues to crack, and the difference of the shear-increased coefficient η among the webs is gradually reduced. When loaded to the bearing capacity of 70500 kN, the shear-increased coefficient η of S-Web, SS-Web, and M-Web is 1.09, 0.95, and 0.93, respectively.

The distribution of the shear-increased coefficient η in each web is a process of redistribution with the change of structural stiffness.

4. Design Measures considering the Laws of Shear Distribution

Through the comparative analysis of the results of the elastic model, we can preliminarily conclude that the shear distribution of webs in a box girder with a single box and multiple chambers is greatly affected by the supporting conditions and number of chambers, but shear-increased coefficient η is not sensitive to web thickness; the shear value shared by the S-Webs is larger than those shared by the SS-Webs and M-Webs. Therefore, the even distribution of shear force of the full section to each web in the design is unsafe to S-Webs. Nonlinear analysis shows that in the elastic phase, shear-increased coefficient η basically does not change, and coefficient η of each web changes with the degree of its cracking in the cracking stage, and S-Webs reach the ultimate load first.

The shear bearing capacity of the box girder is greatly affected by the thickness of its web [11]. Because the S-Webs have absorbed more shear force than the other webs, the shear-increased coefficient η is hardly affected by the web thickness and the main beam in the conventional design is generally in a nonshear crack state; this study proposes the design idea that according to the calculation result of shear-increased coefficient η, the S-Webs are thickened accordingly. This change may make the shear strength of each web match its shear resistance, and the webs almost simultaneously reach their shear capacity.

The formula can be expressed aswhere is the S-Webs thickness adjusted according to ; is the shear-increased coefficient calculated under the average web thickness, and is the average thickness of the webs before adjustment.

Model WebThk-40(60) of the box girder with a single box and four chambers is used as an example. The thickness of the web linearly changes from 40 cm at the middle segment to 60 cm at the widened webs segment. The shear-increased coefficient η of the S-Webs is 1.19, and the η of the SS-Webs and M-Webs is 0.83 and 0.97, respectively. According to the calculation results of η0, the thickness of the S-Webs is adjusted in 40 × η0 = 47.6 cm and 60 × η0 = 71.4 cm. According to the conventional web thickness value, the thickness of the S-Webs is changed to 50 cm and 70 cm, which is 1.25 times of the original thickness. For the webs with shear-increased coefcient η0 less than 1.0, the web thickness is not adjusted.

After the optimal design of S-Web thickness was obtained, the shear-increased coefficient ηs = 1.24 of S-Webs was again extracted, which matches the time value 1.24 of the adjusted thickness to the original thickness, in the same time matches the shear resistance of the S-Webs. On the other hand, ηs for the SS-Webs and M-Webs is 0.80 and 0.93, respectively, which matches their thickness and shear resistance. Therefore, the adjusted webs in the shear distribution are reasonable and safe.

5. Conclusions

In this paper, in the elastic model, the shear distribution of the SB-MC box girder was examined for different support conditions, the number of webs, the stiffness of end diaphragm, and the thickness of the web. In addition, the problem of shear redistribution was also studied in the nonlinear model. Finally, according to the shear-increased coefficient η of the S-Web, the method of structural optimization design is proposed. The main conclusions are as follows:(i)The Midas/FEA software, with great postprocessing capabilities, can easily and quickly extract the internal forces of models for a comparative analysis.(ii)The laws of shear distribution of the elastic model of the SB-MC box girder are as follows:(1)All the calculation conditions show that the S-Webs have more or less shared shear force greater than the average value. This means that it is not safe or even dangerous to design according to all webs bearing the average shear force.(2)When calculating the shear-increased coefcient η using a three-dimensional model, the end diaphragm must be built and included in the three-dimensional model; within the range of thicknesses commonly used for end diaphragms, that is, in the range of 1.0∼1.5 m, the shear-increased coefficient η of the S-Webs is not obvious.(3)Support conditions have an obvious influence on the shear distribution. The less support corresponds to the greater shear-increased coefficient η. In the model SPT-2, the shear-increased coefficient η of the S-Webs reaches 1.64. In addition, the number of the chambers is also an important factor affecting the coefficient η, and the irregularity of the shear distribution increases with the chamber number of the girder; the coefficient η of the S-Webs even reaches 1.65 in model CHB-5.(4)The change in web thickness has little effect on the shear distribution of the box girder.(iii)Nonlinear analysis shows that in the elastic phase, shear-increased coefficient η basically does not change, the coefficient η value of each web changes with the degree of web cracking in the cracking stage, and S-Webs reach the ultimate load first.(iv)The concept of the shear-increased coefficient η of the web was proposed, and the method to adjust the S-Web thickness according to the calculation result of η0 was introduced, which basically makes the shear forces of each web simultaneously reach their shear resistance; thus, the advanced shear failure in webs is prevented.

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 there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors acknowledge the financial support from the project of the National Key R&D Program of China (grant nos. 2018YFC0809600 and 2018YFC0809606) and the MOE Key Lab of Disaster Forecast and Control in Engineering of Jinan University (grant no. 20180930004).