In damage theory, Kachanov regarded microcracks and microvoids in materials as defects that do not bear load (Kachanov, 1992). Consequently, the effective load-bearing area of the material decreases from A to A ′ . If we assume that these microdefects are uniformly distributed in all spatial directions, and that A ′ is independent of the normal direction, an isotropic damage variable δ can be defined as in Equation 1. δ = A − A ′ / A (1) where δ is the isotropic damage variable, A is the nominal load-bearing area, and A ′ is the effective load-bearing area.

When studying structural damage in buildings, both the length and depth of cracks influence the load-bearing capacity and operational performance of the structure (Laxman et al., 2023; Jafarifar et al., 2016; Lin et al., 2012). Existing engineering practices and theoretical research indicate that cracks induced by disturbances from foundation pit excavation typically appear at the junctions of walls and top slabs. Since these components are relatively thin, the resulting macrocracks can be approximated as planes perpendicular to the direction of loading (Li et al., 2022). The area of these cracks can be approximated as the product of the crack length and depth. In this study, the cracked portion of the building is treated as an area that cannot bear tensile forces. Thus, the larger the crack area, the smaller the effective load-bearing area and the greater the degree of structural damage. Conversely, the smaller the crack area, the larger the effective load-bearing area and the lesser the degree of structural damage. Therefore, the crack propagation area c is used to represent the extent of damage in the building, and it is defined as in Equation 2. c = L × B (2) where c is the crack propagation area, L is the crack length, and B is the average crack depth.

In the XFEM model, a crack plane with a specified length and width can be set to be inserted at a specified position in the model as pre-damage of the structure. During XFEM analysis, the mesh that intersects the crack plane is considered to be cracked and cannot withstand the tensile force of the crack plane normal. In this way, the Kachanov damage theory can be applied to XFEM analysis.

A 3D finite element model that integrates the soil, foundation pit, and building was constructed via advanced finite element software (Figure 1). The model dimensions are 400 m × 150 m × 50 m, and a total of 60,384 elements are generated. The soil, diaphragm walls, and building are all modelled with solid elements. The contact between the soil and diaphragm walls, as well as between the building foundation and the soil, is set as surface-to-surface contact, with tangential behavior defined as “penalty friction” and normal behavior as “hard contact.” The diaphragm walls and internal supports are modelled with coupled contact conditions. The soil is represented via the Mohr‒Coulomb material model, whereas the internal supports and building structures are modelled via a linear elastic material model (Hu et al., 2021; Zhou et al., 2024; Deng et al., 2023). Boundary conditions constrain displacements in all three directions at the bottom of the soil and normal displacements along the sides. Dewatering of the foundation pit is simulated by setting pore pressure boundary conditions, and isotropic seepage models are used for groundwater flow (Luo and Li, 2019; Hu et al., 2024). The excavation of the foundation pit is simulated via the *Model Change feature, with excavation steps set according to the conditions listed in Table 1. The specific parameters of the FEM model are shown in Table 2.

Since conventional finite element methods cannot simulate strong discontinuities such as cracks, the extended finite element method (XFEM) is used to simulate crack propagation in buildings under settlement (Yang et al., 2024; Ávila et al., 2024). Referring to the two-stage analysis method in the literature (Liu et al., 2020; Liang et al., 2017; Liu et al., 2021; Korff et al., 2016; Shi et al., 2016), an XFEM model containing only the building was established (Figure 2). The settlement calculated by the FEM model at the building’s foundation is applied as displacement vectors to the corresponding positions in the XFEM model (Figure 3). The cohesive crack initiation criterion is used in the model, and on the basis of reference (Kaklauskas and Ghaboussi, 2001; Liu et al., 2023), the maximum principal stress is set to 1.7 MPa and the fracture energy to 150 N/m. The specific parameters are shown in Table 3.

While XFEM is an effective tool for modeling crack propagation in building structures under foundation pit excavation, it has certain limitations. For example, XFEM assumes that cracks propagate according to predefined failure criteria (such as stress intensity factors or strain criteria). This may not fully capture complex crack behaviors in real-world materials, where cracks can sometimes exhibit more irregular patterns or be influenced by additional factors like soil-structure interaction or dynamic loading. Additionally, XFEM requires a sufficiently refined mesh to accurately capture crack initiation and propagation, which can be computationally expensive. These limitations should be considered when applying XFEM in real-world scenarios, and further research may be needed to improve its applicability to a wider range of structural conditions.

