A Numerical Method for Applying Cohesive Stress on Fracture Process Zone in Concrete Using Nonlinear Spring Element

Abstract

Aiming at the numerical simulation of the entire crack propagation process in concrete, a numerical method is proposed, in which cohesive stress on the fracture process zone (FPZ) is simulated and applied by a nonlinear spring element. Using displacement control, the cohesive stress values on the FPZ are obtained from solving a system of nonlinear equations through an iterative process. According to a crack propagation criterion based on initial fracture toughness, the approach adds the spring elements to finite element analysis when simulating mode I crack propagation in standard three-point bending notched concrete beams with different strengths, initial crack ratios ${(\mathrm{a}}_0/\mathrm{D})$, and depths (D). The simulated load versus displacement (P-Delta) curves are performed to recalculate the fracture energy and verify the accuracy of cohesion in the proposed method. The simulated load versus crack mouth opening displacement (P-CMOD) curves are consistent with the previous experimental results. Subsequently, the variations of the FPZ length and the crack extension resistance (${\mathrm{K}}_{\mathrm{R}})$ curves are studied according to the proposed iterative approach. Compared with the existing methods using a noniterative process, the iterative approach generates a larger maximum FPZ length and ${\mathrm{K}}_{\mathrm{R}}$ curve where the FPZ length is mainly determined by the fracture energy, tensile strength, and geometry shape of the beam, and the ${\mathrm{K}}_{\mathrm{R}}$ curve is primarily determined by the fracture energy and FPZ length. The significant differences in numerical results indicate that the applying cohesion is essential in numerical simulation. It is reasonable to conclude that the proposed nonlinear spring element is more applicable and practical in the numerical simulation of the concrete mode I crack propagation process by improving the accuracy of the cohesion applied on the FPZ.

1. Introduction

It is well known that concrete is a quasi-brittle material. The cracking behavior of notched concrete beams is influenced by size effects and environmental conditions. Fracture tests have been conducted to investigate the effects of these factors [1,2]. In addition, scholars have proposed different fracture criterions to determine the instability of cracks in concrete, such as the maximum circumferential stress criterion [3], the maximum energy release rate criterion [4], and the double-K fracture criterion [5,6]. Specifically, the double-K fracture criterion uses crack initial fracture toughness and unstable fracture toughness to determine the initiation and instability of concrete fracture. Subsequently, researchers have studied all steps of crack development based on crack propagation criterions, which could reflect the association between the propagation force at the crack tip and the material’s own resistance during crack propagation.

To describe the strain softening behavior of concrete fracture, the concept of fracture process zone (FPZ) was established. By considering attractive atomic forces, Barenblatt [7,8] introduced FPZ, which was defined as a confined area near the crack tip. Furthermore, a different theory was proposed by Dugdale [9], which stated that the stress was equivalent to the yield strength of the material that acted across the crack within the plastic zone near the crack tip. This theory is applicable to ductile metal materials. Later, Hillerborg et al. [10] presented a fictitious crack model to represent the FPZ of quasi-brittle materials, where cohesions were distributed on both sides of the fictitious crack according to the settled softening constitutive relationship. The softening constitutive relationship of concrete, such as linear [10], bilinear [11], and exponential curves [12], described the characteristics of cohesive force on the FPZ. Fracture energy ${\mathrm{G}}_{\mathrm{f}}$ represented the envelope area under the softening constitutive curve and concluded that it was a necessary parameter in the numerical analysis of the concrete fracture process. It could be determined through experiments, such as direct tensile tests [13] or three-point bending beam tests [14].

The fictitious crack model has been widely adopted and applied in the crack propagation process simulation. Based on the fictitious crack model proposed by Hillerborg [10], Gerstle and Xie [15] simulated the crack propagation process according to the maximum tensile strength criterion by using a linear softening constitutive relationship. Carpinteri and Massabó [16] introduced a crack propagation criterion characterized by the stress intensity factor for the mode I fracture of cement-based materials, which can be expressed as Equation (1):

$${\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}-{\mathrm{K}}_{\mathrm{I}}^{\mathsf{\sigma }}=0,$$

(1)

where ${\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}$ represents the stress intensity factor generated by external force, and ${\mathrm{K}}_{\mathrm{I}}^{\mathsf{\sigma }}$ represents the stress intensity factor generated by cohesion. The criterion stated that when the crack was in a critical state, the difference between the stress intensity factor caused by external force and cohesion was zero at the crack tip. Ooi and Yang [17,18], Yang and Deeks [19], and Moës and Belytschko [20] applied this criterion to simulate the mode I and mix-mode fracture of reinforced concrete and plain concrete. However, concrete is a quasi-brittle material, and ${\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}{-\mathrm{K}}_{\mathrm{I}}^{\mathsf{\sigma }}$ is supposed to be a finite value. Based on this, Wu et al. [21] proposed an intensity-factor-based fracture propagation criterion where the difference between the stress intensity factor at the crack tip generated from external load and cohesion was greater than the initiation fracture toughness (${\mathrm{K}}_{\mathrm{Ic}}^{\mathrm{ini}}$), and crack started propagating. The propagation criterion can be shown as Equation (2):

