Stochastic Fracture Analysis Using Scaled Boundary Finite Element Methods Accelerated by Proper Orthogonal Decomposition and Radial Basis Functions

This paper presents a stochastic analysis method for linear elastic fracture mechanics using the Monte Carlo simulations (MCs) and the scaled boundary finite element method (SBFEM) based on proper orthogonal decomposition (POD) and radial basis functions (RBF). The semianalytical solutions obtained by the SBFEM enable us to capture the stress intensity factors (SIFs) easily and accurately. The adoption of POD and RBF significantly reduces the model order and increases computation efficiency, while maintaining the versatility and accuracy of MCs. Numerical examples of cracks in homogeneous and bimaterial plates are provided to demonstrate the effectiveness and reliability of the proposed method, where the crack inclination angles are set as uncertain variables. It is also found that the larger the scale of the problem, the more advantageous the proposed method is.


Introduction
Fracture mechanics is central to many engineering applications, but the simulation of fractures poses great challenges to finite element analysis. To address this problem, a number of numerical algorithms were developed, such as meshfree methods [1,2], extended finite element methods [3][4][5], phase field methods [6], boundary element methods [7] and peridynamics methods [8]. Song and Wolf proposed a new semianalytical numerical method, the so-called scaled boundary finite element method (SBFEM) [9,10], which combines the advantages of the finite element method (FEM) and boundary element method (BEM). For the elasticity with simple geometries in which the entire boundary of the geometry can be seen from a scaling center, only the boundary needs to be discretized like BEM. For more complex problems, the interior domain can be divided into subdomains straightforwardly to meet the scaling center requirement. With SBFEM, the singularity of crack tips can be characterized with high accuracy, and the fracture parameters such as stress intensity factors can be directly obtained according to their definitions. In the past two decades, SBFEM has developed rapidly. Song and Wolf [11,12] studied the singular stress field at cracks in anisotropic composite materials. The crack propagation is simulated in [13][14][15][16][17] based on SBFEM. Zhang et al. [18] used SBFEM to solve the crack face contact problem and verified the effectiveness of the method. Tian et al. [19] employ SBFEM to calculate dynamic stress intensity factor and T-stress based on the improved continued fraction formula. Jiang et al. [20] proposed a method of using SBFEM to simulate strong and weak discontinuities by incorporating enrichment methods. Because of the advantages of isogeometric analysis in integrating computer-aided design and numerical analysis [21][22][23][24][25], Natarajan et al. [26] introduced isogeometric analysis to SBFEM for linear fracture mechanics where the physical fields are discretized with the spline basis functions of computer-aided design.
Uncertain problems are ubiquitous in practical engineering applications, and it is of paramount importance to evaluate the responses of the uncertain factors to ensure a reliable analysis result. Since deterministic analysis cannot characterize random fields, stochastic analysis has been extensively studied in recent years [27]. Stochastic analysis techniques can be divided into three categories: Karhunen-Loeve (KL) expanded stochastic spectrum methods [28], perturbation methods (PM) [29,30], and Monte Carlo simulation (MCs) [31][32][33]. As a direct sampling technique, MCs is considered to be a versatile method for uncertainty qualification with the merit of easy implementation [34,35]. Chowdhury et al. employ MCs and SBFEM to investigate the influence of crack geometric uncertainty [36,37] by taking crack lengths and angles as uncertain variables. On the other hand, MCs requires a large number of data samples to achieve high accuracy, leading to a huge computational cost. Alleviating the computational burden is critical for improving the performance of MCs. The model order reduction method based on the combination of proper orthogonal decomposition (POD) and radial basis functions (RBFs) provides an effective way to accelerate MCs [38][39][40][41]. POD can obtain the bases for the initial problems after they are projected to reduced spaces, thereby reducing the degrees of freedom of the governing equations. Ding et al. [42][43][44] coupled POD and MCs to conduct multidimensional uncertainty analysis and verified the effectiveness and efficiency of the method. RBF is used for continuous approximation of the system response that enables fast evaluation of the coefficients of the interpolation in the reduced space [45].
This paper proposes a novel procedure for stochastic analysis of linear fracture mechanics. With this method, SBFEM is used for linear elastic fracture analysis and MCs for uncertainty qualification, which is accelerated by the combination of POD and RBF. The remainder of this paper is structured as follows. Section 2 outlines the MCs in stochastic analysis. Section 3 introduces POD and RBF and how to integrate them with MCs. Section 4 outlines the SBFEM formula in linear elasticity fracture mechanics and the evaluation of stress intensity factors with SBFEM. Section 5 gives numerical examples to test the reliability and accuracy of the proposed method, followed by the conclusions in Section 6.

