Experimental and Numerical Investigations on Crack Intersection and Propagation of Concrete Structures

This study presents experimental and numerical methods to reveal concrete structures’ crack propagation and intersection laws under static and dynamic loads. Firstly, a numerical simulation method was established using the user‑defined material subroutines to solve the free crack surface contact problem of concrete structures. Secondly, three‑point bending tests of concrete beams containing double cracks of different approaching angles were car‑ ried out based on the digital image correlation (DIC) technology to study the intersection and propagation of cracks. Finally, the fracture processes of double‑crack concrete beams and concrete gravity dams with and without a longitu‑ dinal crack were simulated and analyzed under static and dynamic loads. The numerical results were compared with the test results to verify the effectiveness and accuracy of the proposed method in simulating crack intersection and propagation in concrete structures. Results indicate that the crack intersection affects the fracture path of concrete structures, weakens their bearing capacity, and accelerates the failure of the structures. The proposed simulation method provides an effective technical approach for crack propagation prediction and safety evaluation of engineer‑ ing structures.


Introduction
Various cracks, even large-scale longitudinal penetrating cracks (Batta & Pekau, 1996;Horii & Chen, 2003;Zhang et al., 2013), occur in concrete structures under temperature change, foundation settlement, or overload due to their low material tensile strength. In addition, the gravity dams with longitudinal cracks built in earlier years form open longitudinal cracks due to poor grouting or no grouting treatment (Cui et al., 2011). Under the external environment or load, the original cracks in the structure expand and intersect, while the new cracks may also intersect with the original cracks (Shi et al., 2013). Numerous studies have shown that the expansion and intersection of these cracks have an essential impact on the bearing capacity and durability of concrete structures and are the leading causes of engineering structures' ultimate failure (Wen et al., 2020;Zhang et al., 2013). Single crack growth processes in concrete under static and dynamic loadings have been systematically studied (Bhosale et al., 2020;Fan et al., 2020;Wu et al., 2011;Xu & Reinhardt, 1999a, 1999b. However, similar researches on cracks' expansion and interaction (including intersection) are still in its initial stage. The bearing capacity and fracture performance of the concrete structures affected by the interaction between cracks differ from those determined by the single crack model (Ghahremannejad et al., 2018;Shi et al., 2003Shi et al., , 2013. The traditional durability and safety performance evaluation methods that only consider the expansion of a single crack can no longer ensure the safety of the structure (Zhan et al., 2020). Therefore, the studies on the interaction and expansion of cracks and their research methods are of great significance for structural hazard predictions and safety assessments.
In recent years, many researchers have combined experiments and numerical simulations to study the effect of crack propagation and interaction on fracture failure characteristics of concrete. Wang et al. (2018) conducted semi-circular bending (SCB) tests and numerical simulations to study the intersection of the pre-set and cemented natural cracks with non-vertical approaching angles. Wen et al. (2020) proposed a new multi-crack length control method to analyze the propagations of multiple cracks in concrete structures. Shi et al. (2004) conducted a numerical simulation analysis on the fracture mode of plain concrete beams with multiple cracks, focusing on the localization degree and interaction between cracks. Pimanmas and Maekawa (2001) studied the influence of crack intersection on the shear performance of reinforced concrete beams, finding that crack intersection has a predominant impact on the shear bearing capacity, displacement-load curve, failure characteristics, and other fracture characteristics of beams. Bhattacharjee and Léger (1993) and Calayir and Karaton (2005) used the rotating crack model to analyze the dynamic cracking of the Koyna gravity dam under considering the interaction of non-orthogonal cracks. However, previous studies mainly focused on the propagation and interaction of cracks without crossing (such as parallel cracks). The intersection and connection of cracks can cause even faster or more severe damage.
In test methods, the digital image correlation (DIC) technology has become a critical non-contact measurement method to study the crack propagation process and interaction in concrete structures under static and dynamic loads on a laboratory scale (Bu et al., 2020;Fayyad & Lees, 2017;Li et al., 2018Li et al., , 2020 due to its ability to more intuitively analyze the development of micro-cracks and the whole fracture process. In numerical simulation methods (Liu et al., 2020a(Liu et al., , 2020b, the cohesive crack model (CCM) is increasingly popular in applying research on crack propagation and interaction phenomena. The CCM simulates the reduced cohesive force between crack surfaces by describing the tractionseparation relationship between the element surfaces (Choi & Park, 2019;Li et al., 2021). Many researchers combined the CCM with the finite element method FEM by pre-inserting cohesive elements to simulate concrete's complex crack propagation and interaction under static and dynamic loads (Camanho et al., 2003;Maimí et al., 2007aMaimí et al., , 2007bYang et al., 2009). Su et al. (2010) simulated complex crack propagations in concrete wedge splitting tensile and torsional fracture models under static and dynamic loads by batch inserting cohesive elements into initial finite elements. Xu et al. (2014) simulated and analyzed the failure process of the two-dimensional model of the Koyna gravity dam under earthquake load by embedding the cohesive elements in the solid elements. Zhang et al. (2016) established a series of three-point bending fracture models of concrete beams with parallel double cracks by inserting the cohesive elements into the mortar and interface zone as the potential fracture area. Yang and Chen (2005) proposed a fully automatic finite element model to simulate the propagations of multiple interactive cracks in reinforced concrete beams based on the CCM.
In the crack intersection and propagation processes under complex loads, especially the constantly changing dynamic loads, closed-free crack surfaces are very likely to appear. When analyzing the closed-free crack, the contact problem between the crack surfaces is inevitable. That is, it must be considered that the crack surfaces are not allowed to invade each other in the actual contact process, and the normal component of the contact force is not allowed to be tension but only pressure. In addition, the tangential friction behavior between the crack surfaces is also an essential factor affecting the mechanical properties of the crack itself and the structural stability (Wang et al., 2010). Therefore, the contact characteristic of free crack surfaces is a critical factor in crack propagation and intersection simulating. However, the contact problem in the crack propagation process is highly nonlinear; the contact process is usually time-dependent, accompanied by the evolution of material nonlinearity and geometric nonlinearity. At the same time, the contact interface's area, shape, kinematic, and mechanical state are also unknown in advance. In addition, it is difficult to introduce the contact surface conditions into the solution process and track the contact pair's relative position, contact state, and area in finite element software without some particular elements or subroutines. In the iterative calculation process, the contact state may constantly change between adhesion, sliding, and separation, which will cause iterative divergence and failure of calculation, especially when the crack surface presents a sharp serration or under the action of dynamic load. These characteristics determine the complexity and difficulty of the contact problem, so few scholars consider this factor.
This study presents experimental and numerical research methods of crack intersection and propagation in concrete structures under external static and dynamic loads. The study is organized as follows. In Sect. 2, the numerical simulation method is given in detail. Section 3 presents experiments and numerical simulations for Page 3 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 crack intersection and propagation of double-crack concrete beams. In Sect. 4, crack intersection and propagation in the gravity dam under static and dynamic loads are simulated and analyzed. Finally, conclusions are included in Sect. 5.

Subroutine Compilation
This study compiled the user-defined material subroutines (UMAT and VUMAT) to solve the complex contact problem of free crack surfaces in the crack propagation and intersection processes under static and dynamic loads. The material subroutines defined the material constitutive models of cohesive elements before and after forming free crack according to the bilinear traction-separation criterion in CCM and the contact behavior of free crack surfaces, respectively. In a two-dimensional case, the bilinear traction-separation criterion assumes the following: the normal and tangential traction forces ( t n and t s ) increase linearly with the crack opening and sliding displacements ( δ n and δ s ), respectively, before the crack initiation and decrease linearly after the damage begins. Fig. 1 shows the bilinear traction-separation criterion (Yang et al., 2009), where K 0 n and K 0 s represent the initial normal and shear stiffnesses when the material is not cracked. The values of K 0 n and K 0 s should be high enough to describe the uncracked material but not too high to cause numerical illconditioning problems, which can be obtained by referring to Eq. (1) (Qiang et al., 2000): where v, b, and E are Poisson's ratio, the characteristic length, and Young's modulus of the initial elements, respectively, and c is the amplification factor, usually taken as 10 ~ 100.
For I-II mixed-mode crack propagation in the simulation, the quadratic nominal stress criterion was used as the initial critical damage criterion. When the initial thickness of the cohesive element was defined as 1, the criterion can be expressed by where t 0 n and t 0 s represent the maximum nominal stresses of pure mode I and mode II initial damages, respectively; ‹› is the Macaulay bracket, and its meaning is as follows: To describe the damage evolution under I-II mixedmode crack propagation, the effective displacement δ m (Maimí et al., 2007a(Maimí et al., , 2007b was introduced, which can be determined by the following equation: (1) (3) �t n � = t n , t n ≥ 0 0, t n < 0 . Fig. 1 The bilinear traction-separation model: (a) normal constitutive model and (b) tangential constitutive model (Maimí et al., 2007a) Page 4 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 and the degeneration processes of mechanical properties such as stiffness and stress were described by damage variable D, which is defined as where δ max m is the maximum effective displacement in the loading process so far, δ 0 m is the effective displacement at initial damage, and δ f m is the effective displacement under complete failure ( D = 1 ), which can be calculated according to the Benzeggagh-Kenane rule (Benzeggagh & Kenane, 1996) and shown in Eqs. (6) to (11): (4) δ m = �δ n � 2 + δ 2 s , where G n and G s , δ n and δ s , t n and t s are the energy release rates, displacements, and stresses of normal and tangential directions when the material is destroyed; G C n and G C s are mode I and mode II fracture's critical energy release rates; δ f n and δ f s represent the critical displacements when the cohesive element is completely pulled apart and sliding cracking, respectively.
If δ n is negative during iterations, a penalty stiffness k n is introduced to prevent the crack surfaces from kinging into each other. During damage evolution, t n and t s can be obtained from the following equation: where t n and t s are the stress components obtained according to the traction-separation criterion when the material is not damaged.
After the formation of free crack surfaces, the penalty stiffness k n must also be introduced to ensure that the contact surfaces are not embedded. The penalty function friction model with fairly good convergence was used to simulate the friction behavior between free crack surfaces. As shown in Fig. 2a and b, t n and t s of the free crack surfaces can be calculated by Eqs. (14) and (15): where κ is the slope (penalty stiffness) of "elastic slip", δ lit is the critical slipping displacement of "elastic slip", μ is the friction coefficient, and t is the unit vector consistent with the slipping direction. The user subroutine interface of ABAQUS was utilized to assign the defined material properties to the cohesive elements to realize the crack propagation and intersection behavior simulation in the concrete structure under static and dynamic loads.
In VUMAT designed in this paper, owing to the consideration of the stresses transmitted between the free crack surfaces, the cohesive elements between the crack interfaces were not deleted after the formation of the real free crack. The corresponding state variable SDV5 was set in VUMAT as the damage mark value of the element to observe the failure element and fracture path during post-processing, which is defined as The substances and calling flow of UMAT or VUMAT are shown in Fig. 3.

Specimen Preparation and Material Properties
Four groups of concrete beams with pre-set double cracks at different approaching angles θ were designed to reveal the influence of crack intersection on the fracture performance of concrete and verify the effectiveness of the proposed numerical simulation method. Threepoint bending fracture tests of the concrete beams with double cracks and concrete beams with the only edge pre-set crack were conducted under quasi-static and seismic dynamic displacement loadings of 10 -3 mm/s and 1 mm/s. The specially made wooden formwork with grooves was used to embed 2-mm-thick steel sheets before pouring to obtain prefabricated cracks with different θ. Fig. 4a and b shows the pouring formwork and pouring condition, respectively. The specific geometric dimensions of the specimens and the design of the preset crack parameters are shown in Fig. 5 and Table 1. The center of the internal penetrating crack is on the vertical extension line of the edge crack. In this test, the length L, height H, thickness B, and the span S between supports of all specimens are 500 mm, 100 mm, 100 mm, and 400 mm, respectively. According to the θ and loading rate, the specimens were named TPB-Aj-Tp, j represents the value of θ, and p is 1 or 2, representing the loading rate of 10 −3 mm/s and 1 mm/s, respectively.
The cement, coarse aggregate, and fine aggregate used in the test are Portland cement (PO 42.5), gravel limestone with a maximum particle size of 20 mm, and natural river sand, respectively. The mix proportion of the concrete is cement:sand:stone:water = 1:1.37:2.78:0.47. Table 2 lists Young's modulus E, compressive strength f c , tensile strength f t , Poisson's ratio v, and fracture energy G F of concrete cured for 28 days.

Experimental Set-up and Instrumentation
As shown in Fig. 6, the test device consists of the MTS-810NEW hydraulic servo dynamic testing machine and DIC set-up. The MTS testing machine can directly obtain the specimen's mid-span vertical displacement and loading rate. The crack mouth opening displacement (CMOD) of the edge mid-span crack of the specimen is obtained by connecting the large clip extensometer in the acquisition system of the MTS testing machine and the synchronous acquisition of load and CMOD is realized. The three-dimensional DIC technology (Liu et al., 2020a(Liu et al., , 2020bYan & Lin, 2006;Yi-Qiu et al., 2012) is introduced to intuitively monitor cracks' intersection behavior. To accurately track the deformation of the specimen surface, a random matte black-and-white spot pattern must be sprayed on the smooth surface to ensure that the analyzed area has distinct characteristics. The speckled area covers the crack propagation area (approximately 100 × 200 mm in the mid-span) shown in Fig. 7.
Commercial PMLAB software is used for image acquisition and analysis in this test. Two cameras with 2024 × 2024 pixels and 50 mm fixed focal length, and 80 Hz maximum acquisition frequency are used to shoot the specimen's mid-span area of about 0.25 × 0.25 m. Before the test, a standard 6 mm calibration plate is used to calibrate the 3D-DIC setting (the calibration error is 0.02 pixel). The first group of images is collected as the reference images of the no-load undeformed specimen. During the test, digital images of two loading rates are collected at the acquisition frequencies of 1 Hz and 10 Hz, respectively, until the samples failed. The fracture tests of each group are repeated three times.

Finite Element Modeling
A series of two-dimensional beam models with the same geometric dimensions and pre-set crack parameters as the test beams in Sect. 3.1 were established in ABAQUS/ CAE. To ensure the calculation accuracy and efficiency, a focus area with a width of h was set in the mid-span area of the beam based on the distribution of fracture paths observed in the test. When θ is 30°, 45°, 60°, or 90°, h is 40 mm, 50 mm, 50 mm, or 70 mm, respectively.
The model TPB-A30-T1 taken as an example to show the element divisions and boundary conditions is shown in Fig. 8. Both sides of the focus area were discretized by 4-node plane stress (CPS4) elements, and the element length transited from 4 to 10 mm from the mid-span to both ends of the beam. The focus area was discretized by the refined triangular plane stress (CPS3) elements; among them, the element lengths of the edge crack Page 6 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 expansion area and the intersection area with the internal penetrating crack are about 1 mm. Then, the 4-node zero thickness cohesive (COH2D4) elements (Fig. 8c) were inserted into the initial meshes of the focus area (Fig. 8b) in batches by using the algorithm. The bottom supports were fixed, and the tangential friction coefficient between supports and beam was set as 0.3, i.e., the sliding friction coefficient between concrete and steel. The vertical displacement loads with loading rates of 10 -3 mm/s and 1 mm/s consistent with the test conditions were applied on the beam top's mid-span. The material parameters E and ν of the initial elements and the normal critical fracture energy G C n and critical tensile strength t 0 n of the cohesive elements can be obtained from the average values of E, ν, G F , and f t measured by experiments in Sect. 3.1. The tangential  Wang et al. Int J Concr Struct Mater (2022) 16:67 material parameters of the cohesive elements are difficult to be measured through the test and have large discreteness. Therefore, as with the initial stiffnesses K 0 n and K 0 s , they were adjusted by the experimental P-CMOD curves using the inverse analysis method (Hanson et al., 2004;Zhang et al., 2016). K 0 n and K 0 s adjustment ranges also need to refer to Eq. (1). The material parameters of the cohesive elements are listed in Table 3. ABAQUS/Explicit module was used to simulate the whole fracture processes of specimens.

Analysis of Test and Simulation Results
The whole propagation and intersection processes of non-vertical approaching angle cracks and vertical approaching angle cracks in concrete beams are presented by taking models TPB-A45-T2 and TPB-A90-T2 as examples. The principal strain (exx) fields in the measurement areas of different fracture stages calculated by the DIC method are shown in Figs. 9 and 11. The maximum principal stress distribution nephograms of the specimens obtained by numerical simulation are shown in Figs. 10 and 12. Through comparison, it can be concluded that the simulated fracture propagation processes    -T1  20  1  40  30  3   TPB-A30-T2  20  1  40  30  3   TPB-A45-T1  20  1  40  45  3   TPB-A45-T2  20  1  40  45  3   TPB-A60-T1  20  1  40  60  3   TPB-A60-T2  20  1  40  60  3   TPB-A90-T1  20  1  40  90  3   TPB-A90-T2  20  1  40  90  3 Page 8 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 of concrete beams under crack intersection are generally consistent with the experimental results. For the model TPB-A45-T2, the fracture processes can be divided into three stages. As shown in Figs. 9a and 10a, the stress concentrations mainly occurred in the edge pre-set crack tip area and the lower crack tip area of the internal pre-set crack in stage I. With the continuous loading, a micro-crack appeared in the lower crack tip area of the internal pre-set crack. The propagation direction of the edge crack gradually deviated from the original direction to the direction perpendicular to the internal pre-set crack. In stage II, the edge crack propagated near the internal pre-set crack, and the stress concentration and micro-crack appeared in the upper crack tip area of the internal crack, as shown in Figs. 9b and 10b. In stage III, the micro-crack at the upper crack tip of the internal crack gradually expanded towards the loading point. The micro-crack at the lower crack tip of the internal crack gradually disappeared, as shown in Figs. 9c and 10c. Similar crack propagation patterns can be observed in all concrete beams containing non-vertical approaching angle double cracks.
For the model TPB-A90-T2, the fracture processes can also be divided into three stages. The differences with the concrete beams containing non-vertical approaching angle double cracks are as follows: in stage I, the microcrack only occurred in the edge pre-set crack tip area and expanded forward with the loading until it meets the internal transverse crack, as shown in Figs. 11a and 12a; in stage II, when the crack extended to the transverse pre-set crack, stress concentrations and micro-cracks appeared in the crack tip and the middle area of the upper surface of the transverse pre-set crack, as shown in Figs. 11b and 12b; in stage III, the micro-crack in the middle of the upper surface gradually increased and expanded towards the loading point, while the microcrack at the crack tip gradually disappeared due to stress release, as seen in Figs. 11c and 12c. Figs. 13,14,15,16 show the principal strain and SDV5 nephograms of concrete beams containing double cracks with different θ when the loading rates are 10 -3 mm/s and 1 mm/s, respectively. Through comparison, it can be concluded that the fracture paths obtained by numerical simulation are consistent with the test results (due to the randomness of concrete aggregate, the small error can be ignored). A larger θ and a more significant loading rate result in a fewer deviation of the edge crack propagation path. Thus, the accuracy of the numerical simulation method proposed above to predict the concrete fracture processes under the intersection of double cracks at different θ is proved, and the phenomenon that the crack intersection affects the fracture mode is observed. Figs. 17 and 18 show the P-CMOD curves of specimens with different θ at two loading rates obtained by numerical simulations and tests. Through comparison, it is found that the simulation results are consistent with the test results, except that the slopes of the linear elastic stage and descent section are slightly different. In addition, the peak load P max and the critical opening displacement CMOD c corresponding to P max obtained by experiments and simulations are very close. It can also be observed that the numerical simulation method is more likely to observe the fluctuation of the curves caused by the existence of the internal crack.
The variation trends of P max with the increase of θ under the two loading rates obtained by experiments and simulations are plotted in Fig. 19. When θ increases from 30° to 45°, 60°, and 90°, the averages of P max obtained from the test are reduced by 58.7%, 39.4%, 16.6%, and 4.1%, respectively, compared with P max (3.68 kN) of the single crack beam at the loading rate of 10 -3 mm/s; the averages of P max are reduced by 59.3%, 33.0%, 25.3%, and 8.8%, respectively, compared with P max (5.46 kN) of the single crack beam at the loading rate of 1 mm/s. The  Page 10 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 results demonstrate that the intersection of cracks affects the bearing capacity of concrete beams. The smaller θ leads to the more obvious weakening of the bearing capacity of the beams. In addition, it can be observed that P max at the loading rate of 1 mm/s is greater than that of the loading rate of 10 -3 mm/s. To demonstrate the necessity of considering the normal and tangential contact characteristics of free crack surfaces in crack propagation modeling, the numerical results of the model TPB-A30-T1 are given as an example of using the material constitutive model that does not consider crack surface contact and friction. As shown in Fig. 20, the contact surfaces of the closed-free crack are embedded with each other due to the neglect of the normal contact action of free crack surfaces, which is inconsistent with the actual situation. As shown in Fig. 21, TPB-A30-T1Y and TPB-A30-T1N represent P-CMOD curves obtained from numerical simulation results that consider and ignore free crack contact properties, respectively. It can be seen that the P max of TPB-A30-T1N is less than that of TPB-A30-T1Y because the normal contact and tangential friction of the free crack surfaces are not considered.

Establishment of Models
The simulations of the cracking process and crack intersection in dams under earthquake are of great significance for analyzing the dynamic security of dams under earthquake. Therefore, the Koyna gravity dam commonly Page 13 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 used in dynamic analysis and the Koyna gravity dam with a penetrating longitudinal crack were simulated to study the cracking and crack intersection processes under seismic loading. The penetrating longitudinal crack was set according to the construction requirements and the actual fracture path of the dam. Moreover, the impact of crack intersection on the safety performance of the gravity dam was analyzed through comparison. The Koyna dam was built on the Koyna River in 1963, located in Maharashtra state, 240 km southeast of Mumbai, India. The dam section with the largest dam height (103 m) was considered the research object. According to the records, the water level in front of the reservoir was 91.75 m when the dam suffered a strong earthquake. The model dimensions of the Koyna gravity dams with and without the longitudinal crack are shown in Fig. 22a and b, respectively. The intact Koyna dam and the Koyna dam with a longitudinal crack were modeled in ABAQUS/CAE, as shown in Fig. 23. The areas near the slope break of the downstream dam surface were divided by the densified plane strain triangular (CPE3) elements. The parts of the upstream dam surface and other areas were divided by the plane strain quadrilateral reduced integral (CPE4R) elements. To improve the computational efficiency, the cohesive (COH2D4) elements were only inserted in batches between the initial meshes of the regions of interest to simulate potential cracks. The longitudinal crack in Page 14 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 the dam was simulated by inserting cohesive elements. The contact properties between longitudinal crack surfaces were set by defining the constitutive relationship of cohesive elements. The measured time history curves of horizontal and vertical seismic wave acceleration with peak values of 0.474 g and 0.312 g shown in Fig. 24a and b, respectively, were input at the bottom of the dams to simulate the seismic loads (unit: g = 9.81 m/s 2 ). At the same time, water pressures were applied to the upstream dam surface.
The E and v of the dam's initial elements were considered 31,027 MPa and 0.2, and the densities of concrete and water were set as 2463 kg/m 3 and 1000 kg/m 3 , respectively. The Rayleigh damping coefficient β was set as 3.23 × 10 -3 (Sheibany & Ghaemian, 2006). The material properties of the cohesive elements were specified through the compiled UMAT and VUMAT subroutines, where t 0 n and t 0 s were 2.9 MPa, G C n and G C s were 250 N/m, μ was 0.6, and the initial stiffnesses of the cohesive elements K 0 n and K 0 s can be estimated by Eq.
(1), considered as 3.5 × 10 3 MPa/mm in the calculation models in this section. For the Koyna gravity dam without the longitudinal crack: at t = 2.7 s, the crack first sprouted in the area where the slope of the downstream dam surface changes sharply and then expanded obliquely to the upstream of the dam; at t = 2.74 s, the opened real crack closed again, and reopened at 2.9 s; at t = 3.17 s, the crack finally extended to the upstream dam surface, and transverse penetration occurred. As the stresses at the crack tip were caused by the bending and the self-weight of the dam, the crack extended downward bending rather than horizontally upstream, entirely consistent with the actual earthquake damage process and the model test result (National Research Council, 1990) shown in Fig. 26, thus illustrating that the numerical simulation method proposed in this study can effectively simulate dam failure mode under strong earthquake.
For the Koyna dam with the longitudinal crack: at t = 2 s, a crack appeared at the lower end of the longitudinal crack and developed downstream; at t = 2.36 s, a crack appeared near the lower end of the longitudinal crack and developed upstream; at t = 2.4 s, a crack appeared at the upper end of the longitudinal crack and developed upstream. With the changes in seismic load value and direction, the real cracks that have been Page 16 of 21 Wang et al. Int J Concr Struct Mater (2022)   Page 17 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 Fig. 28a and b shows the horizontal and vertical displacement responses of the node at the dams' top without and with the longitudinal crack under earthquake, respectively. It can be seen that the time intervals of displacement peak points of the dam with the longitudinal crack are prolonged, which reflects that the longitudinal crack reduces the overall stiffness of the dam and leads to the advance of crack initiation time. The two vertical displacement time history curves are similar before 2 s (i.e., the time of crack initiation in the dam with the longitudinal crack); they are different but have the same trend after this time. Moreover, for the dam with the longitudinal crack, the time of the entire displacement time history, i.e., the time of complete fracture, is less than that without the longitudinal crack.

Conclusions
Experimental and numerical simulation studies were conducted on the crack intersection and propagation processes for double-crack concrete beams and gravity dams under static and dynamic loads. The numerical results were compared with the test results to verify the effectiveness and accuracy of the proposed method for simulating crack intersection and propagation in concrete structures. The proposed numerical simulation method provides effective technical means for crack propagation prediction and safety evaluation of engineering structures. The main conclusions are concluded as follows: Page 18 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 1. The material constitutive model considering free crack surface contact properties was established using the user-defined material subroutines UMAT or VUMAT. The contact simulation problem of the crack intersection process under complex static and dynamic loads was accurately solved by giving the constitutive model to the cohesive elements inserted between the initial finite elements. 2. The effects of crack approaching angle θ on P-CMOD relation, crack propagation processes, and loadcarrying capacity of double-crack concrete beams were extensively investigated. The crack intersection leads to the change in the fracture processes and the weakening of bearing capacity: the P max of beams and the deviation between the propagation direction of the edge crack and the direction perpendicular to the internal pre-set crack decrease with the decreasing θ. The averages of P max under the loading rate of 1 mm/s are larger than that under the loading rate of 10 -3 mm/s, but both are smaller than the single crack Page 19 of 21 Wang et al. Int J Concr Struct Mater (2022) 16:67 concrete beams without crack intersection under the corresponding loading rate. 3. The major crack in the dam expands upstream in repeated opening and closing with the change of seismic load, and finally, a downward bending transverse crack is formed through the dam. Compared with the intact dam, many micro-cracks have been formed at the tip of the longitudinal crack before forming the major crack. The existence of the longitudinal crack and the intersection of cracks change the failure path and reduce the safety performance of concrete dam structures under static water pressures and dynamic seismic loads.