Analysis of the Fatigue Crack Evolution of Corrugated Web Girders

Based on linear elastic fracture mechanics (LEFM), the fatigue crack evolution process and behavior of corrugated web girders were studied. The global finite element analysis (FEA) model of corrugated web girders was first developed and the equivalent structural stress method was used to reveal the dangerous locations along the weld under the bending load. The weld toe between the tension flange and the web weld, which is near the intersection of the inclined fold and the parallel fold, was determined as the fatigue crack easy-initiating location. Then a small region containing the crack-prone site was extracted as the sub-model for a crack propagating simulation. A semi-circle initial crack with 0.1 mm radius was inserted at the crack easy-initiating location. The stress intensity factors (SIFs; KI, KII, and KIII) of some discrete points along the crack front were calculated by the M-integral method. Based on Nasgro law, the geometry of the new crack front with a given extension length was obtained. Finally, the complete evolution process of the crack propagation was simulated. Results showed that the dominant crack propagating mode is open type (Mode I) and KI is the most important propagating driving force. According to the crack front shape evolution, the whole propagating process was divided into 6 stages. An obvious kink of the crack was found in stage 1, which covered only a very short time. The stages 3, 4 and 5 accounted for the majority of life, among which the stage 3 accounted for as high as 81% of the total life. Therefore, the cycles of the weld toe crack propagating from 0.1 mm to the thickness of the flange can be defined as the prediction life of this kind of structures.


Introduction
Girders with corrugated webs are improved components with some favorable mechanical properties and have been used in bridges and construction industries in recent years. Due to periodically changed corrugation geometry, corrugated web girders have outstanding shear capacity and can carry larger ultimate loads with thinner webs [1][2][3]. The ratio of web height to thickness can be as high as 600 and the buckling resistance can be 2.5-3 times of the flat web girders under the same load conditions [4][5][6][7][8]. With such excellent stability, corrugated web girders can reduce or eliminate the need for transverse stiffeners and the associated fabrication, which can greatly decrease the amount of steel used. Also due to the web corrugation, the longitudinal stress in the flanges is redistributed and results in local compressive stress which causes additional lateral bending [9][10][11][12] and torsion [13][14][15] in some local areas of flanges. However, corrugated web girders can still provide satisfactory bending performance [16][17][18] by using thicker flanges.
From the aspect of manufacturing, corrugated web girders are usually welded structures with plenty of weld seams. Working usually under stochastic cyclic load, the fatigue performance becomes an important issue in the design of corrugated web welded girders. Some investigation about the fatigue life and strength of the weld between the corrugated web and the flange of corrugated web girders have been carried out. Korashy [19] conducted comparative fatigue tests on corrugated web girders and conventional flat web girders, and the results suggested that corrugated web girders provide better fatigue performance. Anami [20,21] and Wang [22][23][24][25] carried out several fatigue tests on corrugated web girders and both thought that fatigue cracks were prone to initiate from the weld toe between the bottom flange and the corrugated web. They also proposed some practical fatigue life prediction methods. Ibrahim [26] theoretically studied the influence of corrugation geometry parameters on fatigue life and proposed that corrugation angle θ and bend radius r c of the transition between the inclined web fold and the longitudinal web fold were the two most important parameters. These studies were limited to a macroscopic perspective and focused on the relationship between loads S and fatigue life N. Of course, some factors such as the welded joints category and corrugation dimensions were also considered. However, the fatigue life in these studies was defined as the number of cycles at the end of the test, which generally meant that the testing girders were totally fractured. Such values usually overestimate the fatigue life and are of bad reliability for the design of corrugated web girders.
In the background of fatigue theory, the evolution of fatigue cracks can generally be divided into three typical stages: Crack initiation, crack propagation, and final fracture. For welded structures, crack-like defects are usually considered to be formed in the weld during the weld process. So the crack initiation stage of welded structures is normally assumed to be short enough to be ignored. Similarly, the final fracture stage is short and sudden and can also be ignored. Thus, the life of welded structures can be defined as the number of cycles of crack propagating from crack-like defects to critical cracks, which can be calculated by fracture mechanics [27]. However, with the rapid development of weld technology, it is controversial whether welded structures, especially those with good weld quality, contain weld defects equivalent to cracks. Some investigations have suggested that high-quality welds contain initial crack-like defects with depth of about 0.01-0.12 mm [28], 0.12-0.15 mm [29] or 0.01-0.4 mm [30]. These cracks are of short crack scale and have different behavior from macroscopic cracks. Life determined by simply integrating the crack propagation rate model, such as Paris-Erdogan law, will be underestimated.
In this work, the global FEA model was developed for a real corrugated web girder. The equivalent structural stress along the weld toe was calculated to determine the weakest fatigue location in the girder, at which the sub-model was extracted to obtain fine mesh and an initial crack was inserted in this sub-model. Based on Nasgro law, the crack front geometries at several given steps were calculated and the complete evolution process of the crack was simulated. According to the simulation results, some behaviors of the crack propagation were discussed.

Conventional Approaches to Predicting Life for Welded Joints
The fatigue cracks of welded structures generally originate from local high-stress areas such as the weld toe or weld root. At present, there are two common methodologies for fatigue life prediction of welded structures. One is the S-N curve methodology which is based on the welded joints categories [31][32][33]. Fatigue life is predicted in terms of various characteristic stress such as nominal stress, hot spot stress, equivalent structural stress, notched stress, peak stress, etc. The other one is the fracture mechanics methodology [34,35] whose theory basis is the relationship between the crack propagation rate and the stress intensity factor (SIF). The most famous model is the Paris-Erdogan law, as shown in Equation (1): where da/dN is the crack propagation rate; ∆K is the range of SIF; C and m are the material parameters. Paris-Erdogan law can be integrated to get the number of stress cycles required for the crack to propagate to a certain size, as shown in Equation (2): where a 0 is the initial crack length; a c is the crack propagating length; N c is the number of cycles when the crack propagates from a 0 to a c . Thus, fatigue life is not only determined by load (∆K), but also related to crack length (a 0 , a c ) so the fracture mechanics methodology is especially fit for the fatigue life prediction problems concerning objects with existing cracks, such as welded structures. Experiments show that crack propagating is a much more complex process which is affected by many factors except for load and crack length. For example, mean stress plays an important role in the crack propagating life, while different crack dimension scale drives the crack to propagate in different behavior. This results in some modified propagation rate models, such as Walker law, Forman law, McEvily law and Nasgro law. Although these models have different forms, they are basically based on the Paris model. After comprehensive consideration of the effects of stress ratio, threshold value, and fracture toughness, etc, they were corrected and improved. In general, the more parameters in the formula, the higher the accuracy of the simulation, but the more difficult it is to control, the greater the limitations. Among these laws, Nasgro law not only considers the influence of stress ratio and crack opening effect on crack propagation rate, but also control of the threshold of SIF range (∆K th ) on the crack beginning to propagate and the material fracture toughness (K C ) on the final fracture, as shown in Equation (3): where ƒ is the crack opening function; R is the stress ratio; K max is the maximum SIF; and C, n, p, q are material-related parameters. Figure 1 shows a schematic diagram of the principle of crack propagation simulation. Firstly, the current crack front is discretized into a series of points, at which SIF are all calculated. Secondly, considering that the crack front extension at discretized points are proportional to the SIFs increase at the same point, the extension length ∆a i is computed according to the given extension length ∆a median at the median point with the corresponding ∆K and ∆K i . Finally, having obtained the extension length at all the discrete points, the predicted crack front can be fitted by spline functions.

