The crack propagation process at the crack tip of solid material is schematically illustrated in Figure 1A. Accordingly, the mechanical behavior of crack propagation can be described using a cohesive interface model, which is usually abstracted as a cohesive element (see Figure 1B). The stress-displacement relationship of a cohesive interface can be described using the bilinear traction-separation constitutive model, as shown in Figure 1C. This constitutive model is divided into linear stress growth and damage softening stages divided by peak strength. In the two-dimensional problem, Eq. 1 can be used to describe the traction-separation behavior in the pre-peak stage [22]: { σ n τ s } = [ E nn 0 0 E ss ] { ε n ε s } = 1 T 0 [ E nn 0 0 E ss ] { δ n δ s } , (1) where σ n and τ s are the nominal traction components: normal and shear stresses, respectively. ε n and ε s are the two components of the nominal strains, respectively. δ n and δ s are the separation components: normal and shear displacements, respectively. T 0 denotes the original thickness of the cohesive element. E nn and E ss are the initial normal and shear stiffness of the interface, respectively.

When the interface displacement increases to δ n 0 and δ s 0 , the interface model is damaged. According to Eq. 1, the displacements of damage initiation can be calculated as follows: δ n 0 = σ nc K nn   a n d   δ s 0 = τ sc K ss , (2) where σ nc and τ sc are normal and shear stresses when the cohesive interface reaches the Mode I and Mode II fracture, respectively.

According to the bilinear traction-separation constitutive law shown in Figure 1C, the complete separation displacements δ n f and δ s f of the cohesive interface can be calculated using the following formula: δ n f = 2 G n σ nc   a n d   δ s f = 2 G s τ sc (3) where G n and G s are the critical fracture energy of Mode I and Mode II fracture at the cohesive interface, respectively.

When the traction nominal stress on the cohesive interface meets the maximum stress criterion, it can be considered that the cohesive interface is damaged. The maximum nominal stress criterion of the cohesive interface is used to represent the damage initiation condition, which can be written as: max { 〈 σ n 〉 σ nc , 〈 τ s 〉 τ sc } = 1 , (4) where 〈 〉 is the Macaulay bracket. When x ≥ 0 , 〈 x 〉 = x . When x < 0 , 〈 x 〉 = 0 .

When the damage condition described in Eq. 4 are met, the traction stress on the cohesive interface enters the post-peak linear softening stage. A scalar damage variable is used to describe the linear reduction of traction stress at the cohesive interface, which is expressed as follows: σ n = ( 1 − D ) σ ¯ n   a n d   τ s = ( 1 − D ) τ ¯ s (5) where σ ¯ n and τ ¯ s are nominal traction normal and shear stresses predicted using the traction-separation law without damage.

For the mixed fracture failure mode, as shown in Figure 1D, the bilinear traction-separation law of the cohesive interface is described by the effective relative displacement δ m , which is defined as: δ m = 〈 δ n 〉 2 + δ s 2 . (6) For linear damage softening, the damage variable can be expressed using the following formula [23]: D = δ m f ( δ m − δ m 0 ) δ m ( δ m f − δ m 0 ) , (7) where δ m 0 and δ m f are the effective displacements at damage initiation and complete debonding, respectively.

For the mixed-mode fracture, the mode combination of the cohesive interface can be quantified by defining the relative ratio of normal fracture and shear fracture energies. The fracture energy G n represents the corresponding traction nominal stress work to its corresponding separation displacement. Thus, the mixed fracture energy G t = G n + G s can be defined. Therefore, the ratios m 1 and m 2 of normal fracture energy G n to shear fracture energy G s can be defined as: m 1 = G n G t   a n d         m 2 = G s G t . (8)

The damage evolution of the cohesive interface can be defined according to the fracture energy, which is equal to the area of the geometric region under the traction-separation curve, which can be expressed by G n and G s , respectively. Based on the power-law fracture criterion, the dependence of fracture energy on mixed-mode can be defined. The power-law fracture criterion is expressed as: { G s G n c } α + { G s G s c } α = 1 , (9) where G n C , G s C and α are the material constants for describing the fracture behavior of the cohesive interface.