$${\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}{-\mathrm{K}}_{\mathrm{I}}^{\mathsf{\sigma }}\ge {\mathrm{K}}_{\mathrm{Ic}}^{\mathrm{ini}},$$

(2)

This modified criterion is recently widely used in simulating the entire process of crack propagation, such as: three-point bending beams [21], infinite slab [22], concrete gravity dam model [23], concrete mode I [24], modes I–II [25] fracture of different strengths, and bimaterial interface crack propagation [26]. Meanwhile, numerical methods for calculating the FPZ length and ${\mathrm{K}}_{\mathrm{R}}$ curve are developed [21,22].

The main difference of the above literature is the utilization of different fracture propagation criterions in finite element numerical simulation to obtain the entire process of fracture propagation. What needs to be emphasized is that whichever fracture propagation criterion is used for numerical simulation, criterions are always based on the fictitious crack model, and the cohesion applied in the fictitious crack is an important link to affect the numerical simulations.

There are mainly two methods to apply cohesion in fictitious crack during crack propagation. The first method is to add the interfacial elements to characterize the cohesion in a separate crack, as shown in Figure 1a. Ingraffea et al. [27,28,29] used six-node interface elements to characterize fictitious crack for the analysis of complex fracture, but the nonlinear calculation efficiency was low. Swenson [30] used the six-node interface elements to simulate the dynamic crack propagation process. In the analysis of Bocca et al. [31] and Gerstle and Xie [15], a four-node linear displacement interface element was used to characterize the cohesion in the fictitious crack. Moreover, none of the above-mentioned finite element models can calculate or apply stress intensity factors.

The second method to apply cohesive stress is considering the calculation results of a finite element and directly applying cohesion to the corresponding crack surface nodes as a boundary condition [21,22,23,24,25,26], as shown in Figure 1b. It is noticeable that the cohesion applied in this process is not accurate enough. Crack opening displacement (COD) is obtained from external load, P, rather than the superposition of external load and the cohesion, P + σ, denoted by ${\mathrm{COD}}_{\mathrm{P}+\mathsf{\sigma }}{ < \mathrm{COD}}_{\mathrm{P}}$, where ${\mathrm{COD}}_{\mathrm{P}+\mathsf{\sigma }}$ is COD induced by P + σ and ${\mathrm{COD}}_{\mathrm{P}}$ is COD induced by P, so that the corresponding cohesion relationship is ${\mathsf{\sigma }}_{\mathrm{P}+\mathsf{\sigma }}{ > \mathsf{\sigma }}_{\mathrm{P}}$, where ${\mathsf{\sigma }}_{\mathrm{P}+\mathsf{\sigma }}$ is a recalculated cohesion under P + σ and ${\mathsf{\sigma }}_{\mathrm{P}}$ is a cohesion calculated from P based on the softening law, which indicates that the cohesion applied on a crack surface according to this method is less than the theoretical value. To obtain a more realistic cohesive stress value, it is plausible to use the multiple iterations of ${\mathrm{COD}}_{\mathrm{P}+\mathsf{\sigma }}$ and σ. However, it is incorrect to solve the nonlinear problem, and this iterative method often gives a nonconvergent result. Although the latter method considered the initial fracture toughness of concrete and used it as a control parameter for a crack propagation criterion, cohesion on the FPZ has not been applied accurately in the past research [21,22,23,24,25,26]. Therefore, other methods need to be explored to effectively solve the nonlinear problem.

In this paper, a numerical method is proposed, in which cohesive stress on FPZ in concrete is simulated and applied by a nonlinear spring element. The crack propagation criterion based on initial fracture toughness, which considers the stress singularity at the crack tip, is used. The proposed method improves the accuracy of cohesion applied on the FPZ and optimizes the numerical simulation process through displacement control. Subsequently, the P-delta curves, the P-CMOD curves, the variation of the FPZ length, and the ${\mathrm{K}}_{\mathrm{R}}$ curves are obtained. The simulated P-delta curves are performed to recalculate fracture energy and verify the accuracy of cohesion in the proposed method. The simulated P-CMOD curves are consistent with the previous experimental results, while the FPZ length and ${\mathrm{K}}_{\mathrm{R}}$ curves are different from the previous numerical results. The reasons that lead to the differences in results are discussed in detail. The comparative results in this paper will lead to a better insight, in which the accurate application of cohesive force on the FPZ has a crucial impact on the results of numerical simulations.

2. Methods and Materials

2.1. Nonlinear Spring Element

The user-defined nonlinear spring element of Ansys finite element software is used to complete the nonlinear finite element iterative solution process to achieve accurate cohesion applied on the concrete FPZ.

The unidirectional spring element with nonlinear generalized force deflection capability, combin39, can be utilized in the structural analysis [32]. This element can achieve longitudinal or torsional force in 1-D, 2-D, and 3-D application. In this paper, the element 1-D mode is used to apply the cohesive stress on the surface of a newly generated separate crack. The generalized force–deflection curve of this element is shown in Figure 2. The first and third quadrants indicate the tension and compression of the element, respectively. The element can output stress according to the user-defined force–deflection relationship curve.