Approach to Simulating Crack Propagation
The formula for computing the propagation extension at discretized points is shown as Equation (4): where ∆a i is the crack front extension length at the i-th discretized point; R i is the stress ratio at the i-th discretized point; R median is the stress ratio at the median point.
The crack front propagation length ∆a i can also be obtained by specifying the load cycles step ∆N i , as shown in Equation (5): Figure 1 shows a schematic diagram of the principle of crack propagation simulation. Firstly, the current crack front is discretized into a series of points, at which SIF are all calculated. Secondly, considering that the crack front extension at discretized points are proportional to the SIFs increase at the same point, the extension length Δai is computed according to the given extension length Δamedian at the median point with the corresponding ΔK and ΔKi. Finally, having obtained the extension length at all the discrete points, the predicted crack front can be fitted by spline functions.

Initial Crack Location
In the weld cross-section, fatigue crack generally initiates from weld toe or weld root. However, the crack sites are affected by weld geometry and load types. Crack initiation sites are usually determined by the stress or strain distribution in the weld, in which principal stress, Von Mises stress, plastic strain, hot spot stress, and equivalent structural stress (ESS) can all be employed. ESS method is a new fatigue life predicting method proposed by Dong [36]. It accounts for a large number of fatigue tests data for welded joints and structures and is proved to be sensitive to the local fatigue damage of the weld toe. In this work, ESS method was employed to predict the crack initiation sites of a corrugated welded girder.
ESS method uses the mesh-insensitive structural stress to characterize the local stress state of the weld detail, which is constructed by the line force and the line moment of the nodes at the weld toe according to structural mechanics, as shown in Equation (6): where σ s is the structural stress; σ m is the membrane stress; σ b is the bending stress; ƒ y is the line force; m x is the line moment; t is the plate thickness. ESS method also accounts for the influence of plate thickness and loading condition on fatigue life. Based on the linear elastic fracture mechanics (LEFM), structural stress can be corrected to equivalent structural stress, as shown in Equation (7): where S S is the equivalent structural stress; m is the slope of the crack propagation curve, here m = 3.6; I(r) is the loading mode correction function; r is the bending ratio, the expression is shown in Equation (8): Using finite element analysis (FEA), ESS along the entire weld toe curve can be calculated and drawn. The location with maximum value can be determined as the potential crack initiation, in which an initial crack can be modelled.
Another advantage of ESS method is that fatigue life can also be estimated from only one S-N curve for various welded joints, instead of judging proper S-N curves for specific welded joints to obtain high prediction accuracy.

Crack Propagation Driving Force
In the theory of LEFM, the commonly used parameters of crack propagation are Stress Intensity Factor (SIF, K), J-integral (J), energy release rate (G) and crack opening displacement (COD). These parameters are equivalent in the stage of linear elasticity. Among those, K is the basic parameter describing the stress field intensity near the crack front. In terms of the stress, strain, and displacement state near the crack front, cracks have three modes: (a) Mode I (opening mode), (b) Mode II (in-plane shear mode), and (c) Mode III (anti-plane shear mode). During the crack propagating, the crack front changes gradually even kink, which leads to a complex shape and propagating path. In this case, it is necessary to synthesize the mixed effect of the three modes. Equivalent SIF K equ is to be used here, as shown in Equation (9): where K equ is the equivalent SIF; K I , K II , K III are the SIFs under corresponding modes; ν is the Poisson's ratio. SIF can generally be solved by analytical methods, numerical methods, and experimental calibration methods. Numerical methods such as the FEA method are often used in practical applications due to the complex structures and the stress-strain state in the weld region. However, the FEA method cannot deal with crack geometry in the normal way and cannot avoid the discontinuity and singularity of the crack front. So, it is difficult to calculate the SIF directly. In practice, the J-integral method is often used to solve SIF indirectly, which is defined as Equation (10): where J is the J-integral; Г is the integral loop around the crack front; σ ij is the stress tensor; u i is the component of the displacement vector; W = σ ij ε ij is the strain energy density; δ 1j is the kronic function; q is the auxiliary function; x j is the integral loop coordinate; ds is a small increment on the integration path. Traditional J-integral definition cannot distinguish K I , K II and K III . So Knowles [37] proposed the M-integral method, by which K I , K II and K III can be separated and acquired by introducing an auxiliary field at the crack front. In the M-integral method, real field is the actual displacement field and stress field at the crack front, while the auxiliary fields are the any possible displacement field and stress field that can satisfy the equilibrium condition and the physical equation. Considering that stress, strain, and displacement are linear elastic and can satisfy the superposition principle, they can be separated as shown in Equation (11), where the superscripts (1) and (2) represent the real field state and the auxiliary field state, respectively.
Note that the third term in Equation (12) is M-integral, which can be shown as Equation (13): Compared with J-integral, M-integral can separate and solve the SIFs of three modes by choosing a special auxiliary field. The relationship between M-integral and K (that is SIF) is shown as Equation 14: where E is the elastic modulus of the material. In order to construct special auxiliary, the following analytic formulas can be used, as shown in Equation (15): Then, K I , K II , and K III can be calculated by introducing the analytical formulas in Equation (15) into Equation (14), as shown in Equation (16):

Crack Kink Criterion
In the crack propagation simulation, the crack propagating direction has an important influence on the crack shape evolution. For Mode I, the crack will propagate along its original crack plane. In the mixed-mode, the crack front is likely to kink and the local propagating direction of the crack is affected by the kink angle θ. Kink angle is usually defined as the angle by which the new crack surface kinks from the previous one. It is worth noting that kink angle θ is closely related to the stress state of the crack front.
There are four kink criteria representing the relationship between the kink angle θ and the SIF K.

MTS (Max Tensile Stress) criterion
This criterion states that the crack will propagate in the direction of maximum circumferential stress. That is, when the circumferential normal stress of a discrete point on the crack front is maximum, the corresponding θ at that discrete point is taken as the crack front kink angle θ kink . The expression is shown as Equation (17):

MSS (Max Shear Stress) criterion
This criterion states that the crack front is determined by shear stress. Then, the effects of shear components of modes II and III are introduced by the weight η II and η III , as shown in Equation (18):

GEN (Generalized Stress) criterion
This criterion presents the competition mechanism of circumferential normal stress and shear stress at the crack front. It assumes that the crack will propagate in the direction with the larger control factor from the above two criteria. The expression is shown in Equation (19):

SERR (Strain Energy Release Rate) criterion
This criterion states that the crack will propagate in the direction with maximum strain energy release rate. In essence, this criterion is based on the energy viewpoint and can account for the comprehensive effect of K I , K II , and K III . The calculation expression is shown as Equation (20):

Case Description
Tests and results in [38] were referred to and discussed in this work, in which six large-scale corrugated web girder specimens were fabricated and fatigue tests were conducted under cyclic loading. The six girder specimens had the same corrugation profile with different loading conditions and weld sizes. The girder specimens and loading equipment can be seen in Figure 2. Among the six specimens, three experienced fatigue failure after a certain number of cycles. Cracks were all observed close to the intersection of the bending radius section and inclined section, which were at the top surface of the lower flange in the middle part of the span.  Among the six specimens, the third specimen was simulated and investigate in this work, whose profile is shown in Figure 3. According to the data given in [38], the total length of the corrugated web girder was 7150 mm and the span was 6750 mm. The web thickness was 6 mm and its depth was 500 mm. The flange thickness was 20 mm and its width was 225 mm. The length of one periodic trapezoidal corrugation was 750 mm, with parallel section length of 210 mm, inclined section length of 165 mm, bending radius of 60 mm and a corrugation angle of 39 • . The top flange, lower flange, and the web were connected by four welds with a weld size of 6 mm. There were three stiffeners between the top flange and the lower flange to improve the local stability at the loading point and the support points. The height of the stiffeners at the loading point and the support points at both ends were 500 mm and 375 mm, respectively, with a thickness of 20mm. The fatigue test was carried out with minimum load F min = 10 kN, maximum load F max = 240 kN and the load ratio R was 0.042. Note that N and S are marked in Figure 3 to indicate the north side and south side of the corrugated web girder, respectively, which will be used in the following simulation analysis.