Monte Carlo Model
Monte Carlo simulation (MCs) is considered to be the simplest, more direct, and the main general-purpose tool in the study of uncertain problems. It is an experiment-based method to study the uncertain output caused by uncertainty of input parameters. Normally, expectations and standard deviations are used to describe this uncertainty. Due to the nature of expectation and variance, it is an integral problem. The Monte Carlo integral can be expressed as where gðxÞ is any integrable function and Ω ⊆ R d is the integration area.x is thed-dimensional random variable in the integration area Ω, which can be expressed as Suppose a set of uniformly distributed and independent random variablesfx i gin interval½a, bwith the density functionf X ðxÞ. Similarly, g * ðxÞis also a random variable with the same distribution and independent of each other, for anyg * ðxÞ that hasg * ðxÞ = g * ðxÞ/f X ðxÞ. So, the expectation of g * ðxÞ can be expressed as According to the theorem of large numbers, when n ⟶ ∞, the mean value of the mutually independent and identically distributed random variable sequence fx i g converges to the mean value of the probability distribution 1, i.e., We define the average value I as follows: By (3) and (4), it can be known that I converges to probability 1 of I, while I can provide an approximation of the expectation of solving I.
For more general problems, supposing that the problem domain is controlled by the single variable element α i and that the probability density function of each individual random variable is Pðα i Þ, then their joint probability density PðαÞ can be expressed as λðαÞ stands for the structural response of random variable α, and the mathematical expectation of any function f related to λðαÞ can be represented as where dV = dα 1 ⋅ dα 2 ⋯ dα n is an infinite small volume element; the expectation of function f ðλðαÞÞ can be obtained by Equation (4), i.e., where M is the number of sampling points; formula (7) can know that the prediction accuracy of Monte Carlo simulation depends on the number of sampling points. The more sampling points, the higher the prediction accuracy, but as the number of sample points increases, it will result in a large amount of computation. In order to solve this problem, consider using singular value decomposition (SVD) and 2 Geofluids radial base function (RBF) method to accelerate the stochastic analysis of MCs.

Proper Orthogonal Decomposition and Radial Basis Functions
It can be seen from the analysis above that direct Monte Carlo simulations require the consideration of computational efficiency. In order to improve computational efficiency, MCs has to be combined with other discrete algorithms. The proper orthogonal decomposition (POD) is usually used to improve computational efficiency. Following the idea of POD, the first step is to introduce design variable α i and its response λðα i Þ as a matrix, which can obtain a solution space Λ after a series of high-fidelity calculations; the solution space Λ can be represented as follows: where Λ ∈ R n×m , n is the number of response functions under any random input variable, and m is the number of design variables; that is, the number of sample points. In Equation (8) matrix, themandnsubscript inλ n,m represents thenth response of the structure under themth design variable, respectively. By the singular value decomposition method (SVD), the solution spaces Λ matrix can be decomposed into where r = min ðm, nÞ; U ∈ R n×n and V ∈ R n×n are all orthogonal matrices; u ij and v ij are the elements in the corresponding orthogonal matrices, respectively. In addition, u j and v j are the eigenvectors of ΛΛ T and Λ T Λ, respectively, which can also be called the left and right singular eigenvectors of matrix Λ; Σ ∈ R n×m is a diagonal matrix in which the diagonal element σ j is arranged in a descending order from the largest to the smallest. By defining φ j = u j and A j ðα i Þ, we obtain where φ j and A j ðα i Þ are the orthogonal basis and the corresponding amplitude, respectively. According to Equation (10), the system response can be simplified and expressed by the linear combination of φ j and A j ðα i Þ, and the degree of freedom of the simplified model is much smaller than that of the initial model. In addition, we use radial basis function (RBF) to interpole the subspace A j ðα i Þ. And we can get an approximation of the system response under any parameter. The approximate expression is as follows: where ϕ is the base function, φ i ðαÞ is taken as the Gauss kernel function, η is the coefficient corresponding to the base function, and m is the number of design variables. The φ i ðαÞ expression is as follows: in which the symbol k⋅k represents the Euclidean norm and the coefficient γ i is the width parameter of the function. Let AðαÞ ≈ÃðαÞ; Equation (11) can be rewritten as the following linear equations: According to the above linear equation system, we can get the coefficient η, and by substituting Equation (11) into Equation (10), we obtain In summary, the system response under any stochastic variable can be approximated by solving the above linear equation system without the need to use a full model to calculate. Using this feature, we can use a small part of the sample data to form a matrix that approximates the full-order system response. It avoids the iteration of complex physical algorithms and quickly obtains the response of variables based on MCs.