In actual construction, different working conditions involve varying construction parameters, meaning that there are countless scenarios. If a numerical model were to be built for each scenario, it would undoubtedly increase computational costs. Therefore, this paper establishes a response surface model to predict the computational results under different construction parameters. The response surface method (RSM) is a statistical approach based on scientifically reasonable experimental designs, which quantifies variables and, through multiple trials, produces a simple mathematical model to describe the unknown relationships between these variables.

Given the high stiffness of the argillaceous siltstone layer, its variability has a relatively small effect on ground settlement and is therefore excluded from the parameter analysis. Additionally, since the plasticity parameters of sandy soil are generally fixed and show minimal variability, they are also excluded. Referring to the parameter sensitivity analysis by Lu (Lu, 2021), this study selects the following influencing factors: the elastic modulus E 1 , density ρ 1 , cohesion c 1 , internal friction angle φ 1 , and Poisson’s ratio μ 1 of the silty clay layer, as well as the elastic modulus E 2 , density ρ 2 , and Poisson’s ratio μ 2 of the sandy soil layer. Furthermore, the distance between the foundation pit and the building D is considered, resulting in a total of nine influencing factors, and the dependent variable is the crack propagation area c in the building. The first-order model formula is as follows Equation 3. c = k 0 + ∑ i = 1 9 k i x i + ε (3) where x i represents the i-th influencing factor ( i = 1 ∼ 9 ), k i denotes the coefficient of x i , and ε represents the error term.

The third-order model formula includes square and cubic terms, which allows consideration of the strong nonlinear influence of the factors on the predicted variable. The formula is as follows Equation 4. c = k 0 + ∑ i = 1 9 k i x i + ∑ i = 1 9 k i + 9 x i 2 + ∑ i = 1 9 k i + 18 x i 3 + ε (4)

In this equation, the nine independent variables lead to 28 coefficients, which requires 28 equations to solve. To eliminate the effects of dimensional differences between parameters and ensure the comparability of the data, the orthogonal test results were normalized via the Equation 5 (Han et al., 2020; Cao et al., 2014): x ′ = x − x min x max − x min (5) where x ′ represents the normalized data, x represents the original data, x max ​ represents the maximum value, and x min represents the minimum value.

When the response surface function is fit via the least squares method, high-order functions can easily lead to overfitting. Therefore, ridge regression was used to fit the re-sponse surface function in this study. Ridge regression, a regularized linear regression model, helps mitigate multicollinearity and overfitting issues. The penalty parameter for ridge regression was set to 1.

Since the number of influencing factors in this study is large, analysing the effect of each factor on crack propagation individually would require many experiments. To optimize computational resources, the orthogonal experimental method is used to investigate the impact of various factors on crack propagation. In orthogonal experiments, too few levels for each factor can affect the prediction accuracy, whereas too many levels can result in excessive computational demands. On the basis of the coefficient of variation for each parameter, four levels were selected for each factor. Ultimately, the L32(4⁹) orthogonal experiment table was chosen, and the resulting orthogonal table is shown in Table 4. In this table, Group 1 corresponds to the average values of all the parameters.

The normalized parameters calculated by Equation 5 are as follows Equation 6. E 1 ′ = 0.119 E 1 − 0.746 , E 2 ′ = 0.037 E 2 − 0.444 , ρ 1 ′ = 5.556 ρ 1 − 10.444 , ρ 2 ′ = 5.64 ρ 2 − 10.444 , c 1 ′ = 0.033 c 1 − 0.444 , φ 1 ′ = 0.176 φ 1 − 3 , μ 1 ′ = 11.111 μ 1 − 3 , μ 2 ′ = 11.111 μ 2 − 2.667 , D ′ = 0.059 D − 0.176 (6)

As shown in Figure 7, the crack propagation area at the top of the shear wall was extracted, and the results for each group are displayed in Figure 8. The analysis reveals the following: Figure 9.