Proposed Simulation Procedure
Firstly, the FEA model of the corrugated web girder with shell elements was established, in which the weld was modeled by a layer of shell elements whose dimensions were determined according to the ESS definition. Secondly, the ESS at the lower flange weld toe was calculated and the crack initiation location was analyzed in terms of the ESS distribution. Thirdly, the global 3D FEA model of the same specimen was established, with which the crack initiation region was extracted as the sub-model. The sub-model was associated with the global model through boundary displacement conditions and a semi-circle initial surface crack was inserted into the sub-model at the crack initiation location with the major diameter along the weld toe direction and the minor diameter along the thickness direction. Fifthly, extension length at every discrete point of the crack front was calculated on the basis of the principle described in Section 2.2 and the new crack front curve was fitted. Then a new crack with the new front curve was inserted into the sub-model again, and the same computing procedure was carried out again. In this way, the crack propagation process was simulated until the specimen could not bear the testing load any longer and finally fractured. The complete crack propagation simulation flow chart was shown in Figure 4. In the simulation,

Proposed Simulation Procedure
Firstly, the FEA model of the corrugated web girder with shell elements was established, in which the weld was modeled by a layer of shell elements whose dimensions were determined according to the ESS definition. Secondly, the ESS at the lower flange weld toe was calculated and the crack initiation location was analyzed in terms of the ESS distribution. Thirdly, the global 3D FEA model of the same specimen was established, with which the crack initiation region was extracted as the sub-model. The sub-model was associated with the global model through boundary displacement conditions and a semi-circle initial surface crack was inserted into the sub-model at the crack initiation location with the major diameter along the weld toe direction and the minor diameter along the thickness direction. Fifthly, extension length at every discrete point of the crack front was calculated on the basis of the principle described in Section 2.2 and the new crack front curve was fitted. Then a new crack with the new front curve was inserted into the sub-model again, and the same computing procedure was carried out again. In this way, the crack propagation process was simulated until the specimen could not bear the testing load any longer and finally fractured. The complete crack propagation simulation flow chart was shown in Figure 4. In the simulation, coordination, kink angle, and SIF value of the crack front discretized points generated in every step were extracted and were used in the following analysis and discussion. coordination, kink angle, and SIF value of the crack front discretized points generated in every step were extracted and were used in the following analysis and discussion.

Predicted Crack Initiation Location
The ESS method was used to predict the crack initiation location on the weld toe. The FEA model of the corrugated web girder specimen, in which both flanges, corrugated web, stiffeners, and four welds were modeled as shell elements of S4 (in Abaqus), was established. The total elements number was 65069 and the total nodes number was 63078. The FEA model was constrained as simple supported girders and the concentrated force of F = 240 kN was applied in the middle of the girder span. It is worth noting that the weld was modeled as one layer of inclined shells between the flange and the corrugated web, whose nodes were exactly at the weld toe on the flange side and the web side, respectively, as shown in Figure 5. Modelling like this is in order to extract the node forces and the node moments at the weld toe nodes.

Predicted Crack Initiation Location
The ESS method was used to predict the crack initiation location on the weld toe. The FEA model of the corrugated web girder specimen, in which both flanges, corrugated web, stiffeners, and four welds were modeled as shell elements of S4 (in Abaqus), was established. The total elements number was 65069 and the total nodes number was 63078. The FEA model was constrained as simple supported girders and the concentrated force of F = 240 kN was applied in the middle of the girder span. It is worth noting that the weld was modeled as one layer of inclined shells between the flange and the corrugated web, whose nodes were exactly at the weld toe on the flange side and the web side, respectively, as shown in Figure 5. Modelling like this is in order to extract the node forces and the node moments at the weld toe nodes. coordination, kink angle, and SIF value of the crack front discretized points generated in every step were extracted and were used in the following analysis and discussion.

Predicted Crack Initiation Location
The ESS method was used to predict the crack initiation location on the weld toe. The FEA model of the corrugated web girder specimen, in which both flanges, corrugated web, stiffeners, and four welds were modeled as shell elements of S4 (in Abaqus), was established. The total elements number was 65069 and the total nodes number was 63078. The FEA model was constrained as simple supported girders and the concentrated force of F = 240 kN was applied in the middle of the girder span. It is worth noting that the weld was modeled as one layer of inclined shells between the flange and the corrugated web, whose nodes were exactly at the weld toe on the flange side and the web side, respectively, as shown in Figure 5. Modelling like this is in order to extract the node forces and the node moments at the weld toe nodes.   Figure 6 shows the ESS curve along the weld toes of the north side and south side on the lower flange, in which the blue solid line indicates the ESS of south side weld toe (ESS_S), the green solid line indicates the ESS of north side weld toe (ESS_N), the black solid line presents the corrugated shape of the web, and the vertical black dash line marks the middle of the girder span. It can be seen that, due to the folding effect of the corrugated web, ESS on both north and south weld toes exhibit periodic change symmetric to the middle line. No matter on the south or north side, the ESS on the inclined section is larger than that on the parallel section. The stress in the middle of the span is the largest. The farther away it is from the middle span, the smaller the maximum value of the stress is. ESS reaches the large values almost at the intersections of the inclined section and the bending radius section in the middle part of the span, which are all the potential fatigue weakest points of the specimen.
Metals 2019, 9, x FOR PEER REVIEW 10 of 22 Figure 6 shows the ESS curve along the weld toes of the north side and south side on the lower flange, in which the blue solid line indicates the ESS of south side weld toe (ESS_S), the green solid line indicates the ESS of north side weld toe (ESS_N), the black solid line presents the corrugated shape of the web, and the vertical black dash line marks the middle of the girder span. It can be seen that, due to the folding effect of the corrugated web, ESS on both north and south weld toes exhibit periodic change symmetric to the middle line. No matter on the south or north side, the ESS on the inclined section is larger than that on the parallel section. The stress in the middle of the span is the largest. The farther away it is from the middle span, the smaller the maximum value of the stress is. ESS reaches the large values almost at the intersections of the inclined section and the bending radius section in the middle part of the span, which are all the potential fatigue weakest points of the specimen. The ESS curve of the mid-span region is exaggeratedly shown and compared with the test results in [38]. Three specimens were found with fatigue cracks among the six tests, in which cracks occurred on the north side in the second bending radius section in specimens 3 and 4, and on the south side in the first bending radius section in specimen 5, as shown in Figure 7. It is obvious that the ESS at these locations were all at high levels, indicating that the ESS method has advantage of capturing the crack initiation location for welded structures. Note that cracks will not always occur in the location where the ESS is maximum. At some positions with larger ESS but lower than the maximum value, cracks also probably initiate according to the weld quality at the initiation sites, where the weld toe morphology should be paid attention to. In Figure 7, the ESS value corresponding to the fatigue failure location of specimen 3 is 106.84 MPa. Considering the stress ratio R is 0.042, fatigue life at this point can be predicted with different survival probabilities using the master S-N curve provided by Dong [39], which is expressed as Equation 21: The ESS curve of the mid-span region is exaggeratedly shown and compared with the test results in [38]. Three specimens were found with fatigue cracks among the six tests, in which cracks occurred on the north side in the second bending radius section in specimens 3 and 4, and on the south side in the first bending radius section in specimen 5, as shown in Figure 7. It is obvious that the ESS at these locations were all at high levels, indicating that the ESS method has advantage of capturing the crack initiation location for welded structures. Note that cracks will not always occur in the location where the ESS is maximum. At some positions with larger ESS but lower than the maximum value, cracks also probably initiate according to the weld quality at the initiation sites, where the weld toe morphology should be paid attention to.
Metals 2019, 9, x FOR PEER REVIEW 10 of 22 Figure 6 shows the ESS curve along the weld toes of the north side and south side on the lower flange, in which the blue solid line indicates the ESS of south side weld toe (ESS_S), the green solid line indicates the ESS of north side weld toe (ESS_N), the black solid line presents the corrugated shape of the web, and the vertical black dash line marks the middle of the girder span. It can be seen that, due to the folding effect of the corrugated web, ESS on both north and south weld toes exhibit periodic change symmetric to the middle line. No matter on the south or north side, the ESS on the inclined section is larger than that on the parallel section. The stress in the middle of the span is the largest. The farther away it is from the middle span, the smaller the maximum value of the stress is. ESS reaches the large values almost at the intersections of the inclined section and the bending radius section in the middle part of the span, which are all the potential fatigue weakest points of the specimen. The ESS curve of the mid-span region is exaggeratedly shown and compared with the test results in [38]. Three specimens were found with fatigue cracks among the six tests, in which cracks occurred on the north side in the second bending radius section in specimens 3 and 4, and on the south side in the first bending radius section in specimen 5, as shown in Figure 7. It is obvious that the ESS at these locations were all at high levels, indicating that the ESS method has advantage of capturing the crack initiation location for welded structures. Note that cracks will not always occur in the location where the ESS is maximum. At some positions with larger ESS but lower than the maximum value, cracks also probably initiate according to the weld quality at the initiation sites, where the weld toe morphology should be paid attention to. In Figure 7, the ESS value corresponding to the fatigue failure location of specimen 3 is 106.84 MPa. Considering the stress ratio R is 0.042, fatigue life at this point can be predicted with different survival probabilities using the master S-N curve provided by Dong [39], which is expressed as Equation 21: In Figure 7, the ESS value corresponding to the fatigue failure location of specimen 3 is 106.84 MPa. Considering the stress ratio R is 0.042, fatigue life at this point can be predicted with different survival probabilities using the master S-N curve provided by Dong [39], which is expressed as Equation (21): where N is the fatigue life; A and B are the property parameters, which are shown with different failure probabilities in Table 1.  Table 1 also lists the calculated fatigue life with different failure probabilities, among which the life with 50% failure probability is 968,100 cycles. The predicted life is much lower than the test life of 1,305,500 cycles given in [38], with an error about 26.07%. It can be included that the predicted life is not satisfyingly accurate, with at least 50% failure probability. In fact, if we take the test life as the calculated life, its failure probability is about 61.79%, which is absolutely unacceptable in engineering applications.
Therefore, it is suggested in this paper that, although the ESS method can capture the fatigue crack initiation location accurately, but it cannot be used to predict the fatigue life due to bad accuracy. Simulation based fracture mechanics method should be employed which can characterize the fatigue propagating behavior physically.