The Scaled Boundary Finite Element Method (SBFEM)
The scaled boundary finite element method (SBFEM) is a new semianalytical numerical calculation method. When performing numerical analysis, SBFEM is similar to FEM. Their basic idea is to discretize the calculation area into one or more subdomains and then follow a certain order. Assemble the subdomains to form an overall solution equation to solve the unknown parameters. As shown in Figure 1(a), an ordinary SBFEM arbitrary polygon unit is described. A so-called scaling center O is chosen in a region from which the total boundary S is visible. As shown in Figure 1(b), in the unit containing cracks, in addition to 3 Geofluids satisfying the condition of being visible to all boundaries, the scaling center must be placed at the crack tip, and the crack surface is described by two "side faces."

Scaled Boundary Finite Element Formulation.
During the derivation of the SBFEM equilibrium equation, the local coordinate system ðξ, sÞ should be established at the scaling center of the element, which is located at the crack tip O ðx 0 , y 0 Þ. The local coordinate system for a 2-D problem is defined in the local circumferential direction (s axis) and radial direction (ξ axis). The radial coordinate ξ is selected as ξ = 0 at the scaling center and ξ = 1 on the boundary. The coordinates ξ and s are called the scaled boundary coordinates. As shown in Figure 1, the coordinates of any point about Cartesian system on the boundary of the element can be obtained by interpolation using the local coordinates of the element: where x h and y h are the node coordinates of the element and ½NðsÞ is the shape function matrix, that is, ½NðsÞ = ½N 1 ðsÞ N 2 ðsÞ ⋯ N n ðsÞ. The coordinates of any point inside the element can be obtained by scaling the corresponding coordinates on the boundary in the radial direction ξ. Therefore, the coordinates of any point inside the element can be represented by the local coordinate system ðξ, sÞ as where ðx 0 , y 0 Þ is the coordinate of the scaling center O. Similarly, if the polar coordinate system ðr, θÞ is established at the scaling center, then the polar coordinate r and θ can be, respectively, expressed as In SBFEM, the displacement field of S element is expressed in a semianalytical way. Then, in the Cartesian coordinate system, the displacement of point ðξ, sÞ in the S element can be expressed as As can be seen from Equation (18), when the displacement of any point in the region is expressed, SBFEM only adopts interpolation function in the circumferential direction and is analytical in the radial direction. This is different from FEM which adopts interpolation function in both directions. Therefore, SBFEM is a semianalytical numerical method. The strain field and stress field under linear elastic conditions can be expressed as where ½D is the elasticity matrix and ½B 1 ðsÞ and ½B 2 ðsÞ describe the strain-displacement relationship [10,46]. At present, three different theories can be used to derive the governing equation of SBFEM, which are, respectively, Galerkin's weighted residual method [46], similarity principle [47], and virtual work principle [48]. The SBFEM static equilibrium equation is derived by using Galerkin's weighted residual method. The static equilibrium equation of SBFEM in two dimensions can be expressed as   4 Geofluids where ½E 0 , ½E 1 , and ½E 2 are coefficient matrices. It can be proved that ½E 0 is a positive definite matrix, ½E 1 is a positive semidefinite matrix, and ½E 2 is an asymmetric matrix. The internal node forces along radial lines can be expressed as [49] q By solving the static equation of scaled boundary finite element, the displacement field and stress field in each subdomain are, respectively, expressed as where and ½V 11 are eigenvalue matrices and eigenvector matrices obtained by block diagonal Schur decomposition in the process of solving equations. Meanwhile, ½S 11 contains all nonpositive eigenvalues, ½V 11 is the corresponding eigenvector, and fcg are integration constants.