The heat conduction between the upper and lower surfaces of the cohesive interface can be considered as a function of the temperature difference between the upper and lower interfaces, which is defined by the following formula: q = k cz ( θ + − θ − ) , (10) where q is the heat flux per unit area crossing the cohesive interface. θ + and θ − are the temperatures on the top and bottom surfaces of the cohesive interface. k cz is the gap conductance of the cohesive zone. Its value depends on the normal separation displacement, i.e. k cz = k cz ( δ n ) . In this study, it is considered that the relationship between gap conductance k cz and normal separation displacement δ n can be described using the following piecewise function. k cz = { k cz0 − k δ n ,             δ n < k cz0 / k 0   ,                 δ n ≥ k cz0 / k , (11) where k cz0 is the thermal conductivity of cohesive interface without any normal separation displacement ( δ n = 0), which can be considered to have the same value as the thermal conductivity of surrounding solid material. k is the coefficient that the thermal conductivity of the cohesive region decreases with the increase of separation displacement δ n . When the separation displacement δ n exceeds the maximum crack width k cz0 / k , k cz = 0 . It is considered that there is no heat conduction between the open cohesive interfaces.

Concrete is a heterogeneous composite material. When studying its mechanical behavior, the heterogeneity of concrete can be simplified as a two-phase medium constituting aggregate and mortar matrix. Figure 2A shows the geometric plane model of the concrete specimen, a square area with a side length of 150 mm. The model’s right horizontal and vertical-up directions are the positive x-axis and y-axis, respectively.

In this plane area of the geometric model, as shown in Figure 2A, the computer reconstruction technology proposed in the articles [26, 27] is used to generate aggregates randomly. This study uses polygons to reconstruct the geometric model of aggregates. In order to control the calculation scale, the equivalent diameter of the simulated minimum aggregate particle is greater than 2.36 mm. Aggregates smaller than this size are considered part of the mortar matrix.

In this model, the aggregates are divided into five groups according to the size range in the generated geometric model. The total volume fraction of aggregates in the whole model area is 30%. The aggregate size ranges are (26.5,19.0) with volume fraction of 2.37% (19.0,16.0) with volume fraction of 6.78% (16.0,9.5) with volume fraction of 12.96% (9.5,4.75) with volume fraction of 6.66% and (4.75, 2.36) with volume fraction of 1.23%.

In this study, Abaqus is used to carry out coupled thermal-displacement analysis to simulate the thermal-induced spalling failure of the concrete specimen in the high-temperature furnace. The geometric plane model is meshed using the first-order linear plane strain quadrilateral element with a temperature degree of freedom, as shown in Figures 2B,C. The element type is CPE3T. This numerical model has 4244 CPE3T elements and 10,708 nodes.

This study uses the first-order cohesive element, including temperature degrees of freedom, to simulate the spalling failure behavior caused by the thermal stress of concrete under high temperatures. The aggregates are rock particles whose strengths are significantly higher than the mortar matrix. Thus, it is reasonable to consider that the high-temperature spalling failure of concrete occurs only in the cement mortar matrix. Therefore, the cohesive elements are embedded into the area of the cement mortar matrix, as shown in Figure 2D. This model embeds 5,193 cohesive elements, including 10,517 nodes. The cohesive element type is COH2D4T.

In this numerical model, three materials are involved. They are aggregate, cement mortar and cohesive interfaces. The thermal and mechanical behaviors of aggregate and cement mortar are described using a thermo-elastic model, which contains six material parameters: density ρ , elastic modulus E , Poisson’s ratio ν , thermal expansion coefficient α , thermal conductivity λ and specific heat capacity c p .

In this study, the high-temperature experiments of concrete materials were not carried out. By referring to the research results of the articles [16, 17, 19, 28], the thermal and mechanical parameters for the numerical simulation are obtained. For aggregate, they are ρ = 2650 kg ⋅ m - 3 , E = 45 GPa , ν = 0.2 , α = 2.0 × 10 − 5 , λ = 3.0 W ⋅ m - 1 ⋅ K - 1 and c p = 750 J ⋅ kg - 1 ⋅ K - 1 . For cement mortar, they are ρ = 2200 kg ⋅ m - 3 , E = 22 GPa , ν = 0.22 , α = 1.5 × 10 − 5 , λ = 1.5 W ⋅ m - 1 ⋅ K - 1 and c p = 1200 J ⋅ kg - 1 ⋅ K - 1 .

The mechanical behavior of the cohesive interface is described using the traction-separation constitutive model. This model has two elastic parameters: normal traction modulus E nn = 0.5 GPa and tangential traction modulus E ss = 0.2 GPa , and five damage parameters: nominal traction normal stress σ nc = 1.8 MPa , nominal traction shear stress τ sc = 5.2 MPa , critical fracture energy G n = 150 N / m for Mode I fracture, critical fracture energy G s = 1500 N / m for Mode II fracture, and the power exponent η = 0.97 for the power-law fracture criterion. The cohesive interface model has two thermal parameters: initial gap thermal conductivity k cz0 = 1.5 W ⋅ m - 1 ⋅ K - 1 and gap thermal conductivity reduction coefficient k = 0.3 .

In order to simulate the thermal expansion deformation, two nodes, A and B, are selected in the central area of the numerical model. Constrain the displacement in the x-axis direction of point A and the displacement in the y-axis direction of point B.

In this study, ISO834 standard fire curve is used as the temperature load, which is described using the following formula: θ = θ 0 + 345 ⁡ log 10 ( 8 t / 60 + 1 ) (12) where θ 0 is the initial ambient temperature. In this study, θ 0 = 20 C ∘ . t is the heating time in second. θ is the temperature at the current time.

An initial temperature field θ 0 = 20 C ∘ is applied on the nodes of all elements in the numerical model. The temperature load of the standard fire curve is applied on the nodes in the four boundaries of the model.

In this study, the simulation time is 1800 s. That is, the concrete specimen is heated for 30 min. The numerical calculation cannot converge when the heating time reaches 1,399 s. This shows that the concrete specimen exhibits a severe thermal-induced explosive spalling failure at this time.

Figure 3 shows the deformations of the concrete specimen under high temperatures and the damage states of cohesive elements. It can be seen that the deformation of the concrete material increases gradually. At the same time, the damage of cohesive elements in mortar increases, and the thermal-induced damage gradually extends from the outer boundaries to the interior of the model. The damage evolution of the cohesive elements shows that the thermal cracks originate at the aggregate boundaries. Thermal cracks are easier to sprout at aggregates-mortar interfaces due to different thermal expansion coefficients. The initiated thermal cracks propagate along with the interfaces.

When t = 360 s, as shown in Figure 3A-1, the displacement field of the concrete material is roughly centrosymmetric. Nevertheless, the aggregate distribution affects the centrosymmetry of the displacement field. At the four corners, the thermal-induced displacement values are in a range of 0.38–0.45 mm. Figure 3B-1 shows that, in the local areas about 2–5 cm away from the four corners, the values of damage variable D are about 0.57–0.86. It shows that the cohesive elements in these local areas have been damaged and deteriorated, and thermal cracks begin to propagate.

Figure 3A-2 shows that, when t = 720 s, the centrosymmetry of the displacement field is weakened by the distribution of aggregates. The displacement values at the four corners increased to 0.64–0.77 mm. The displacements in the upper-left and lower-left corners increase rapidly, related to the damage states of the cohesive elements, as shown in Figure 3B-2. It can be seen that the damaged area gradually extends from the four corners towards the depths of about 4–6.5 cm. The fractured cohesive elements were observed at three positions: 6.6 cm-depth from the upper-left corner, 5.8 cm-depth from the upper-right corner and 6.5 cm-depth from the lower-right corner. The damage variable values of cohesive elements at these three positions are 0.90, 0.91 and 0.98, respectively. It can be seen that the cohesive elements at these three positions have been completely destroyed to form thermal cracks. These crack lengths are about 13.6, 8.9 and 39.8 cm, respectively.

When t = 1080 s, as shown in Figure 3A-3, the displacements are in a range of 0.85–1.00 mm observed in the upper-left, lower-left and upper-right corners. The maximum displacement appears in the lower-right corner, and its value is about 2.0–2.4 mm. The thermal-induced displacement of the concrete specimen has completely lost its centrosymmetry. The reason is related to the propagation and evolution of thermal cracks, as shown in Figure 3B-3. It can be seen that when the time reaches 1080 s, some thermal cracks gradually form in the model. The longest thermal crack appears in the lower-right corner, and its middle is about 6.5 cm away from the lower-right corner. The two ends of this crack are 4.5 cm away from the lower-left corner and 4.0 cm away from the upper-right corner, respectively. This crack length is about 14.8 cm, and its maximum width is about 1.7 mm. The damage variable of cohesive elements on this thermal-induced spalling crack ranges from 0.91 to 1.00.

When t = 1399 s, as shown in Figure 3A-3, the thermal deformations increase to 0.93–1.35 mm in the upper-left, lower-left and upper-right corners. The maximum displacement appears in the lower-right corner, and its value is about 3.0 mm. It can be seen from Figure 3B-3 that when the time reaches 1399 s, the high-temperature spalling damage of the concrete material is more severe. The numerical calculation has been unable to converge. It can be observed that many large thermal cracks have penetrated into the middle of the specimen. An obvious longitudinal thermal crack penetrates from the position 1.4 cm away from the upper boundary to the center of the specimen. This longitudinal thermal crack is about 7.8 cm long, and the maximum width in the middle part is about 0.4 mm. The damage variable on this thermal crack ranges from 0.80 to 1.00. It shows that the cohesive elements on this damage zone have been destroyed. These destroyed cohesive elements formed a thermal crack penetrating the specimen. At this time, the length of the largest thermal crack that appeared in the lower-right corner increases to 15.5 cm, and its maximum width is about 2.2 mm. The high temperature has caused the lower-right corner to be completely destroyed and spalled from the concrete specimen.

Figure 4A shows the temperature field in the concrete material at different times. It can be seen that the temperature continues to rise with increasing of time. Heat is gradually transferred from the outer boundaries to the interior of the concrete specimen. The temperature growth at the outer boundaries meets the predictions of the standard fire curve. The minor and major principal thermal stress fields are shown in Figures 4B,C. It is found that the thermal stresses change with increasing temperature. However, the overall distributions of the thermal stress fields do not change significantly with increasing temperature.

When t = 360 s, as shown in Figure 4A-1, the temperature at the outer boundaries of the model reaches 605.4°C. The temperature in the central area is slightly higher than the initial temperature of 20°C. The corresponding minor and major principal thermal stresses are shown in Figures 4B-1,C-1. The thermal stresses caused by high temperature are in a tensile state. Large aggregates greatly influence the thermal stress fields. In the areas 0.5–2.0 cm away from the boundaries, the temperature ranges from 210 to 370°C, where some large aggregates exhibit thermal stress concentrations. The minor principal thermal stresses are in a compressive state, ranging from −25.34 MPa to −30.87 MPa. The major principal thermal stresses are in a tensile state, ranging from 19.91 to 25.49 MPa.

As shown in Figure 4A-2, when t = 720 s, the temperature reaches 705.5 and 31.63°C at the outer boundaries and in the central area, respectively. The temperature field is roughly centrosymmetric. The aggregate distribution has little effect on the temperature field. Figures 4B-2,C-2 show that, at the same time, the thermal stresses are in a tensile state. At a distance of 0.5–2.0 cm away from the boundaries, the temperature rises to a range of 290–510°C. It is observed that thermal stresses concentrate around some large aggregates. Compared with t = 360 s, the values of major principal thermal stresses increase to a range from −25.81 MPa to −31.33 MPa. However, on the same aggregates, the values of major principal thermal stresses decrease to a certain extent due to the initiation and propagation of thermal cracks.

Figure 4A-3 shows that when t = 1080 s, the temperature at the boundaries increases to 964.9°C, and the temperature in the central area increases to 65.3°C. By referring to Figure 3, it is found that the initiation and propagation of a large number of thermal cracks cause the temperature field at t = 1080 s to be no longer symmetrical. There is an obvious temperature difference on both sides of the thermal crack in the lower-right corner. This indicates that large thermal cracks weaken the heat conductance towards the central region. The thermal stress fields at this time are shown in Figures 4B-3,C-3. It can be seen that, compared with t = 360 s and t = 720 s, the minor and major principal thermal stresses are still in compressive and tensile states, respectively. However, their values are reduced to a certain extent. This indicates that the initiation and propagation of thermal-induced spalling cracks lead to the damage of concrete material and then reduce the values of thermal stresses.

Figure 4A-4 shows that t = 1399 s, the concrete specimen is destroyed. On the boundaries and in the central areas, the temperature rises to 804.1 and 98.7°C, respectively. Compared with t = 1080 s, the high-temperature spalling failure is more severe. By referring to Figure 3, it is found that many thermal-induced spalling cracks propagate to the center of the specimen, which greatly affects the temperature distribution. Figure 4B-2,C-2 show that the propagation of thermal cracks reduces the values of minor and major principal thermal stresses. However, thermal-induced explosive spalling damage is relatively lower in the middle area, about 0.5–1.5 cm away from the left boundary. Some large aggregates in this area still exhibit a certain degree of thermal stress concentration.