Modeling of Initial Crack
Initial crack size is the most important parameter for fatigue crack propagation simulation. Generally, initial defects near the weld toe in the weld such as cavities and pits can be regarded as short cracks, while some visible defects on the weld surface slag together with the starting and ending location of the weld can be regarded as micro cracks. All these defects have very complicated geometry and their scale involve several magnitude orders of 10 µm, 100 µm, and 1 mm, which directly impact on crack propagating behavior and fatigue life. Some investigations suggest that the weld toe surface crack usually evolves into a semi-elliptic shape when it propagates to a certain scale, no matter what the shape of the original defect. Considering this fact, the initial crack was set as a semicircular shape in this paper. Otherwise, IIW [40] recommends the initial weld toe crack depth to be 0.1 mm, which covers about 3-5 grains for common structural steel materials and can characterize some propagation behaviors for short cracks to some extent. So the initial crack in this paper was set with depth of 0.1 mm. Due to the very small crack size, a very fine mesh for the crack front with 3D solid elements was needed. If all of elements are of same size, the FEA model with a large number of elements will be generated, resulting in significant increase in computation time. Thus, sub-model technology was used to achieve faster calculation efficiency together with higher solution accuracy.
Some primary procedures in the simulation are shown in Figure 8, in which Figure 8a presents the global FEA model of the total girder. The elements are all C3D20R with the common size of 15 mm. The total elements number is 131012 and the total nodes number is 644018. Figure 8b shows the sub-model including a small part of lower flange, corrugated web, and the connecting weld, which is in the second corrugation section right to the mid-line of the girder. The sub-model covers the bending radius section and the inclined section, where an initial crack will be inserted at the potential crack site which is determined by the ESS method. The inserted crack is of semi-elliptic shape with length of 0.1 mm, which is parallel to the weld toe and perpendicular to the upper surface of the lower flange, as shown in Figure 8c. In order to make the element size of the initial crack front coincide with the element size of other regions, the global element size of the sub model is set to 1.5 mm, while the local element size near the crack front is set to 0.2 mm. Elements with gradually varied size were utilized at the transition part between the global model and the sub-model, as well as between the crack front and the region far away from the crack front in the sub-model, as shown in Figure 8d. It is worth noting that in every step in the simulation, the sub-model should be re-meshed in terms of the current crack front, and the mesh of the crack front must follow some modelling rules in order to eliminate the stress singularity at the crack front, as shown in Figure 8f. Three layers of elements were used to simulate the crack front, among which singular wedge elements with 15 nodes (C3D15) were used for the first layer elements to model the innermost circle of the crack front, and hexahedral elements with 20 nodes (C3D20) were used for modelling the second and third layers as two rings. In addition, pyramid elements with 13 nodes (C3D13) were used to connect the crack front region with the rest part of the sub-model. utilized at the transition part between the global model and the sub-model, as well as between the crack front and the region far away from the crack front in the sub-model, as shown in Figure 8d. It is worth noting that in every step in the simulation, the sub-model should be re-meshed in terms of the current crack front, and the mesh of the crack front must follow some modelling rules in order to eliminate the stress singularity at the crack front, as shown in Figure 8f. Three layers of elements were used to simulate the crack front, among which singular wedge elements with 15 nodes (C3D15) were used for the first layer elements to model the innermost circle of the crack front, and hexahedral elements with 20 nodes (C3D20) were used for modelling the second and third layers as two rings.
In addition, pyramid elements with 13 nodes (C3D13) were used to connect the crack front region with the rest part of the sub-model.

Simulation Parameters
In terms of [38], the material of the web and flanges of the corrugated web girder specimen was S355, which was a kind of low-alloy structural steel with basis mechanical properties shown in Table  2. In the crack propagation simulation, the crack extension length Δamedian is a vital parameter which

Simulation Parameters
In terms of [38], the material of the web and flanges of the corrugated web girder specimen was S355, which was a kind of low-alloy structural steel with basis mechanical properties shown in Table 2. Moreover, according to American Society for Testing and Materials code (ASTM), corresponding material parameters of grade S355 can be referred to grade A588, which has Nasgro parameters such as C = 7.3080 × 10 −14 , n = 3.3, p = 0.5, q = 0.5, ∆K th = 295.375 MPa·mm 1/2 , K c = 3475 MPa·mm 1/2 .
In the crack propagation simulation, the crack extension length ∆a median is a vital parameter which controls the simulation process. In the beginning of the simulation, the crack grew slowly, so the crack extension length was configured as 0.02 mm. In the subsequent process, the extension length was gradually increased according to the crack size and configured as 0.04 mm, 0.08 mm, 0.1 mm, 0.2 mm, 0.4 mm, 0.8 mm, and 1 mm respectively. When the crack reached the boundary of different geometry regions such as from flange to weld, special extension length should be configured in real time according to actual situation.
Several simulation schemes were considered to analyze the effects of crack propagation driving forces and crack kink criteria, as shown in Table 3. The simulation results show that the crack propagation patterns using MTS criterion and SERR criterion are both consistent with the test result, while those using MSS criterion and GEN criterion are different from the test, which shows that K I plays a leading role in crack propagation of corrugated web girders under three-point bending. Considering the mixed effect of three crack modes, the results of three schemes of K I + MTS, K I + SERR, and K equ + SERR are discussed in the following, as indicated by * in Table 3.