Computation of Generalized Stress Intensity Factors
Using SBFEM. The definition of generalized stress intensity factors (SIFs) is not only to eliminate the stress singularity in the asymptotic solution of the crack tip but also to characterize the strength of the stress field near the crack tip. Therefore, generalized stress intensity factor is a very important physical quantity, which is often used as one of the important conditions for judging whether the crack is propagating. Otherwise, even the smallest external load can make its value tend to infinity, which is obviously not suitable for the criteria and basis of structural design. Generally speaking, the asymptotic solution is obtained from the eigenvalue problem. The eigenvalue and eigenfunction, respectively, correspond to the singular order and angle change of the stress field. Whether they are real or complex numbers is determined by the material properties and geometric characteristics. Because the mathematical functions of singular real numbers and complex numbers are different, the traditional definitions of generalized stress intensity factor in these two forms are independent of each other. In SBFEM, a polygon element as shown in Figure 1 can be used to simulate the crack, and only the boundary and crack surface need not be discretized.
In order to explain that the stress field obtained by SBFEM can explicitly characterize the singularity of the crack tip, eigenvalue matrix ½S 11 and eigenvalue vector matrix ½V 11 in Equation (23) are further divided into blocks as [50] where the eigenvalues of ½S ðsÞ (superscript (s) for singular) and ½S n1 need to satisfy −1 < λð½S ðsÞ Þ < 0 and λð½S n1 Þ < −1, respectively. If −1 < λð½S 11 Þ < 0 and ξ tends to 0, the value of matrix function ξ −½S 11 −½I in Equation (23) tends to infinity; then, singularity exists. Therefore, ½S ðsÞ and its corresponding eigenvector ½ψ ðsÞ can be directly used to express the stress singularity of the crack tip [51]. Therefore, according to the partition of eigenvalue matrix in Equation (24), Equation (23) can be rewritten as where where ½ψ ðsÞ σ ðsÞ is stress singular term, ½ψ ðTÞ σ ðsÞ is T stress term, and fc ðsÞ g and fc ðTÞ g are integration constants.
For convenience in fracture analysis, the stress modes are transformed to polar coordinates ðr, θÞ. By introducing the polar coordinate system ðr, θÞ at the crack tip as shown in Figure 2, the definition of the generalized stress intensity factor is further proposed. The original scale boundary coordinate ξ can be expressed as [50] where rðθÞ is the distance from the center of the scale to the boundary along the radial coordinate at the angle θ. The purpose of introducing the characteristic length L is to put forward the definition of the stress intensity factor independent of the unit. In Equation (25), the matrix power function of ξ is rewritten in the polar coordinates as [49] By substituting Equation (28) into Equation (26), we obtain where fσ ðsÞ ðr∧, θÞg and ½ψ ðsÞ σL ðθÞ are singular stress in polar coordinates and stress modes at the characteristic length L, respectively. At present, there are two definitions of stress intensity factor commonly used.
When materials 1 and 2 are the same homogeneous material as shown in Figure 2, the stress intensity factor of crack can be defined as whereK Ι and K II are stress intensity factors and σ ðsÞ θθ ðr, 0Þ and σ ðsÞ rθ ðr, 0Þ are two stress components in the front of the crack tip θ = 0. In two different isotropic materials, the cracks are distributed at the interface of the two materials; the stress intensity factor can be defined as c r∧ ð Þ = cos ε ln r∧ L , where L is the characteristic length. In this definition, the characteristic length L can be taken to any value, and the amplitude of the stress intensity factor does not change with the change of L. In general, L is the length of the crack. ε is the oscillation coefficient of the bimaterials, whose magnitude is only determined by the material parameters of the two materials. According to Equations (31) and (32), Equation (29) can be written as where fKðθÞgis the definition of generalized stress intensity factor. In the process of solving the stress intensity factor and simulating the crack growth, the stress value of the intersection points along the crack surface and the boundary is generally used to calculate the stress intensity factor, where θ = 0. In order to explain the relationship between the generalized stress intensity factor and the traditional stress intensity factor, auxiliary variables are introduced: Considering Equation (36), Equation (34) can be rewritten as In the homogeneous material plate with cracks, the singularity order obtained by SBFEM is ½S ðSÞ = 0:5½I, and ½I is the identity matrix. By substituting singularity order into Equation (37), take point A as an example, we obtain the classical definition: Note that the characteristic length L vanishes from the definition.
In an isotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is According to the matrix change, we can obtain Substituting Equation (40) into Equation (37), take point A as an example, we can obtain  6 Geofluids In an anisotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is 0:5 ± iε [50]. Cho et al. [52] gave the formula for calculating the generalized stress intensity factor. However, there was an error in the original formula, which was corrected by Song et al. [49]. The specific formula is as follows: where W 1 and W 2 can be determined by the material constant of the anisotropic plate. Considering Equation (34), Equation (42) can be rewritten as Therefore, the generalized stress intensity factor is unified with the traditional stress intensity factor, and there is no requirement for special modification when used, and when dealing with these three types of singular problems, a unified expression can be used for calculation, making the solution process more concise and clear.