(1) In groups 12 and 16, no crack propagation was observed. In contrast, Group 2 experienced the most severe crack propagation, with a crack length of 9.6 m and an average depth of 0.156 m. The other groups displayed varying degrees of damage propagation. Notably, in several cases, cracks extended to the top slab of the building, with the total crack length in Group 22 exceeding 14 m.

The response surface model was constructed based on the orthogonal test results, and the regression coefficients were calculated using ridge regression in MATLAB. Table 6 presents the regression coefficients obtained from the model (the intercept is 1.896), which describe the influence of each factor on crack propagation.

The comparison between the numerical simulation results of 32 working conditions and the values predicted from the response surface function is shown in Figure 10. In civil engineering and fracture mechanics, data typically exhibit high variability and uncertainty, and if the mean squared error (MSE) is less than 0.4, the fit is considered good (Güçlüer et al., 2021; Chai et al., 2023; Salami et al., 2021). In this study, the ridge regression MSE was 0.15, confirming the reliability of the predictions for crack propagation in existing buildings.

To verify the predictive effectiveness of the above model, four random working conditions, as shown in Table 7, were selected. Finite element calculations were performed to obtain the damage development for each condition (Figure 11), specifically the crack propagation area, and the results were compared with those predicted by the model, as shown in Figure 12. The analysis reveals that the damage trends predicted by the model closely match the numerical simulation results, indicating that the established response surface model can effectively predict damage propagation in buildings affected by foundation pit excavation.

The FEM-XFEM model provides detailed, high-accuracy predictions of crack initiation and propagation, but it is computationally expensive, particularly when simulating a large number of scenarios. In contrast, the RSM model offers a computationally efficient alternative, though it sacrifices some precision for speed. The comparison between the two models reveals that both models are effective, but the choice between them depends on the specific application—FEM-XFEM is preferable for detailed analysis, while RSM is more suitable for rapid predictions across multiple scenarios.

To analyse the sensitivity of building damage to various influencing factors, further calculations were conducted to determine the Pearson Product-Moment Correlation Coefficient (PPMCC) and the Maximal Information Coefficient (MIC) between the total crack propagation area and each parameter. The PPMCC measures the linear relationship between variables, with values ranging from −1 to 1, whereas the MIC measures the nonlinear relationship between variables, with values ranging from 0 to 1. The results are shown in Figure 13, and the analysis reveals the following:

(1) In terms of the MIC, the distance between the foundation pit and the building ( D ) has the strongest correlation (0.828), followed by the elastic modulus of the sandy soil layer ( E 2 , 0.596) and the Poisson’s ratio of the silty clay layer ( μ 1 , 0.419). These parameters significantly affect the crack propagation area in a nonlinear manner. In terms of linear correlation (PPMCC), the elastic modulus of the silty clay layer ( E 1 , −0.45) and the distance between the foundation pit and the building ( D , 0.27) show strong linear correlations, with both exhibiting a negative correlation.

To analyse the influence of each factor on building crack propagation individually, all other factors were set to their average values while the target factor was varied. For each target factor, standardized values were selected at intervals of 0.01. This study involves nine influencing factors, resulting in a total of 900 scenarios. The prediction model from Section 4.2 was used to calculate the crack propagation area for each scenario. Figure 14 illustrates the impact of each of the nine influencing factors on crack propagation. The analysis reveals:

(1) The trends for E 1 ′ , E 2 ′ , μ 1 ′ , and D ′ are quite similar, with the crack propagation area decreasing or stabilizing as the independent variable increases. Notably, E 1 ′ and E 2 ′ show the most significant impact, with the crack propagation area decreasing from approximately 1.6 m² to approximately 0.4 m². A higher elastic modulus means less soil deformation under the same construction load, resulting in a smaller impact on existing buildings. Additionally, D ′ indicated that attention should be given to crack development in buildings closer to foundation pits. On the other hand, c 1 ′ , ρ 2 ′ , and μ 2 ′ significantly increase in the crack propagation area as the factor increases, with a “turning point” ( c 1 ′ = 0.36 , ρ 2 ′ = 0.68 , μ 2 ′ = 0.31 ) where the crack area reaches its maximum.