Crack Front Shape Evolution
It is noticeable that all the K I + MTS, K I + SERR and K equ + SERR schemes have similar simulation results. Therefore, the K equ + SERR simulation result is selected here and Figure 9a shows its simulated crack shape on the upper surface on the lower flange, while Figure 9b shows the corresponding test results. It can be seen that they are in good consistence. One end of the crack propagates toward the nearby edge of the lower flange until goes through it, while the other end propagates toward the weld and then the web and finally the other side of the flange. Viewed from a micro scale, the crack surface is almost perpendicular to the longitudinal direction of the girder, which shows the principal stress on the lower flange is the dominant factor driving the crack propagating.
It is worth noting that in the very early stage of the crack propagation simulation, distinct crack front kinks are observed in all three schemes. Figure 10 shows the first six simulation steps under the K equ + SERR scheme. It can be seen that the long diameter of the initial crack is originally parallel to the weld toe direction. However, after six steps of propagation, it gradually rotates to the direction almost parallel to the width of the flange. At the 6th step, the crack size scale is very small, with the depth in the thickness direction of about 0.09 mm and the length on the surface of about 0.35 mm. The shape ratio a/c varies from 1 to 0.7 and the variation of kink angle is about 49 • . After that, the crack propagates almost in the same plane. The above phenomenon was also observed in other investigations [20,22], which indicated that the macro crack of the corrugated web girders was almost 3~5 • with the longitudinal direction of the girders. Considering the corrugation angle of the specimens in this work is 39 • , the crack kink angle ranges from 46 • to 48 • , which is very close to the simulation result.
It is noticeable that all the KI + MTS, KI + SERR and Kequ + SERR schemes have similar simulation results. Therefore, the Kequ + SERR simulation result is selected here and Figure 9a shows its simulated crack shape on the upper surface on the lower flange, while Figure 9b shows the corresponding test results. It can be seen that they are in good consistence. One end of the crack propagates toward the nearby edge of the lower flange until goes through it, while the other end propagates toward the weld and then the web and finally the other side of the flange. Viewed from a micro scale, the crack surface is almost perpendicular to the longitudinal direction of the girder, which shows the principal stress on the lower flange is the dominant factor driving the crack propagating.  Kequ + SERR scheme. It can be seen that the long diameter of the initial crack is originally parallel to the weld toe direction. However, after six steps of propagation, it gradually rotates to the direction almost parallel to the width of the flange. At the 6th step, the crack size scale is very small, with the depth in the thickness direction of about 0.09 mm and the length on the surface of about 0.35 mm.
The shape ratio a/c varies from 1 to 0.7 and the variation of kink angle is about 49°. After that, the crack propagates almost in the same plane. The above phenomenon was also observed in other investigations [20,22], which indicated that the macro crack of the corrugated web girders was almost 3~5° with the longitudinal direction of the girders. Considering the corrugation angle of the specimens in this work is 39°, the crack kink angle ranges from 46° to 48°, which is very close to the simulation result. By extracting the geometry information of the crack front under all the steps and drawing the critical crack front, the crack front evolution can be divided into six stages, as shown in Figure 11 and Table 4.  By extracting the geometry information of the crack front under all the steps and drawing the critical crack front, the crack front evolution can be divided into six stages, as shown in Figure 11 and Table 4. Figure 11 shows that in stage 1 (in red), the crack propagates simultaneously in the weld toe direction and the depth direction, with the crack front shape gradually changing from semi-circle to semi-ellipse and the crack rapidly kinking from parallel to the weld toe direction to almost perpendicular to the longitudinal direction of the girder. In stage 2 (in orange), one end of the crack front near the weld approaches the weld and the other end continues to propagate on the upper surface of the lower flange to the north side. In stage 3 (in yellow), the end point near the weld enters the weld and subsequently the web, and the other end continues to propagate in the same manner as before. The crack shape of this stage is approximately semi-ellipse but slightly larger. In stage 4 (in green), the crack approaches the lower surface of the lower flange in the depth direction and is divided into two parts. It is worth noting that during this stage the crack front end near the weld enters the web and subsequently the south side weld. In stage 5 (in blue), the two cracks formed in the previous stage propagate mainly on the bottom flange as the edge cracks at the respective sides simultaneously, with both shapes change from approximately semi-ellipse to approximately line. Stage 6 is very short, in which the two edge cracks propagate in high rate until the north side edge crack tears the flange. Then the specimen completely loses bearing capacity and the simulation stops. By extracting the geometry information of the crack front under all the steps and drawing the critical crack front, the crack front evolution can be divided into six stages, as shown in Figure 11 and Table 4.

Characteristic Points along the Crack Front
During the crack propagation, the shape and the length of the crack front are changing constantly. In order to describe the distribution of the crack front angle and key parameters such as SIF, the crack front curves in all the simulation steps were calibrated in terms of curve coordinates, in which two endpoints are defined as 0 (the south end) and 1 (the north end), respectively, and the rest points are set coordinates linearly from 0 to 1, according to the ratio location of the point on the entire curve. It should be noted that when the crack front is just broken into two parts, the curve coordinate of the breakpoint is 0.57 according to the simulation results. This point can be regarded as the control point describing the crack propagating in the depth direction. All these three points are defined as characteristic points, as shown in Figure 12. Hereby, the crack front segment of 0-0.57 is defined as Crack Front 1, and the crack front segment of 0.57-1 is defined as Crack Front 2. In Figure 12

Kink Angle
Crack kinking lasts a very short time and only occurs in stage 1. The absolute kink angle of all discrete points along the crack front of different simulation steps are extracted and the kink angle distribution curves are drawn, as shown in Figure 13, which is under scheme Kequ + SERR. In this figure, macro crack surface location is defined as the reference plane with a kink angle of 0°. Results of the six steps in stage 1 and the first two steps in stage 2 are chosen for comparison. It can be seen that the kink angles of the discrete points along the crack front are symmetrically distributed to the central point (0.5) in terms of curve coordinates. The kink angles of the two endpoints are opposite in direction and almost the same in magnitude, which demonstrates the whole crack front is almost kept in the same plane. The kink angles of different steps change gradually until the crack turns to the macro crack direction. The calculated kink angles seem to fluctuate slightly on both sides of the crack, which probably result from the numerical interpolation computation method.

KI
Considering that different crack propagation driving forces and kink angle criteria may affect the initial behavior of the crack, the values of KI are extracted under the three simulation schemes of KI + MTS, KI + SERR, and Kequ + SERR in stage 1 and the distribution of KI is presented, as shown in Figure 14, in which the abscissa is the crack propagation depth. It can be seen that the overall variation trend of the same point is similar, with slight differences under different simulation schemes. The value of KI of CP_0 (location at the south end) is slightly greater than the other two characteristic points, which means that the propagating rate at this point is the largest. Meanwhile, CP_0.57 (location at the break point) has the lowest propagating rate, which means the crack may be hindered in the depth direction. It is the difference of the propagation rate at different points on the crack front that leads to the crack shape changing from semi-circle to semi-ellipse. In addition, from Figure 14, it can be seen that the curves of Kequ + SERR usually lie between those of KI + MTS and KI + SERR, so the results of Kequ + SERR are to be analyzed and discussed in the following.