Angled Crack in Rectangular Orthotropic Body.
A homogeneous plate with an angled crack is considered in this section, as shown in Figure 3(a). The homogeneous plate is represented by rectangular orthotropic body with the same material properties. The dimensions are chosen as b/h = 1 and a/h = 0:5. Length of crack is2a = 2; edge length is 2b = 2 h = 4. P = 1 is uniformly applied in the vertical direction. The elastic properties of the material are as follows: E 11 = E 22 = E 33 = 15:4 × 10 6 psi, G 11 = G 22 = G 33 = 15:7 × 10 6 psi, and ν 12 = ν 23 = ν 13 = 0:4009. Here, the structure is considered as a plane strain model. The plate is divided into two subdomains, each of which contains a crack tip, as shown in Figure 3(b). Each side of the subdomain is divided into 14 high-order (linear) elements with Gauss-Lobatto-Legendre shape functions, and the overall DOFs is 280. Herein, take the crack angle α ∈ ½0, π/2 as the uncertain parameter.
The numerical solution of this paper the SIFs is compared with the reference solution obtained by Banks-Sills et al. [53] using M-integral and displacement extrapolations, as shown in Figure 4, where they are normalized with P ffiffiffiffiffi ffi πa p . It can be seen from Figure 4 that the numerical solution based on SBFEM has a good agreement with that of Banks-Sills et al. [53], which verifies the correctness and effectiveness of solving SIFs in fracture problem based on SBFEM. Meanwhile, the relationship between SIFs and crack angle with different order elements is also considered, as shown in Figure 5. It can be seen that with the increase of crack angle, the value of normalized K 1 gradually decreases, and that of K 2 increases and then gradually decreases. However, when the crack angle is the same, the values of SIFs under different order elements are very close to almost the same, indicating that the order element has little effect on SIFs.
Since the mesh generation of SBFEM is more convenient and flexible than that of the FEM, only the position of the scaling center at the crack tip needs to be changed in this example without needing remeshing procedure.
We randomly selected 40 samples as the total number of input variables. In addition, we, respectively, construct several of small-scale response matrices containing 10, 20, and 30 samples from the total number of samples to form a low-order sample space vector. Then, perform SVD decomposition and RBF interpolation approximate calculation on the low-order matrix Λ formed by the three sets of small sample data. Figure 6 shows the normalized values of SIFs, which consist of original solution calculated using SBFEM and approximate solution using POD and RBF.
From Figure 6, it can be seen that when the sample numbers of a reduced model are 20, 30, and 40, the results of the reduced model are very close to the original solution, but when the sample number is 10, the result fluctuates slightly. It indicates that the combination of POD and RBF can accurately predict the response results of any input variable.
In order to further illustrate the reliability and accuracy of the algorithm. We introduce R 2 to evaluate the accuracy of the interpolation result and the formula is as follows: where n is the number of predictor variables when using RBF interpolation approximation and y i is the predicted value under the ith input variable; y∧ i is the analytical solution, which is approximated by the numerical solution; y i is the average value of the analytical solution. The evaluation coefficient R 2 represents the relative error between the numerical solution and the approximate solution, and its value range is (0-1). The closer R 2 is to 1, the higher the interpolation accuracy. Table 1 shows the comparison of the accuracy results of normalized K 1 and K 2 under four different schemes. In Table 1, the number of prediction points is 40, which means that the original 40 samples are selected as the prediction points. The number of prediction points is not equal to 40 7 Geofluids in the table, which means that some sample points are uniformly selected from the original 40 samples as new sample points, while the rest of the sample points are the number of predicted points. When the number of prediction points is 40, the accuracy evaluation coefficient R 2 is 1, which indicates that the approximate result of the RBF interpolation is consistent with the given numerical solution. As the number of samples decreases, the number of prediction points increases, and the value of the evaluation coefficient is gradually decreasing, which indicates that the interpolation accuracy of RBF is decreasing. However, when the number of samples is 10, the evaluation coefficient values of normalized K 1 and K 2 are both above 0.987, which can indicate that RBF can obtain high-precision numerical results by using small sample interpolation and verifies the reliability of the combination of POD and RBF to predict the response Disp. [53] Disp. [53] K 2 [50] (b) Figure 4: Normalized SIFs of angled crack in a homogeneous rectangular plate: (a) normalized SIFs K 1 /PðπaÞ 0:5 at different crack angles; (b) normalized SIFs K 2 /PðπaÞ 0:5 at different crack angles. 8 Geofluids results. Although the calculation accuracy will increase as the number of sample points increases, as the complexity of the problem increases, the calculation cost will also increase greatly. Therefore, it is necessary to select an appropriate number of sample points. Next, we will use POD and RBF to accelerate MCs to obtain the expectations and standard deviation of SIFs. Adopt the above physical model and take the crack angles α ðϵ ð0, π/2Þ as the uncertain parameter which satisfies the Gaussian distribution (μ = 0:785,σ = 0:262). The value of SIFs and the displacement values at 140 points are selected to form the response vector. Using 10, 20, 30, and 40 training data, respectively, we get 40 samples for random analysis  through the reduced-order fast algorithm and RBF interpolation. The dimensions of the different sample snapshot matrix are 142 × 10, 142 × 20, 142 × 30, and 142 × 40, which can be reduced by the SVD method in POD. Then, the fullorder 142 × 40 matrix is obtained by using RBF interpolation. In addition, before the RBF interpolation, we also truncated the decomposition matrix of SVD to varying degrees according to the size of the singular value. Different schemes are used as follows: case 1: direct Monte Carlo simulation is used to perform full-order simulation calculations using SBFEM; case 2: using SVD decomposition without truncating the singular value matrix after decomposition; and case 3: using SVD decomposition, the singular value matrix after decomposition is truncated according to the smallest row and column dimension. We compare the accuracy of statistical characteristics after different operations in Figure 7.
From Figure 7, we can see that the accuracy of the reduced-order fast algorithm increases with the increase of training samples which indicates that the accuracy of POD and RBF is directly related to the number of samples. Case 1 has 40 training samples, which can realize a full-level Monte Carlo simulation. But the number of samples is small, the deviation of cases 2, 3, and 1 are large. It is because the data space of cases 2 and 3 are smaller than that of case 1, where the degree of reduction of case 3 is the largest. When the number of samples is less than 20, cases 2 and 3 have the same calculation effect, but the dimension of case 3 is much smaller than that of case 2, so the prediction accuracy is higher. It indicates that this method can not only reduce the dimensionality of the data but also obtain the response results accurately and quickly. Besides, it is particularly important to select an appropriate number of samples for order reduction operations.  The plate is divided into four subdomains, as shown in Figure 8(b). The scaling center of the subdomains is indicated by the marker "⊕." The crack tips are at the scaling center of subdomains 1 and 2. Each side of the subdomain is divided into 10. We consider the main material angle θ ð1Þ ∈ ½0, π/2 to be the only uncertain parameter. The principal material axis 1 of material 1 is rotated from the  11 Geofluids horizontal x-axis by an angle θ ð1Þ to simulate anisotropic material behavior with reference to the crack. The step size of the design variable is set as 1, so we can get 90 samples. The two component values of SIFs and 158 nodal displacements values are taken as the response vector. The approximate effect of POD combined with RBF was tested here. Sampling points are selected according to step sizes 2, 5, and 10, respectively, to form 160 × 45, 160 × 18, and 160 × 9 matrices in sequence, and 160 × 90 matrix is obtained by SVD and RBF approximate interpolation. The comparison between interpolation result and numerical solution is shown in Figure 9.
From Figure 9, it can be seen that the results under different interpolation steps are in good agreement with the numerical solution, which indicates that the POD combined with the RBF method has a good approximation and prediction effect on the system response. Different steps correspond to different amounts of data, through interpolation computation, which can fit the full-order160 × 90matrix, which indicates that this method has a certain potential in accelerating computation. Table 2 shows the relative error between the approximate value and the numerical solution after POD order reduction and RBF interpolation under different steps. By comparison, it is found that the finer the interpolation steps, the closer the predicted value to the numerical solution, but the more interpolation points are needed. In general, the prediction accuracy is positively correlated with the number of interpolation points.
Next, we consider that when the input variable θ ð1Þ obeys the Gaussian distribution, the normalized K 1 expectation and standard deviation can be obtained after accelerating MCs through POD and RBF reduction. The mean value of Gaussian distribution function is μ = 45°, the standard deviation is σ = 15°, and the other parameters remain unchanged. The specific computation process is to stochastically select 100 sample points (satisfy the Gaussian distribution), obtain the numerical solution of K 1 based on SBFEM, and then directly use MCs to get expectations and standard deviation as a comparison. In addition, it can be seen from the previous example that the accuracy of the fast algorithm is positively correlated with the number of training samples. Therefore, in this example, we select 20, 40, 60, and 80  samples from the original samples as the training samples and then combine POD and RBF to get the expectations and standard deviations of these 100 sample points. It is worth noting that in this example, we all take case 3 reduced-order operation. The expectations and standard deviations of the full-order numerical solution and the reduced-order are shown in Figure 10. From Figure 10, it can be seen that the reliability of the reduced-order fast algorithm improves with the increase in the number of sample points. When the number of samples is close to 40% of the total number of samples, the results of the order reduction is consistent with the numerical solution. It shows that it is necessary to select the appropriate sample points for order reduction computation. Besides, it is also shown that the combination of POD and RBF is effective in model reduction and accelerating MCs computation, which achieves the goal of reducing computing cost and improving computing efficiency.

Conclusions
This paper presents an efficient stochastic analysis method for linear elastic fracture problems. The MCs is used to capture the statistical characteristics of the structural responses under the impact of uncertain variables. SBFEM is applied for fracture analysis which is able to evaluate the stress intensity factors directly with high accuracy. The model order reduction method based on proper orthogonal decomposition (POD) and radial basis function (RBF) is employed to accelerate Monte Carlo simulation (MCs). Numerical examples analyze the angled cracks in the orthotropic body and the interface crack structure between orthotropicisotropic materials. The results show that the full-order simulation computation results based on SBFEM are in good agreement with the results based on the POD and RBF model order reduction method, which verifies the effectiveness and efficiency of the proposed method. In the future, we will incorporate cohesive crack models for simulating fracture behaviors of quasibrittle materials that are commonly found in rock mechanics [54]. Additionally, the present method can also be used for accelerating structural optimization with large-scale design variables [55,56].

Data Availability
All date used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.