The Combin39 element controls the mechanical behavior through key options. KEYOPT(1) controls the unloading path of the element. Since there is no COD reduction in the FPZ of concrete mode I cracks, 0 is assigned to KEYOPT(1) to unload according to the F-D loading curve. KEYOPT(2) controls the deformation behavior of the element under compression. The Combin39 elements used in this simulation provide tensile force on the surface of the newly generated separate crack to simulate cohesion that hinders the crack from opening and expanding. The cohesion does not transmit compressive stress, for KEYOPT(2), and assigns a value of 1 so that the element does not provide compressive stress under compression, as shown in Figure 3. KEYOPT(3) = 0 makes the element provide stress along the x-axis direction of the element itself. KEYOPT(4) = 0 makes the element provide 1-D F-D relationships.

A Combin39 element defines multiple points to represent the F-D relationship, so using this element can achieve stress softening behavior. The bilinear softening constitutive proposed by Petersson [11] is used in this paper to represent the relationship between cohesion and COD in the FPZ. This constitutive relationship is consistent with that used in previous research studies [21,22,23,24,25,26], so the results of this paper can be compared with those obtained before for discussion. The region under the softening constitutive curve represents the material’s fracture energy. The softening constitutive relationship is shown in Figure 4.

The corresponding formula is as follows:

$$\mathsf{\sigma }\left(\mathrm{w}\right)=\{\begin{array}{c}{ \mathrm{f}}_{\mathrm{t}}{-(\mathrm{f}}_{\mathrm{t}}{-\mathsf{\sigma }}_{\mathrm{s}}{\left)\right(\mathrm{w}/\mathrm{w}}_{\mathrm{s}}), 0 \le \mathrm{w} \le { \mathrm{w}}_{\mathrm{s}},\\ {\mathsf{\sigma }}_{\mathrm{s}}{(\mathrm{w}}_0{-\mathrm{w})/(\mathrm{w}}_0{-\mathrm{w}}_{\mathrm{s}}{), \mathrm{w}}_0 \le \mathrm{w} \le { \mathrm{w}}_{\mathrm{s}}\\ \hspace{1em}0\hspace{1em},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{w} \ge { \mathrm{w}}_0,\end{array},$$

(3)

where ${\mathrm{f}}_{\mathrm{t}}$ is the tensile strength of concrete, w is the opening displacement at any position between the newly generated crack surface, ${\mathrm{w}}_0$ is critical crack opening displacement, ${\mathsf{\sigma }}_{\mathrm{s}}$ and ${\mathrm{w}}_{\mathrm{s}}$ separately represent the cohesive force and crack opening displacement of the turning point in a bilinear softening constitutive relationship. The concrete tensile strength ${\mathrm{f}}_{\mathrm{t}}$ and fracture energy ${\mathrm{G}}_{\mathrm{f}}$ determine ${\mathsf{\sigma }}_{\mathrm{s}}$, ${\mathrm{w}}_0$, and ${\mathrm{w}}_{\mathrm{s}}$. The expressions are as follows:

$${\mathsf{\sigma }}_{\mathrm{s}}{=\mathrm{f}}_{\mathrm{t}}/3,$$

(4)

$${\mathrm{w}}_0=3{.6\mathrm{G}}_{\mathrm{f}}{/\mathrm{f}}_{\mathrm{t}},$$

(5)

$${\mathrm{w}}_{\mathrm{s}}=0{.8\mathrm{G}}_{\mathrm{f}}{/\mathrm{f}}_{\mathrm{t}}.$$

(6)

A point should be noted that the Combin39 element in Ansys requires the stiffness characterized by the first section in the first quadrant to be positive. To this end, a positive value must be assigned to the stiffness of the spring element. In addition, it is necessary to give the spring element an infinite initial stiffness, which is consistent with the softening constitutive relationship as in Figure 5. Therefore, the value corresponding to 0.1% of the critical crack opening displacement ${\mathrm{w}}_0$ on the softening constitutive curve is used to calculate the initial stiffness. In fact, during the nonlinear solution process, the relative displacements of the nodes at both ends of the spring element with an initial length of zero will be greater than 0.1% ${\mathrm{w}}_0$, and the effect of the first elastic section in the spring element constitutive relationship on the overall numerical simulation can be ignored.

2.2. Methods of Simulation

In this finite element analysis, the element in use is six-node triangular elements, the crack propagation length is set to 2 mm per step, and the gird near the crack needs to be encrypted. Different from the existing cohesive stress application method, after re-establishing and remeshing the model, the Combin39 nonlinear spring element should be added to connect the corresponding node on the crack surface, as shown in Figure 6. Due to the fine meshing of the model, a trapezoidal formula is used to apply nodal force at nodes of the triangular element instead of the nonuniform surface load of the element.

Then enter the solver, apply the load, and use the Newton–Raphson method to perform a nonlinear solution using line search. For nonlinear problems, Ansys uses Gaussian point results for calculation by default. At this time, the cohesive stress calculated is more accurate than the existing cohesion application method [21,22,23,24,25,26] under the same load constraint.

Equation (2) is applied as the crack propagate criterion [21] in this paper. The entire crack propagation process numerical simulation algorithm is described as follows:

• According to the test of geometric parameters and material parameters, a finite element model is established and meshed. The crack tip grid element should be performed to meet a singularity of −1/2.
• Apply displacement constraints to the finite element model, and calculate the stress intensity at the crack tip. When it reaches the initial fracture toughness, the crack starts to propagate.
• Crack spreads forward for a certain unit length Δa, and the finite element model is re-established and remeshed. Add the spring element between the newly generated interface.
• Use the displacement results from the previous step as the initial displacement conditions for this step.
• Determine whether the crack propagation criterion is satisfied. If not, increase the loading point displacement until the criterion is satisfied.
• Repeat steps 4 to 5 until the crack tip extends to the edge of the model, and the simulation of the entire crack propagation process ends. The calculation flow diagram is shown in Figure 7.

In the simulation, different from the load control method used in previous literatures [21,22,23,24,25,26], a displacement control method consistent with the real test is used. In a practical test represented by a three-point bending beam test, displacement-controlled loading is used, and the external load at the loading point is measured by a load sensor. Compared with the load control, the displacement control can quickly get a convergent solution from the finite element as a boundary condition. In contrast, load-controlled loading will lose the calculation accuracy in the nonlinear calculation process, which requires more iterations and even does not converge. In this paper, the proposed spring element is used to apply cohesive stress on the concrete FPZ to simulate the mode I crack propagation process of three-point bending beams.

2.3. Materials for Simulation

One set simulates concrete with different strengths, which is named C-series; the experimental data are from Dong [24] and Wang [33]; and concrete properties are shown in

Table 1. Material properties.

Concrete	${\mathbf{f}}_{\mathbf{c}}$ $\left(\mathbf{MPa}\right)$	${\mathbf{f}}_{\mathbf{t}}$ $\left(\mathbf{MPa}\right)$	$\mathbf{E}$ $\left(\mathbf{GPa}\right)$	${\mathbf{K}}_{\mathbf{Ic}}^{\mathbf{ini}}$ $(\mathbf{MPa}\mathbf{·}{\mathbf{m}}^{1/2})$	${\mathbf{G}}_{\mathbf{f}}$ $(\mathbf{N}/\mathbf{m})$
C20	32.8	3.05	29.9	0.461	117.1
C40	48.9	3.74	33.2	0.616	124.5
C60	69.9	4.43	35.7	0.632	114.9
C80	84.1	5.01	38.1	0.667	120.5
C100	115.8	5.71	41.4	0.917	115.4

(where ${\mathrm{f}}_{\mathrm{c}}$ represents the compressive strength, and E represents the elastic modulus). The size of the specimen in this set is S × D × B = 480 mm × 120 mm × 60 mm, and the initial crack ratio ${\mathrm{a}}_0/\mathrm{D}$ is set to 0.3, as shown in Figure 8.

Another set of data is from Dong [21]. This set is used to simulate beams with a different initial crack ratio ${\mathrm{a}}_0/\mathrm{D}$ and depth D. The specimens in the literature [21] were divided into two series, named B-series and L-series, and the same series names were still used in this paper. For the B-series, the size of specimens is retained at S × D × B = 600 mm × 150 mm × 40 mm, but the initial crack ratio ${\mathrm{a}}_0/\mathrm{D}$ varies from 0.2 to 0.6 (

Table 2. Size and material parameters for B-series beams.

Specimen	$$\begin{gathered} \mathbf{S} \times \mathbf{D} \times \mathbf{B} \\ \left(\mathbf{mm}\right) \end{gathered}$$	${\mathbf{a}}_0/\mathbf{D}$	${\mathbf{K}}_{\mathbf{IC}}^{\mathbf{ini}}$ $(\mathbf{MPa}\mathbf{·}{\mathbf{m}}^{1/2})$	${\mathbf{G}}_{\mathbf{f}}$ $(\mathbf{N}/\mathbf{m})$
B2-1	$600\times 150\times 40$	0.2	0.60	96
B2-2			0.59	100
B2-3			0.62	92
Avg.			0.60	96
B3-1	$600\times 150\times 40$	0.3	0.65	100
B3-2			0.63	105
B3-3			0.53	88
Avg.			0.60	98
B4-2	$600\times 150\times 40$	0.4	0.58	85
B4-3			0.61	100
Avg.			0.60	92.5
B6-1	$600\times 150\times 40$	0.6	0.60	105
B6-2			0.62	120
B6-3			0.58	100
Avg.			0.60	108

), and the material mechanical properties are ${\mathrm{f}}_{\mathrm{t} }$ = 2.4 MPa and E = 28 GPa. For L-series specimens, the span-to-depth ratio S/D is set to 4, the initial crack ratio ${\mathrm{a}}_0/\mathrm{D}$ is set to 0.4, the depth D varies from 100 mm to 300 mm (

Table 3. Size and material parameters for L-series beams.

Specimen	$$\begin{gathered} \mathbf{S} \times \mathbf{D} \times \mathbf{B} \\ \left(\mathbf{mm}\right) \end{gathered}$$	${\mathbf{a}}_0/\mathbf{D}$	${\mathbf{K}}_{\mathbf{Ic}}^{\mathbf{ini}}$ $(\mathbf{MPa}\mathbf{·}{\mathbf{m}}^{1/2})$	${\mathbf{G}}_{\mathbf{f}}$ $(\mathbf{N}/\mathbf{m})$
L1-1	$400\times 100\times 100$	0.4	0.52	104
L1-2			0.52	90
Avg.			0.52	97
L2-1	$800\times 200\times 100$	0.4	0.71	155
L2-2			0.62	151
Avg.			0.67	153
L3-1	$1200\times 300\times 100$	0.4	0.78	123
L3-2			0.75	146
L3-3			0.75	153
Avg.			0.76	141

), and the material mechanical properties are ${\mathrm{f}}_{\mathrm{t}}$= 2.3 MPa and E = 24 GPa.

3. Results and Discussion

The simulation results of the crack propagation process based on the proposed method are compared with the previous experimental results and previously simulated results. In the same series, the crack propagation criterion, the softening constitutive, the material mechanical properties, and the specimen size are all the same except for the cohesion application method. In this section, the numerical method based on applying cohesion by the Combin39 nonlinear spring element is called iteration, and the previous numerical method that directly applies cohesion is called no iteration. The results will be discussed through four aspects: P-delta curve, P-CMOD curve, FPZ length, and ${\mathrm{K}}_{\mathrm{R}}$ curve.

3.1. P-Delta Curve

The P-delta curves were not provided in the previous study, so the comparison between simulation results and experimental results cannot be performed. However, simulation results obtained from different methods can be compared. Fracture energy is represented by cohesion on the FPZ, and the fracture energy can be calculated by Equation (7):

$${\mathrm{G}}_{\mathrm{f}}=\frac{{\mathrm{W}}_0+\mathrm{m}\mathrm{g}{\mathrm{\delta }}_0}{{\mathrm{A}}_{\mathrm{l}\mathrm{i}\mathrm{g}}},$$

(7)

where ${\mathrm{W}}_0$ represents the external load work, which is equal to the enveloped area under the P-delta curve; $\mathrm{m}\mathrm{g}{\mathrm{\delta }}_0$ represents the work done by specimen self-weight; and ${\mathrm{A}}_{\mathrm{l}\mathrm{i}\mathrm{g}}=\left(\mathrm{D}-{\mathrm{a}}_0\right)\times \mathrm{B}$ represents the ligament area. Therefore, fracture energy can be recalculated from the P-delta curves obtained by numerical simulation and compared with the input values to verify the accuracy of the applied cohesion with iteration. The P-delta curves obtained by the iterative and noniterative methods are plotted in Figure 9, and the plots enable us to find that the enveloped areas under the P-delta curves are significantly different under the two different methods. The fracture energy is recalculated according to Equation (7), and comparison results are shown in

Table 4. Recalculated fracture energy comparison table.

Group	$\mathbf{m}\mathbf{g}$ $\left(\mathbf{N}\right)$	${\mathbf{\delta }}_0$ $\left(\mathbf{mm}\right)$	${\mathbf{W}}_0$ $(\mathbf{N}\mathbf{·}\mathbf{mm})$	${\mathbf{A}}_{\mathbf{l}\mathbf{i}\mathbf{g}}$ $(\mathbf{mm} \times \mathbf{mm})$	${\mathbf{G}}_{\mathbf{f}}-\mathbf{Recalc}$ $(\mathbf{N}/\mathbf{m})$	${\mathbf{G}}_{\mathbf{f}}$ $(\mathbf{N}/\mathbf{m})$	%
C20-Iteration	81.37	0.60	538.5	$84 \times 60$	116.5	117.1	−0.5%
C20-No Iteration		0.22	224.2		48.0		−59.0%
C40-Iteration	81.37	0.55	592.5	$84 \times 60$	126.4	124.5	1.5%
C40-No Iteration		0.24	286.0		60.6		−51.3%
C60-Iteration	81.37	0.48	553.0	$84 \times 60$	117.5	114.9	2.2%
C60-No Iteration		0.23	283.5		60.0		−47.8%
C80-Iteration	81.37	0.45	581.3	$84 \times 60$	122.6	120.5	1.7%
C80-No Iteration		0.23	303.0		63.8		−47.0%
C100-Iteration	81.37	0.30	571.6	$84 \times 60$	118.3	115.4	2.5%
C100-No Iteration		0.26	374.4		78.5		−32.0%
B2-Iteration	84.76	0.46	438.2	$120 \times 40$	99.5	96.0	3.6%
B2-No Iteration		0.26	235.8		53.7		−44.0%
B3-Iteration	84.76	0.52	385.9	$105 \times 40$	102.4	98.0	4.5%
B3-No Iteration		0.27	203.2		53.8		−45.1%
B4-Iteration	84.76	0.46	300.2	$90 \times 40$	94.1	92.5	1.8%
B4-No Iteration		0.29	164.3		52.5		−43.3%
B6-Iteration	84.76	0.48	207.1	$60 \times 40$	106.8	108.0	−1.1%
B6-No Iteration		0.30	104.9		54.3		−49.7%
L1-Iteration	94.18	0.45	523.0	$60 \times 100$	94.2	97.0	−2.9%
L1-No Iteration		0.22	256.7		46.2		−52.3%
L2-Iteration	376.7	0.66	1555.2	$120 \times 100$	150.3	153.0	−1.8%
L2-No Iteration		0.44	831.3		82.1		−45.7%
L3-Iteration	847.6	0.75	2178.4	$180 \times 100$	156.2	141.0	10.8%
L3-No Iteration		0.50	1308.6		96.2		−31.7%

* ${\mathrm{G}}_{\mathrm{f}}$-recalc means recalculated fracture energy; % means the difference percentage of ${\mathrm{G}}_{\mathrm{f}}$-recalc from the material’s ${\mathrm{G}}_{\mathrm{f}}$.

.

In Table 4, it can be concluded that the fracture energy calculated by the iterative method is close to the input values, while the fracture energy calculated by the noniterative method is much smaller than the input values. This conclusion is also consistent with the theoretical analysis, indicating that directly applying the cohesion approach results in a smaller value than the true cohesion value. In other words, the accuracy of cohesion can be improved by the proposed method, where cohesion is applied by the spring element.

3.2. P-CMOD Curve

The test and simulated P-CMOD results of the C-series are plotted in Figure 10.

It is worth noting from Figure 10 that, except for the C20 group, the simulated curves by iteration fit the experimental curves perfectly. In each set of curves, the peak load of the simulation obtained by no iteration is less than that obtained by iteration. This is because without iteration, the COD result from external load is larger than the COD obtained from the superposition of external load and cohesion so that the corresponding cohesive force is less based on the softening constitutive. In contrast, the iteration method can get a larger but more accurate cohesion, so it can get a larger peak load according to the crack propagation criterion in use. The descending section of the curve obtained by the iteration method agrees better with test curves compared with the no iteration method. Previous simulation results [24] are plotted in Figure 10 in green. From the comparison of the results of the P-CMOD curves, it can be found that previous simulated curves given by Dong are close to the result of no iteration.

The experimental and simulated curves of the B-series and L-series are plotted in Figure 11 and Figure 12. P-CMOD curves of no iteration are consistent with previous simulated results, and only the P-CMOD curves from [21] are drawn in these figures for comparison.

The first linear segment in Figure 11 and Figure 12 represents the elastic deformation behavior of concrete. Using the point on the linear segment on the test curve, the elastic modulus of concrete can be calculated by Equation (8) [34]:

$$\mathrm{CMOD}=\frac{24\mathrm{P}\mathsf{\lambda }}{\mathrm{EB}}\left[0.76-2.28\mathsf{\lambda }+3.87{\mathsf{\lambda }}^2-2.04{\mathsf{\lambda }}^3+\frac{0.66}{{\left(1-\mathsf{\lambda }\right)}^2}\right],$$

(8)

where λ is equal to ${(\mathrm{a}+\mathrm{H}}_0{)/(\mathrm{d}+\mathrm{H}}_0)$, ${\mathrm{H}}_0$ is the thickness of the knife edge holding the clip gauges used to measure CMOD and equal to 2 mm [24], P and CMOD are data value chosen on the linear segment in the test curve, and B is the thickness of the test specimen. The difference of the linear segment between the numerical simulation and the test indicates that the elastic modulus of the concrete measured by the test is not accurate enough. However, material parameters from Dong [21] are still used in the simulation for comparison. It can be seen from figures that the simulation has the same trend as the previous set.

Conclusions can be drawn that improving the accuracy of cohesion in a fictitious crack has a positive effect on the numerical simulation, and a curve obtained by the iterative method gives a larger peak load. The criterion based on material initial fracture toughness with iterative cohesion is suitable for simulating the P-CMOD curves of concrete with a different strength, ${\mathrm{a}}_0/\mathrm{D}$, and size.

3.3. FPZ Length

The simulated FPZ length of the C-series as the crack propagation process is shown in Figure 13, where the horizontal axis is the ratio of the fictitious crack propagation length to the ligament length, expressed by ${\mathsf{\Delta }\mathrm{a}/(\mathrm{D} - \mathrm{a}}_0)$. Only the C20, C60, and C100 FPZ length results are given by Dong [24]. For the B-series and L-series, the results are plotted in Figure 14, the horizontal axis represents the ratio of the crack length to the depth of the beam expressed by a/D, and the vertical axis represents the ratio of the sum of the initial crack length and FPZ length to the beam depth expressed by ${(\mathrm{a}}_0{+\mathrm{a}}_{\mathsf{\sigma }})/\mathrm{D}$.

The reason for choosing different coordinate systems in the comparison is to facilitate comparison with the original literature. From the comparison results in Figure 13, it can be found that after reaching the maximum load corresponding to the P-CMOD curve, the FPZ length simulated by a different method keeps growing with the crack propagation. The iteration method obtains a greater maximum length value compared with the no iteration method, as shown in Figure 13 and Figure 14. The cohesion will appear at the region with COD less than ${\mathrm{w}}_0$ after multiple iterations, while cohesion will not exist at a part of the same region whose COD is greater than the critical value without iteration. After FPZ hits the maximum value, the FPZ length with iteration decreases more significantly. The numerical simulated results simulated by Dong et al. [21,24] were compared with the numerical simulated results in this paper. The variation of the C20 group is close to the iteration result, but the FPZ length of the C60 and C100 groups exceeds the iteration results. It may suggest that a certain iterative process during simulation process, “repeatedly solving COD and σ”, is conducted when calculating the FPZ length. However, the results are unstable and not accurate enough for the C-series.

For the B-series and L-series, FPZ lengths obtained from Dong et al. [21] are smaller than the results of the iterative ones. In the process of the iteration method, after adding the Combin39 nonlinear spring element, only boundary conditions of the model are needed. The load P, COD, and σ are solutions of nonlinear equations so that the cohesion applied in the FPZ can be characterized more accurately by the iteration method. Next, factors that affect the FPZ length based on the iterative results are figured out as follows.

For the C-series, it is found that with the increase in concrete strength, the maximum FPZ length decreases. Based on the softening law, fracture energy and tensile strength are the main factors in determining critical crack opening displacement. Small critical crack opening displacement results in small FPZ length. The fracture energy for the C-series is close, so the concrete with higher tensile strength results in a shorter length of the FPZ. For the B-series, as ${\mathrm{a}}_0/\mathrm{D}$ decreases, the maximum FPZ length increases since the beam with smaller ${\mathrm{a}}_0/\mathrm{D}$ has a longer ligament length, which can make the FPZ fully develop. For the L-series, as the depth of the beam increases, the maximum FPZ length increases, and the principle is the same as the B-series. From another perspective, it can be concluded that the geometry shape of beams affects the FPZ length.

3.4. ${\mathrm{K}}_{\mathrm{R}}$ Curve

The ${\mathrm{K}}_{\mathrm{R}}$ curves are used to represent the change in fracture toughness in the process of crack propagation. Foote et al. [35] proposed a theoretical analysis model of the ${\mathrm{K}}_{\mathrm{R}}$ curve before the FPZ fully developed. The theoretical model for analyzing the ${\mathrm{K}}_{\mathrm{R}}$ curve proposed by Hu and Wittmann [36] accurately restored the wedge opening loaded mortar specimens’ test results. Xu and Reinhardt [37] proposed a ${\mathrm{K}}_{\mathrm{R}}$ calculation model based on cohesion in the FPZ where ${\mathrm{K}}_{\mathrm{R}}$ was the stress intensity factor generated by initial fracture toughness ${\mathrm{K}}_{\mathrm{IC}}^{\mathrm{ini}}$ and cohesion in the FPZ together. However, an assumption was made that the FPZ length and the cohesive stress distribution stayed unchanged after the FPZ fully developed. Therefore, the ${\mathrm{K}}_{\mathrm{R}}$ curve increased as the ratio of the effective crack length to the beam depth increased. Lutz and Swain [38,39] found that the ${\mathrm{K}}_{\mathrm{R}}$ resistance curve for ceramic brittle materials increased as the increase in the FPZ length. After the FPZ hit the maximum value, ${\mathrm{K}}_{\mathrm{R}}$ remained a stable value. This result was consistent with the conclusion explored by Xu et al. [40]. Dong et al. [21] proposed another method for numerically calculating the ${\mathrm{K}}_{\mathrm{R}}$ curve, in which ${\mathrm{K}}_{\mathrm{R}}$ is calculated by:

$${\mathrm{K}}_{\mathrm{R}}\left(\mathsf{\Delta }\mathrm{a}\right){=\mathrm{K}}_{\mathrm{Ic}}^{\mathrm{ini}}{+\mathrm{K}}_{\mathrm{I}}^{\mathsf{\sigma }}{=\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}\left(\mathrm{P},\mathsf{\Delta }\mathrm{a}\right),$$

(9)

where $\mathsf{\Delta }\mathrm{a}$ means the crack propagation length in numerical simulation, and ${\mathrm{K}}_{\mathrm{R}}\left(\mathsf{\Delta }\mathrm{a}\right)$ and ${\mathrm{K}}_{\mathrm{I}}^{\mathrm{P}}\left(\mathrm{P},\mathsf{\Delta }\mathrm{a}\right)$ represent the stress intensity factor of the crack extension resistance and external load, respectively. In other words, ${\mathrm{K}}_{\mathrm{R}}$ is equivalent to the stress intensity factor corresponding to the external load during crack propagation. This paper uses the proposed spring element model to calculate ${\mathrm{K}}_{\mathrm{R}}$ curves based on Equation (9). The simulated ${\mathrm{K}}_{\mathrm{R}}$ curves’ results of the B-series and L-series from this paper and previous simulation results [21] are plotted in Figure 15 for comparison.

In the conclusion of a previous research [21], there was a plateau in the ${\mathrm{K}}_{\mathrm{R}}$ curve with the FPZ length variation. However, it can be found that there is no plateau existing in the ${\mathrm{K}}_{\mathrm{R}}$ curve by iterative approach in Figure 15. In comparison with Figure 13 and Figure 14, the ${\mathrm{K}}_{\mathrm{R}}$ curve still exhibits an increasing trend after the FPZ length has fully developed, and only the raising rate becomes slightly decreased. To further explore the relationship between the FPZ length and the ${\mathrm{K}}_{\mathrm{R}}$ curve, this paper draws the FPZ length and ${\mathrm{K}}_{\mathrm{R}}$ curve of the L3 specimen together in Figure 16. From Figure 16, it can still be inferred that the FPZ length affects the shape of the ${\mathrm{K}}_{\mathrm{R}}$ curve. After the FPZ length has fully developed, the ${\mathrm{K}}_{\mathrm{R}}$ curve growth rate decreases. The value of the ${\mathrm{K}}_{\mathrm{R}}$ curve simulated in this paper is larger than that from previous simulations [21], and at the end of the curves, the ratio of the ${\mathrm{K}}_{\mathrm{R}}$ values obtained by two methods is over two times. The difference in results verifies that it is essential to apply cohesive force in FPZ considering iteration.

It is necessary to reiterate that for the B-series and L-series, each set shares the same material mechanical parameter except for ${\mathrm{G}}_{\mathrm{f}}$ and ${\mathrm{K}}_{\mathrm{Ic}}^{\mathrm{ini}}$. It is more appropriate to explore the factors that affect the ${\mathrm{K}}_{\mathrm{R}}$ curve from these two series first. According to Figure 14b and Figure 15b and Table 3, it can be found that the fracture energy ${\mathrm{G}}_{\mathrm{f}}$ affects the shape of the ${\mathrm{K}}_{\mathrm{R}}$ curve. Although the depths of beams are different from the L-series, the FPZ length variation of the L1 and L2 specimens is similar from overall view in Figure 15b. The reason why the ${\mathrm{K}}_{\mathrm{R}}$ curve value of L2 is greater than L1 is that ${\mathrm{G}}_{\mathrm{f}}$ of L2 is larger than that of L1. Comparing the ${\mathrm{K}}_{\mathrm{R}}$ curve of L1 with L3, though ${\mathrm{G}}_{\mathrm{f}}$ of L3 is greater than that of L1, the FPZ length variation is the main reason for L3 getting a smaller ${\mathrm{K}}_{\mathrm{R}}$ value at the end of the curve at this moment. For the B-series as in Figure 14a, the initial crack ratio ${\mathrm{a}}_0/\mathrm{D}$ of specimens is different so the ${\mathrm{K}}_{\mathrm{R}}$ curve starts at different points, but at the end of the curve, ${\mathrm{G}}_{\mathrm{f}}$ is the main reason for the difference in the ${\mathrm{K}}_{\mathrm{R}}$ value.

The ${\mathrm{K}}_{\mathrm{R}}$ curves of the C-series are plotted in Figure 17. The ${\mathrm{K}}_{\mathrm{R}}$ value of higher-strength concrete is the highest at first, but then it becomes the smallest. It seems like the tensile strength affects the ${\mathrm{K}}_{\mathrm{R}}$ value, and the concrete with a higher tensile strength will get a lower ${\mathrm{K}}_{\mathrm{R}}$ value at the end of the curve. According to the data in Table 1, it is found that the fracture energy of concrete with different strengths obtained from the experiment is relatively close. However, as discussed above about the influencing factors for the FPZ length of the C-series, higher tensile strength results in smaller maximum FPZ length under the same fracture energy. Therefore, it can be concluded that the FPZ length is the main reason for the variation in the ${\mathrm{K}}_{\mathrm{R}}$ value of the C-series.

4. Conclusions

Aiming at the numerical simulation of the concrete mode I crack propagation process, a numerical method using the Combin39 nonlinear spring element to apply cohesive stress on the FPZ is proposed in this paper. Based on displacement control, the values of cohesion on the FPZ are nonlinear solutions according to the Newton–Raphson method using line search. According to the comparison between the above numerical calculation results and test results, conclusions are drawn as follows:

• The P-delta curves obtained from numerical simulation were used to recalculate fracture energy and compared with the input value. The direct application of cohesion produced results that were much smaller than the true value, while the proposed method of applying cohesion by the spring element can improve the accuracy of the applied cohesion.
• Using the proposed method to apply the cohesive stress on the FPZ for numerical calculation, though the peak load obtained is larger than the no-iteration result, P-CMOD curves obtained are still fitted the experimental results, indicating that the cohesion is applied in a reasonable manner, and the crack propagation criterion based on material initial fracture toughness is suitable for numerical simulation.
• The FPZ length obtained by the iteration method can reach a larger maximum value compared with no-iteration, and the decline gradient of the FPZ length also becomes larger after the FPZ is fully developed. The fracture energy, tensile strength, and geometry shape of the beam are main reasons for deciding the FPZ length. According to the bilinear softening constitutive, fracture energy and tensile strength determine the critical crack opening displacement. A larger critical crack opening displacement or a longer ligament makes it easier to obtain a larger FPZ length.
• The ${\mathrm{K}}_{\mathrm{R}}$ curve obtained with iteration is significantly different from the noniterative curve. The ${\mathrm{K}}_{\mathrm{R}}$ value is continually rising with the crack propagation. After the FPZ length has fully developed, the rising rate of the ${\mathrm{K}}_{\mathrm{R}}$ curve has become slow. Through the synergistic comparison with the change of the FPZ length, it is proved that the ${\mathrm{K}}_{\mathrm{R}}$ value is mainly affected by the fracture energy and FPZ length. The iterative method is more suitable for simulating the ${\mathrm{K}}_{\mathrm{R}}$ curves by improving the accuracy of the cohesion on the FPZ.

The significant difference in numerical results indicates that applying cohesion is essential for numerical simulation, and the accuracy of cohesion has a great influence on the results of numerical simulation. The method has been proved to be well applicable to numerical simulation for the concrete mode I crack propagation process by improving the accuracy of cohesion applied on concrete FPZ. According to the crack propagation criterion based on initial fracture toughness, more numerical simulations are needed to explore the effect of the cohesion application in the modes I–II crack propagation process, which is one of our future studies.