Kink Angle
Crack kinking lasts a very short time and only occurs in stage 1. The absolute kink angle of all discrete points along the crack front of different simulation steps are extracted and the kink angle distribution curves are drawn, as shown in Figure 13, which is under scheme K equ + SERR. In this figure, macro crack surface location is defined as the reference plane with a kink angle of 0 • . Results of the six steps in stage 1 and the first two steps in stage 2 are chosen for comparison. It can be seen that the kink angles of the discrete points along the crack front are symmetrically distributed to the central point (0.5) in terms of curve coordinates. The kink angles of the two endpoints are opposite in direction and almost the same in magnitude, which demonstrates the whole crack front is almost kept in the same plane. The kink angles of different steps change gradually until the crack turns to the macro crack direction. The calculated kink angles seem to fluctuate slightly on both sides of the crack, which probably result from the numerical interpolation computation method.

Kink Angle
Crack kinking lasts a very short time and only occurs in stage 1. The absolute kink angle of all discrete points along the crack front of different simulation steps are extracted and the kink angle distribution curves are drawn, as shown in Figure 13, which is under scheme Kequ + SERR. In this figure, macro crack surface location is defined as the reference plane with a kink angle of 0°. Results of the six steps in stage 1 and the first two steps in stage 2 are chosen for comparison. It can be seen that the kink angles of the discrete points along the crack front are symmetrically distributed to the central point (0.5) in terms of curve coordinates. The kink angles of the two endpoints are opposite in direction and almost the same in magnitude, which demonstrates the whole crack front is almost kept in the same plane. The kink angles of different steps change gradually until the crack turns to the macro crack direction. The calculated kink angles seem to fluctuate slightly on both sides of the crack, which probably result from the numerical interpolation computation method.

KI
Considering that different crack propagation driving forces and kink angle criteria may affect the initial behavior of the crack, the values of KI are extracted under the three simulation schemes of KI + MTS, KI + SERR, and Kequ + SERR in stage 1 and the distribution of KI is presented, as shown in Figure 14, in which the abscissa is the crack propagation depth. It can be seen that the overall variation trend of the same point is similar, with slight differences under different simulation schemes. The value of KI of CP_0 (location at the south end) is slightly greater than the other two characteristic points, which means that the propagating rate at this point is the largest. Meanwhile, CP_0.57 (location at the break point) has the lowest propagating rate, which means the crack may be hindered in the depth direction. It is the difference of the propagation rate at different points on the crack front that leads to the crack shape changing from semi-circle to semi-ellipse. In addition, from Figure 14, it can be seen that the curves of Kequ + SERR usually lie between those of KI + MTS and KI + SERR, so the results of Kequ + SERR are to be analyzed and discussed in the following.

K I
Considering that different crack propagation driving forces and kink angle criteria may affect the initial behavior of the crack, the values of K I are extracted under the three simulation schemes of K I + MTS, K I + SERR, and K equ + SERR in stage 1 and the distribution of K I is presented, as shown in Figure 14, in which the abscissa is the crack propagation depth. It can be seen that the overall variation trend of the same point is similar, with slight differences under different simulation schemes. The value of K I of CP_0 (location at the south end) is slightly greater than the other two characteristic points, which means that the propagating rate at this point is the largest. Meanwhile, CP_0.57 (location at the break point) has the lowest propagating rate, which means the crack may be hindered in the depth direction. It is the difference of the propagation rate at different points on the crack front that leads to the crack shape changing from semi-circle to semi-ellipse. In addition, from Figure 14, it can be seen that the curves of K equ + SERR usually lie between those of K I + MTS and K I + SERR, so the results of K equ + SERR are to be analyzed and discussed in the following.  Figure 15 shows the distribution of KI of all the six stages with the abscissa of normalized coordinates of the crack front, in which several curves concerning different simulation steps may be of one stage with the same color. It can be seen that in stage 1and 2, the KI of the entire crack front is almost constant and at a low level, which is mainly due to the short crack situation at these two stages. When the crack grows into stage 3, the KI of CP_1 gradually increases, resulting in a KI larger than those of CP_0.57 and CP_0. This is mainly due to the fact that CP_0 encounters the hindrance of the weld and subsequently the web during the propagation process, which decreases the propagation rate. In stage 4, the curve form is similar to that of stage 3, except that the trend becomes more obvious. The maximum KI of the entire crack front appears at the characteristic point of 0.96 with values of 696 MPa·mm 1/2 and 870MPa·mm 1/2 , according to the two simulation steps (step 53 and step 56), while the minimum values appear at CP_0 with value of 306 MPa·mm 1/2 and 316MPa·mm 1/2 respectively, which means that that the propagation rate of CP_0 is nearly unchanged. In stage 5, the crack breaks into two parts with CP_0.57 corresponding to the two endpoints of the newly generated crack fronts. So, a jump of KI occurs at the breakpoint, indicating that the propagation rates of the two endpoints on the lower surface of the flange are different. In the beginning of this stage, these two crack fronts still appear nearly elliptical, and their distribution curves of KI exhibit U-shape, which is mainly due to the fact that the propagation rates of the two endpoints both located on the surfaces are greater than those in the flange material. In stage 6, the two cracks evolved into line-shape gradually, resulting in line-shaped distribution curves of KI. It should be noted that the crack on the north side only propagates on the flange, while the crack on the south side propagates on the weld, web, and flange, reflecting the distribution curve of the right part is higher than that of the left part.   Figure 15 shows the distribution of K I of all the six stages with the abscissa of normalized coordinates of the crack front, in which several curves concerning different simulation steps may be of one stage with the same color. It can be seen that in stage 1and 2, the K I of the entire crack front is almost constant and at a low level, which is mainly due to the short crack situation at these two stages. When the crack grows into stage 3, the K I of CP_1 gradually increases, resulting in a K I larger than those of CP_0.57 and CP_0. This is mainly due to the fact that CP_0 encounters the hindrance of the weld and subsequently the web during the propagation process, which decreases the propagation rate. In stage 4, the curve form is similar to that of stage 3, except that the trend becomes more obvious. The maximum K I of the entire crack front appears at the characteristic point of 0.96 with values of 696 MPa·mm 1/2 and 870MPa·mm 1/2 , according to the two simulation steps (step 53 and step 56), while the minimum values appear at CP_0 with value of 306 MPa·mm 1/2 and 316MPa·mm 1/2 respectively, which means that that the propagation rate of CP_0 is nearly unchanged. In stage 5, the crack breaks into two parts with CP_0.57 corresponding to the two endpoints of the newly generated crack fronts. So, a jump of K I occurs at the breakpoint, indicating that the propagation rates of the two endpoints on the lower surface of the flange are different. In the beginning of this stage, these two crack fronts still appear nearly elliptical, and their distribution curves of K I exhibit U-shape, which is mainly due to the fact that the propagation rates of the two endpoints both located on the surfaces are greater than those in the flange material. In stage 6, the two cracks evolved into line-shape gradually, resulting in line-shaped distribution curves of K I . It should be noted that the crack on the north side only propagates on the flange, while the crack on the south side propagates on the weld, web, and flange, reflecting the distribution curve of the right part is higher than that of the left part.  Figure 15 shows the distribution of KI of all the six stages with the abscissa of normalized coordinates of the crack front, in which several curves concerning different simulation steps may be of one stage with the same color. It can be seen that in stage 1and 2, the KI of the entire crack front is almost constant and at a low level, which is mainly due to the short crack situation at these two stages. When the crack grows into stage 3, the KI of CP_1 gradually increases, resulting in a KI larger than those of CP_0.57 and CP_0. This is mainly due to the fact that CP_0 encounters the hindrance of the weld and subsequently the web during the propagation process, which decreases the propagation rate. In stage 4, the curve form is similar to that of stage 3, except that the trend becomes more obvious. The maximum KI of the entire crack front appears at the characteristic point of 0.96 with values of 696 MPa·mm 1/2 and 870MPa·mm 1/2 , according to the two simulation steps (step 53 and step 56), while the minimum values appear at CP_0 with value of 306 MPa·mm 1/2 and 316MPa·mm 1/2 respectively, which means that that the propagation rate of CP_0 is nearly unchanged. In stage 5, the crack breaks into two parts with CP_0.57 corresponding to the two endpoints of the newly generated crack fronts. So, a jump of KI occurs at the breakpoint, indicating that the propagation rates of the two endpoints on the lower surface of the flange are different. In the beginning of this stage, these two crack fronts still appear nearly elliptical, and their distribution curves of KI exhibit U-shape, which is mainly due to the fact that the propagation rates of the two endpoints both located on the surfaces are greater than those in the flange material. In stage 6, the two cracks evolved into line-shape gradually, resulting in line-shaped distribution curves of KI. It should be noted that the crack on the north side only propagates on the flange, while the crack on the south side propagates on the weld, web, and flange, reflecting the distribution curve of the right part is higher than that of the left part.   Figure 16 shows the variation of SIF K I with the crack propagation depth of the characteristic points 0, 0.57, and 1 in the first four stages under the K equ + SERR simulation scheme. In these four stages, K I of the characteristic points 0.57 and 1 increase linearly with the increase of the propagation depth, in which the K I value of CP_0.57 is slightly greater than that of CP_1. The K I of CP_0 increases continuously in stage 1 and stage 2, and a steep drop occurs when it enters the north side weld. After that, it increases again until entering the web with another obvious drop in stage 4. Then, it grows slowly in a nearly constant rate, with a K I value far lower than those of the characteristic points 0.57 and 1. The two drops of the K I value of CP_0 correspond to the fact of being hindered by the weld and the web, respectively. It is interesting that after the hindering, CP_0 will always catch up with the other two characteristic points in order to propagate with almost the same rate.
Metals 2019, 9, x FOR PEER REVIEW 18 of 22 Figure 16 shows the variation of SIF KI with the crack propagation depth of the characteristic points 0, 0.57, and 1 in the first four stages under the Kequ + SERR simulation scheme. In these four stages, KI of the characteristic points 0.57 and 1 increase linearly with the increase of the propagation depth, in which the KI value of CP_0.57 is slightly greater than that of CP_1. The KI of CP_0 increases continuously in stage 1 and stage 2, and a steep drop occurs when it enters the north side weld. After that, it increases again until entering the web with another obvious drop in stage 4. Then, it grows slowly in a nearly constant rate, with a KI value far lower than those of the characteristic points 0.57 and 1. The two drops of the KI value of CP_0 correspond to the fact of being hindered by the weld and the web, respectively. It is interesting that after the hindering, CP_0 will always catch up with the other two characteristic points in order to propagate with almost the same rate. After entering stage 5, two new crack fronts are formed, which are marked as Crack Front 1 and Crack Front 2. Therefore, the characteristic points are marked as CP1_0, CP1_0.57, CP2_0.57, and CP2_1, respectively, among which CP1_0.57 and CP2_0.57 are on the lower surface of the flange. KI curves of the four characteristic points are drawn in Figure 17. It can be obviously seen that in stage 5, the propagation of CP1_0.57 and CP2_0.57, which are on the lower surface of the flange, decrease at first and rise again in stage 6. At the end of stage 6, KI of CP2_0.57 and CP2_1, which are endpoints of the newly generated Crack Front 2, are significantly greater than those of CP1_0 and CP1_0.57, resulting in the rapid instability fracture of the north side flange.  Figure 18 and Figure 19 are the evolution curves of KII and KIII under the Kequ + SERR simulation scheme. Overall, the values of KII and KIII are much smaller than KI. According to Figure 18, KII basically fluctuates between ±20 MPa·mm 1/2 , except a small region near the endpoint (CP_0) in stage 4 (step 53-step 56), with the minimum value of -93 MPa·mm 1/2 . At that situation, CP_0 begins to After entering stage 5, two new crack fronts are formed, which are marked as Crack Front 1 and Crack Front 2. Therefore, the characteristic points are marked as CP1_0, CP1_0.57, CP2_0.57, and CP2_1, respectively, among which CP1_0.57 and CP2_0.57 are on the lower surface of the flange. K I curves of the four characteristic points are drawn in Figure 17. It can be obviously seen that in stage 5, the propagation of CP1_0.57 and CP2_0.57, which are on the lower surface of the flange, decrease at first and rise again in stage 6. At the end of stage 6, K I of CP2_0.57 and CP2_1, which are endpoints of the newly generated Crack Front 2, are significantly greater than those of CP1_0 and CP1_0.57, resulting in the rapid instability fracture of the north side flange.

KII and KIII
Metals 2019, 9, x FOR PEER REVIEW 18 of 22 Figure 16 shows the variation of SIF KI with the crack propagation depth of the characteristic points 0, 0.57, and 1 in the first four stages under the Kequ + SERR simulation scheme. In these four stages, KI of the characteristic points 0.57 and 1 increase linearly with the increase of the propagation depth, in which the KI value of CP_0.57 is slightly greater than that of CP_1. The KI of CP_0 increases continuously in stage 1 and stage 2, and a steep drop occurs when it enters the north side weld. After that, it increases again until entering the web with another obvious drop in stage 4. Then, it grows slowly in a nearly constant rate, with a KI value far lower than those of the characteristic points 0.57 and 1. The two drops of the KI value of CP_0 correspond to the fact of being hindered by the weld and the web, respectively. It is interesting that after the hindering, CP_0 will always catch up with the other two characteristic points in order to propagate with almost the same rate. After entering stage 5, two new crack fronts are formed, which are marked as Crack Front 1 and Crack Front 2. Therefore, the characteristic points are marked as CP1_0, CP1_0.57, CP2_0.57, and CP2_1, respectively, among which CP1_0.57 and CP2_0.57 are on the lower surface of the flange. KI curves of the four characteristic points are drawn in Figure 17. It can be obviously seen that in stage 5, the propagation of CP1_0.57 and CP2_0.57, which are on the lower surface of the flange, decrease at first and rise again in stage 6. At the end of stage 6, KI of CP2_0.57 and CP2_1, which are endpoints of the newly generated Crack Front 2, are significantly greater than those of CP1_0 and CP1_0.57, resulting in the rapid instability fracture of the north side flange.  Figure 18 and Figure 19 are the evolution curves of KII and KIII under the Kequ + SERR simulation scheme. Overall, the values of KII and KIII are much smaller than KI. According to Figure 18, KII basically fluctuates between ±20 MPa·mm 1/2 , except a small region near the endpoint (CP_0) in stage 4 (step 53-step 56), with the minimum value of -93 MPa·mm 1/2 . At that situation, CP_0 begins to  Figures 18 and 19 are the evolution curves of K II and K III under the K equ + SERR simulation scheme. Overall, the values of K II and K III are much smaller than K I . According to Figure 18, K II basically fluctuates between ±20 MPa·mm 1/2 , except a small region near the endpoint (CP_0) in stage 4 (step 53-step 56), with the minimum value of -93 MPa·mm 1/2 . At that situation, CP_0 begins to propagate into the surface of the web, which seems to have some special influences on K II . According to the Figure 19, the variation of K III is more obvious. In the first three stages, K III fluctuates between ±20 MPa·mm 1/2 . From stage 4 on, the variation of K III increases significantly with a maximum value of no more than 100 MPa·mm 1/2 . However, although it seems that K II and K III varies greatly, their magnitudes are far lower than K I , which again demonstrates that K I is the dominant driving force controlling the crack propagation of corrugated web girders. propagate into the surface of the web, which seems to have some special influences on KII. According to the Figure 19, the variation of KIII is more obvious. In the first three stages, KIII fluctuates between ±20 MPa·mm 1/2 . From stage 4 on, the variation of KIII increases significantly with a maximum value of no more than 100 MPa·mm 1/2 . However, although it seems that KII and KIII varies greatly, their magnitudes are far lower than KI, which again demonstrates that KI is the dominant driving force controlling the crack propagation of corrugated web girders.

Fatigue life
Based on Nasgro law, fatigue life, life ratio, and the propagation depth of every stage are calculated under KI + MTS, KI + SERR, and Kequ + SERR schemes, as shown in Tables 5, 6, and 7. It can be seen that the fatigue life under these three schemes are similar, especially the total life which are all more than 120 × 10 4 cycles and the errors with the test life all less than 7%. So, the calculation accuracy of this method is much higher than that of the ESS method. Moreover, the calculation life is less than the test life, which indicates that it is even safer in engineering applications. Table 6 lists the life ratio of every stage. It can be found that the life of stage 1 accounts for only a small part of the whole life, with the value from 2.7% to 3.2% under three schemes, indicating that the crack kink behavior lasts a very short time. Compared with the other stages, the life of stage 3 accounts for the vast majority of the entire life with the value of more than 100 × 10 4 cycles. Stage 2 covers the second longer life with about 12 × 10 4 cycles and stage 4 covers the third longest life with 7-8 × 10 4 cycles. It is worth noting that the crack propagates in the depth direction during stages 2, 3, and 4 with their life sums under three simulation schemes all exceeding 95% of the total life. Table 7 lists the propagation depth of the crack in every stage. It can be seen that the propagation depths of stage 1 and 2 are extremely small, with both on the order of 0.1 mm. The propagation depths of stage 3 and 4 are comparatively larger-about 11 mm in stage 3 and about 8 mm in stage 4-with the sum of these two stages of about 19 mm. Considering that the thickness of the lower flange is 20 mm, so the propagation ratio of these two stages is almost 95%. It should be noted that although the propagate into the surface of the web, which seems to have some special influences on KII. According to the Figure 19, the variation of KIII is more obvious. In the first three stages, KIII fluctuates between ±20 MPa·mm 1/2 . From stage 4 on, the variation of KIII increases significantly with a maximum value of no more than 100 MPa·mm 1/2 . However, although it seems that KII and KIII varies greatly, their magnitudes are far lower than KI, which again demonstrates that KI is the dominant driving force controlling the crack propagation of corrugated web girders.

Fatigue life
Based on Nasgro law, fatigue life, life ratio, and the propagation depth of every stage are calculated under KI + MTS, KI + SERR, and Kequ + SERR schemes, as shown in Tables 5, 6, and 7. It can be seen that the fatigue life under these three schemes are similar, especially the total life which are all more than 120 × 10 4 cycles and the errors with the test life all less than 7%. So, the calculation accuracy of this method is much higher than that of the ESS method. Moreover, the calculation life is less than the test life, which indicates that it is even safer in engineering applications. Table 6 lists the life ratio of every stage. It can be found that the life of stage 1 accounts for only a small part of the whole life, with the value from 2.7% to 3.2% under three schemes, indicating that the crack kink behavior lasts a very short time. Compared with the other stages, the life of stage 3 accounts for the vast majority of the entire life with the value of more than 100 × 10 4 cycles. Stage 2 covers the second longer life with about 12 × 10 4 cycles and stage 4 covers the third longest life with 7-8 × 10 4 cycles. It is worth noting that the crack propagates in the depth direction during stages 2, 3, and 4 with their life sums under three simulation schemes all exceeding 95% of the total life. Table 7 lists the propagation depth of the crack in every stage. It can be seen that the propagation depths of stage 1 and 2 are extremely small, with both on the order of 0.1 mm. The propagation depths of stage 3 and 4 are comparatively larger-about 11 mm in stage 3 and about 8 mm in stage 4-with the sum of these two stages of about 19 mm. Considering that the thickness of the lower flange is 20 mm, so the propagation ratio of these two stages is almost 95%. It should be noted that although the

Fatigue Life
Based on Nasgro law, fatigue life, life ratio, and the propagation depth of every stage are calculated under K I + MTS, K I + SERR, and K equ + SERR schemes, as shown in Table 5, Table 6, and Table 7. It can be seen that the fatigue life under these three schemes are similar, especially the total life which are all more than 120 × 10 4 cycles and the errors with the test life all less than 7%. So, the calculation accuracy of this method is much higher than that of the ESS method. Moreover, the calculation life is less than the test life, which indicates that it is even safer in engineering applications.   Table 6 lists the life ratio of every stage. It can be found that the life of stage 1 accounts for only a small part of the whole life, with the value from 2.7% to 3.2% under three schemes, indicating that the crack kink behavior lasts a very short time. Compared with the other stages, the life of stage 3 accounts for the vast majority of the entire life with the value of more than 100 × 10 4 cycles. Stage 2 covers the second longer life with about 12 × 10 4 cycles and stage 4 covers the third longest life with 7-8 × 10 4 cycles. It is worth noting that the crack propagates in the depth direction during stages 2, 3, and 4 with their life sums under three simulation schemes all exceeding 95% of the total life. Table 7 lists the propagation depth of the crack in every stage. It can be seen that the propagation depths of stage 1 and 2 are extremely small, with both on the order of 0.1 mm. The propagation depths of stage 3 and 4 are comparatively larger-about 11 mm in stage 3 and about 8 mm in stage 4-with the sum of these two stages of about 19 mm. Considering that the thickness of the lower flange is 20 mm, so the propagation ratio of these two stages is almost 95%. It should be noted that although the propagation depth of stage 4 is larger than that of stage 2, but its life is only nearly 10% of stage 3, and even less than the life of stage 2.
According to the analysis above, it can be suggested that the propagating life of an initial semi-circle crack with depth of 0.1 mm propagating through the thickness of the flange, which covers stages 2, 3, and 4, can be utilized to predict the life of corrugated web girders under three-point bending with sufficient accuracy.

Conclusions
In this paper, the entire process of the weld toe surface crack of a corrugated web girder specimen under cyclic three-point bending was simulated for different crack propagation driving forces and crack kink criteria. The evolution of crack front shape, kink angle, K I , K II , K III , and propagation life in different stages were analyzed. The main conclusions are: (1) Under the cyclic bending load, the fatigue crack of the corrugated web girder appears mainly in the weld toe of the flange side, which is at the intersection of the inclined section and the bending radius section in the middle part of the girder span. The ESS method can accurately capture the initiation location of the weld toe crack, but the predicted fatigue life is unsatisfying. (2) A semi-circle initial crack with a = 0.1 mm was inserted at the crack initiation location and the evolution processes under different schemes were simulated. The simulation results show that the MTS and SERR kink criteria are consistent with the test results. Therefore, it is considered that the crack in this investigation is of opening mode, in which K I is the dominating propagation driving force. (3) In terms of the crack front shape and the location of the end points, the entire crack propagation process can be divided into six stages. In stage 1, the crack kinked and the kink angle was related to the inclined angle of the corrugation profile and the principal stress direction of the lower flange. In all the following stages, the crack propagated approximately in the same plane.
(4) In a very short time, the crack with semi-circle shape quickly evolved into a semi-ellipse crack in stage 1 and finally evolved into two line-shape cracks in the last two stages. The weld and the web are the geometry factors influencing the crack shape and its propagation behavior. (5) SIF of three characteristic points at the crack front are analyzed. K I is the key driving force and its distribution directly affects the propagation rates of characteristic points. Compared with K I , the values of K II and K III are lower, and thus can be ignored. Funding: This research was funded by National Natural Science Foundation of China (51575408).

Conflicts of Interest:
The authors declare no conflict of interest.