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

Authors: Xiaowei Shen; Haowen Hu; Zhongwang Wang; Xiuyun Chen; Chengbin Du
Journal: Geofluids
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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–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–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–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–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–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–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. Section2 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.


## 2. 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(1)q=∫Ωgxdx,where gx is any integrable function and Ω⊆Rd is the integration area.x is thed-dimensional random variable in the integration area Ω, which can be expressed as x=x1,x2,⋯xd. Suppose a set of uniformly distributed and independent random variablesxiin intervala,bwith the density functionfXx. Similarly, g∗xis also a random variable with the same distribution and independent of each other, for anyg∗x that hasg∗x=g∗x/fXx. So, the expectation of g∗x can be expressed as
(2)Eg∗x=∫abg∗xfXxdx=∫abgxdx=I.According to the theorem of large numbers, whenn⟶∞, the mean value of the mutually independent and identically distributed random variable sequence xi converges to the mean value of the probability distribution 1, i.e.,
(3)Prlimn⟶∞1n∑i=1ng∗xi=I=1.We define the average valueI¯ as follows:
(4)I¯=1n∑i=1ng∗xi.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
(5)Pα=Pα1Pα2⋯Pαn=∏i=1nPαi.λα stands for the structural response of random variable α, and the mathematical expectation of any function f related to λα can be represented as
(6)Efλα=∫ΩfλαpαdV,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.,
(7)Efλα≈1M∑i=1Mfλαi,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 radial base function (RBF) method to accelerate the stochastic analysis of MCs.


## 3. 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:
(8)Λ=λα1,λα2,⋯,λαm=λ11λ12⋯.λ1mλ21λ22⋯.λ2m⋯λn1λn2⋯.λnm,where Λ∈Rn×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,mrepresents 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
(9)Λ=UΣVT=∑j=1rujσjvjΤ,where r=minm,n; U∈Rn×n and V∈Rn×n are all orthogonal matrices; uij and vij are the elements in the corresponding orthogonal matrices, respectively. In addition, ujand vj are the eigenvectors of ΛΛT and ΛTΛ, respectively, which can also be called the left and right singular eigenvectors of matrix Λ; Σ∈Rn×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=uj and Ajαi, we obtain
(10)λαi=∑j=1rφjAjαi,where φj and Ajα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 Ajα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 Ajαi. And we can get an approximation of the system response under any parameter. The approximate expression is as follows:
(11)Aα≈A~α=∑i=1mηiϕiα,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:
(12)φiα=e−1/γi2α−αi,in which the symbol ⋅ represents the Euclidean norm and the coefficient γi is the width parameter of the function. Let Aα≈A~α; Equation (11) can be rewritten as the following linear equations:
(13)ϕ1α1ϕ2α1⋯ϕmα1ϕ1α2ϕ2α2⋯ϕmα2⋯⋯⋯⋯ϕ1αmϕ1αm⋯ϕmαmη1η2⋯ηm=Aα1Aα2⋯Aαm.According to the above linear equation system, we can get the coefficientη, and by substituting Equation (11) into Equation (10), we obtain
(14)λα=∑j=1rφjA~jα.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.


## 4. 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 Figure1(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 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.”Figure 1
Polygon representation of SBFEM: (a) polygonS-domain; (b) polygon S-domain with cracks.
(a)(b)


### 4.1. 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 Ox0,y0. 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:
(15)xs=Nsxh,ys=Nsyh,where xh and yh are the node coordinates of the element and Ns is the shape function matrix, that is, Ns=N1sNs2⋯Nns.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
(16)xξ,s=x0+ξxs=x0+ξNsxh,yξ,s=y0+ξys=y0+ξNsyh,where x0,y0 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
(17)rξ,s=ξrs=ξx2s+y2s,θs=arctanysxs.In SBFEM, the displacement field ofS 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
(18)u=uξ,s=Nsuξ.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
(19)ε=B1suξ,ξ+1ξB2suξ,σξ,s=DB1suξ,ξ+1ξB2suξ,where D is the elasticity matrix and B1s and B2s 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
(20)E0ξ2uξ,ξξ+E0−E1+E1Tξuξ,ξ−E2uξ=0,where E0, E1, and E2 are coefficient matrices. It can be proved that E0 is a positive definite matrix, E1 is a positive semidefinite matrix, and E2 is an asymmetric matrix. The internal node forces along radial lines can be expressed as [49]
(21)qξ=E0ξuξ,ξ+E1Tuξ.By solving the static equation of scaled boundary finite element, the displacement field and stress field in each subdomain are, respectively, expressed as(22)uξ,s=NsV11ξ−S11c,(23)σξ,s=−DB1sV11S11ξ−S11−Ic+DB2sV11ξ−S11−Ic,where and V11 are eigenvalue matrices and eigenvector matrices obtained by block diagonal Schur decomposition in the process of solving equations. Meanwhile, S11 contains all nonpositive eigenvalues, V11 is the corresponding eigenvector, and c are integration constants.


### 4.2. 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 Figure1 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 matrixS11 and eigenvalue vector matrix V11 in Equation (23) are further divided into blocks as [50]
(24)S11=diagSn1,−I,Ss,0,V11=ψn1,ψT,ψs,ψr,where the eigenvalues of Ss (superscript (s) for singular) and Sn1 need to satisfy −1<λSs<0 and λSn1<−1, respectively. If −1<λS11<0 and ξ tends to 0, the value of matrix function ξ−S11−I in Equation (23) tends to infinity; then, singularity exists. Therefore, Ss 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
(25)σξ,s=ψσssξ−Ss−Ics+ψσTscT+o1,where
(26)ψσss=D−B1sψsSs+B2sψs,ψσTs=DB1sψTSs+B2sψT,where ψσss is stress singular term, ψσTs is T stress term, and cs and cT are integration constants.For convenience in fracture analysis, the stress modes are transformed to polar coordinatesr,θ. 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]
(27)ξ=r∧rθ=Lrθ×r∧L,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]
(28)ξ−Ss−I=Lrθ−Ss−Ir∧L−Ss−I.Figure 2
Polygon representation of SBFEM [17].By substituting Equation (28) into Equation (26), we obtain
(29)σsr∧,θ=ψσLsθr∧L−Ss−Ics,(30)ψσLsθ=ψσssLrθ−Ss−I,where σsr∧,θ and ψσLsθ 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 Figure2, the stress intensity factor of crack can be defined as
(31)σθθsr,0τrθsr,0=12πrKΙKΙΙ,whereKΙ and KII are stress intensity factors and σθθsr,0 and σrθsr,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
(32)σθθsr∧,0τrθsr∧,0=12πr∧cr∧−sr∧sr∧cr∧KΙKΙΙ,(33)cr∧=cosεlnr∧L,sr∧=sinεlnr∧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
(34)σsr∧,θ=12πLψσLsθr∧L−Ss−IψσLsθ−1Kθ,(35)Kθ=2πLψσLsθcs,where Kθis 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:
(36)S~Sθ=ψσLsθSs+IψσLsθ−1.Considering Equation (36), Equation (34) can be rewritten as
(37)σsr∧,θ=12πLr∧L−S~SθKθ.In the homogeneous material plate with cracks, the singularity order obtained by SBFEM isS~S=0.5I, and I is the identity matrix. By substituting singularity order into Equation (37), take point A as an example, we obtain the classical definition:
(38)KIKII=2πrσθθsr,0τrθsr,0.Note that the characteristic lengthL vanishes from the definition.In an isotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is(39)S~S=0.5+ε−ε0.5.According to the matrix change, we can obtain(40)r∧L−S~S=r∧L−0.5I−0+ε−ε0=r∧L−0.5cr∧−sr∧sr∧cr∧.Substituting Equation (40) into Equation (37), take point A as an example, we can obtain
(41)KΙKΙΙ=2πr∧cr∧−sr∧sr∧cr∧−1σθθsr∧,0τrθsr∧,0.In an anisotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is0.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:
(42)σθθsr∧,0τrθsr∧,0=12πr∧W1W201cr∧−sr∧sr∧cr∧KIKII,where W1 and W2 can be determined by the material constant of the anisotropic plate. Considering Equation (34), Equation (42) can be rewritten as
(43)σθθsr∧,0τrθsr∧,0=12πLW1W201r∧L0.5+ε‐ε0.5W1W201−1KIKII.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.


## 4.1. 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 Ox0,y0. 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:
(15)xs=Nsxh,ys=Nsyh,where xh and yh are the node coordinates of the element and Ns is the shape function matrix, that is, Ns=N1sNs2⋯Nns.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
(16)xξ,s=x0+ξxs=x0+ξNsxh,yξ,s=y0+ξys=y0+ξNsyh,where x0,y0 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
(17)rξ,s=ξrs=ξx2s+y2s,θs=arctanysxs.In SBFEM, the displacement field ofS 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
(18)u=uξ,s=Nsuξ.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
(19)ε=B1suξ,ξ+1ξB2suξ,σξ,s=DB1suξ,ξ+1ξB2suξ,where D is the elasticity matrix and B1s and B2s 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
(20)E0ξ2uξ,ξξ+E0−E1+E1Tξuξ,ξ−E2uξ=0,where E0, E1, and E2 are coefficient matrices. It can be proved that E0 is a positive definite matrix, E1 is a positive semidefinite matrix, and E2 is an asymmetric matrix. The internal node forces along radial lines can be expressed as [49]
(21)qξ=E0ξuξ,ξ+E1Tuξ.By solving the static equation of scaled boundary finite element, the displacement field and stress field in each subdomain are, respectively, expressed as(22)uξ,s=NsV11ξ−S11c,(23)σξ,s=−DB1sV11S11ξ−S11−Ic+DB2sV11ξ−S11−Ic,where and V11 are eigenvalue matrices and eigenvector matrices obtained by block diagonal Schur decomposition in the process of solving equations. Meanwhile, S11 contains all nonpositive eigenvalues, V11 is the corresponding eigenvector, and c are integration constants.


## 4.2. 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 Figure1 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 matrixS11 and eigenvalue vector matrix V11 in Equation (23) are further divided into blocks as [50]
(24)S11=diagSn1,−I,Ss,0,V11=ψn1,ψT,ψs,ψr,where the eigenvalues of Ss (superscript (s) for singular) and Sn1 need to satisfy −1<λSs<0 and λSn1<−1, respectively. If −1<λS11<0 and ξ tends to 0, the value of matrix function ξ−S11−I in Equation (23) tends to infinity; then, singularity exists. Therefore, Ss 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
(25)σξ,s=ψσssξ−Ss−Ics+ψσTscT+o1,where
(26)ψσss=D−B1sψsSs+B2sψs,ψσTs=DB1sψTSs+B2sψT,where ψσss is stress singular term, ψσTs is T stress term, and cs and cT are integration constants.For convenience in fracture analysis, the stress modes are transformed to polar coordinatesr,θ. 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]
(27)ξ=r∧rθ=Lrθ×r∧L,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]
(28)ξ−Ss−I=Lrθ−Ss−Ir∧L−Ss−I.Figure 2
Polygon representation of SBFEM [17].By substituting Equation (28) into Equation (26), we obtain
(29)σsr∧,θ=ψσLsθr∧L−Ss−Ics,(30)ψσLsθ=ψσssLrθ−Ss−I,where σsr∧,θ and ψσLsθ 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 Figure2, the stress intensity factor of crack can be defined as
(31)σθθsr,0τrθsr,0=12πrKΙKΙΙ,whereKΙ and KII are stress intensity factors and σθθsr,0 and σrθsr,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
(32)σθθsr∧,0τrθsr∧,0=12πr∧cr∧−sr∧sr∧cr∧KΙKΙΙ,(33)cr∧=cosεlnr∧L,sr∧=sinεlnr∧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
(34)σsr∧,θ=12πLψσLsθr∧L−Ss−IψσLsθ−1Kθ,(35)Kθ=2πLψσLsθcs,where Kθis 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:
(36)S~Sθ=ψσLsθSs+IψσLsθ−1.Considering Equation (36), Equation (34) can be rewritten as
(37)σsr∧,θ=12πLr∧L−S~SθKθ.In the homogeneous material plate with cracks, the singularity order obtained by SBFEM isS~S=0.5I, and I is the identity matrix. By substituting singularity order into Equation (37), take point A as an example, we obtain the classical definition:
(38)KIKII=2πrσθθsr,0τrθsr,0.Note that the characteristic lengthL vanishes from the definition.In an isotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is(39)S~S=0.5+ε−ε0.5.According to the matrix change, we can obtain(40)r∧L−S~S=r∧L−0.5I−0+ε−ε0=r∧L−0.5cr∧−sr∧sr∧cr∧.Substituting Equation (40) into Equation (37), take point A as an example, we can obtain
(41)KΙKΙΙ=2πr∧cr∧−sr∧sr∧cr∧−1σθθsr∧,0τrθsr∧,0.In an anisotropic bimaterial plate with interface cracks, the singular order obtained by SBFEM is0.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:
(42)σθθsr∧,0τrθsr∧,0=12πr∧W1W201cr∧−sr∧sr∧cr∧KIKII,where W1 and W2 can be determined by the material constant of the anisotropic plate. Considering Equation (34), Equation (42) can be rewritten as
(43)σθθsr∧,0τrθsr∧,0=12πLW1W201r∧L0.5+ε‐ε0.5W1W201−1KIKII.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.


## 5. Numerical Examples

### 5.1. Angled Crack in Rectangular Orthotropic Body

A homogeneous plate with an angled crack is considered in this section, as shown in Figure3(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=2h=4. P=1 is uniformly applied in the vertical direction. The elastic properties of the material are as follows: E11=E22=E33=15.4×106psi, G11=G22=G33=15.7×106psi, 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.Figure 3
A rectangular plate with an angled crack [50]: (a) geometry; (b) mesh of 14th-order elements.
(a)(b)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πa. 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 K1 gradually decreases, and that of K2 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.Figure 4
Normalized SIFs of angled crack in a homogeneous rectangular plate: (a) normalized SIFsK1/Pπa0.5 at different crack angles; (b) normalized SIFs K2/Pπa0.5 at different crack angles.
(a)(b)Figure 5
The relationship between SIF and crack angle under different element orders. (a) The normalized SIFsK1/Pπa0.5 varies with the crack angle under different element orders. (b) The normalized SIFs K2/Pπa0.5 varies with the crack angle under different element orders.
(a)(b)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.Figure 6
Curve graph of normalized SIFs and numerical solution under different training points. (a) The normalized SIFsK1/Pπa0.5 varies with the crack angle under different training points. (b) The normalized SIFs K2/Pπa0.5 varies with the crack angle under different training points.
(a)(b)From Figure6, 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 introduceR2 to evaluate the accuracy of the interpolation result and the formula is as follows:
(44)R2=1‐∑i=1nyi−y∧i2∑i=1nyi−y¯i2,where n is the number of predictor variables when using RBF interpolation approximation and yi 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 R2 represents the relative error between the numerical solution and the approximate solution, and its value range is (0-1). The closer R2 is to 1, the higher the interpolation accuracy.Table1 shows the comparison of the accuracy results of normalized K1 and K2 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 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 R2 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 K1 and K2 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 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.Table 1
Comparison of the accuracy results of four different schemes.
Number of samplesNumber of prediction pointsR2 for normalized K1/Pπa0.5R2 for normalized K2/Pπa0.540401130101120201110300.99460.9875Next, 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 full-order 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.Figure 7
Expectations and standard deviation of different degrees of reduction: (a) expectations of the normalized SIFsK1/Pπa0.5; (b) expectations of the normalized SIFs K2/Pπa0.5; (c) standard deviation of the normalized SIFs K1/Pπa0.5; (d) standard deviation of the normalized SIFs K2/Pπa0.5.
(a)(b)(c)(d)From Figure7, 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.


### 5.2. Interfacial Central Crack between Isotropic-Orthotropic Materials

The rectangular bimaterial plate composed of isotropic-orthotropic materials is considered in this section, as shown in Figure8(a). The height and the width of the rectangular plate are 2H and H, respectively. The length of the central crack is 2a. Uniform tension P is applied in the vertical direction. The elastic properties are E11=100GPa, E22=30GPa, G12=23.5212GPa, and ν12=0.3 for orthotropic material 1 and E=100GPa and ν=0.3 for isotropic material 2. Plane stress conditions are assumed. The characteristic length is chosen as L=2a.Figure 8
Rectangular bimaterial plate with an interfacial central crack [50]: (a) geometry; (b) mesh of 10th-order elements.
(a)(b)The plate is divided into four subdomains, as shown in Figure8(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 horizontal x-axis by an angle θ1to 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.Figure 9
Different length interpolation results: (a) normalized SIFsK1/Pπa0.5 under different angles θ1; (b) normalized SIFs K2/Pπa0.5 under different angles θ1.
(a)(b)From Figure9, 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.Table 2
Relative error table of approximate value and numerical solution under different interpolation steps.
Angleθ1|RBF-exact|/exactNormalized SIFsK1/Pπa0.5Normalized SIFsK2/Pπa0.5Step=2Step=5Step=10Step=2Step=5Step=107°0.006%0.040%2.753%0.035%0.251%0.800%13°0.004%0.025%1.931%0.044%0.166%0.492%27°0.002%0.004%0.945%0.004%0.005%0.287%37°0.001%0.004%0.604%0.024%0.039%0.212%53°0.004%0.004%0.360%0.057%0.057%0.668%67°0.054%0.066%0.379%0.067%0.131%0.555%73°0.129%0.377%1.217%0.064%0.101%0.601%87°0.051%0.604%2.601%0.071%0.518%3.570%Next, we consider that when the input variableθ1 obeys the Gaussian distribution, the normalized K1 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 K1 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.Figure 10
Expectations and standard deviation of different degrees of reduction: (a) expectations of the normalized SIFsK1/Pπa0.5 under different sample points; (b) standard deviation of the normalized SIFs K1/Pπa0.5 under different sample points.
(a)(b)


## 5.1. Angled Crack in Rectangular Orthotropic Body

A homogeneous plate with an angled crack is considered in this section, as shown in Figure3(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=2h=4. P=1 is uniformly applied in the vertical direction. The elastic properties of the material are as follows: E11=E22=E33=15.4×106psi, G11=G22=G33=15.7×106psi, 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.Figure 3
A rectangular plate with an angled crack [50]: (a) geometry; (b) mesh of 14th-order elements.
(a)(b)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πa. 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 K1 gradually decreases, and that of K2 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.Figure 4
Normalized SIFs of angled crack in a homogeneous rectangular plate: (a) normalized SIFsK1/Pπa0.5 at different crack angles; (b) normalized SIFs K2/Pπa0.5 at different crack angles.
(a)(b)Figure 5
The relationship between SIF and crack angle under different element orders. (a) The normalized SIFsK1/Pπa0.5 varies with the crack angle under different element orders. (b) The normalized SIFs K2/Pπa0.5 varies with the crack angle under different element orders.
(a)(b)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.Figure 6
Curve graph of normalized SIFs and numerical solution under different training points. (a) The normalized SIFsK1/Pπa0.5 varies with the crack angle under different training points. (b) The normalized SIFs K2/Pπa0.5 varies with the crack angle under different training points.
(a)(b)From Figure6, 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 introduceR2 to evaluate the accuracy of the interpolation result and the formula is as follows:
(44)R2=1‐∑i=1nyi−y∧i2∑i=1nyi−y¯i2,where n is the number of predictor variables when using RBF interpolation approximation and yi 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 R2 represents the relative error between the numerical solution and the approximate solution, and its value range is (0-1). The closer R2 is to 1, the higher the interpolation accuracy.Table1 shows the comparison of the accuracy results of normalized K1 and K2 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 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 R2 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 K1 and K2 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 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.Table 1
Comparison of the accuracy results of four different schemes.
Number of samplesNumber of prediction pointsR2 for normalized K1/Pπa0.5R2 for normalized K2/Pπa0.540401130101120201110300.99460.9875Next, 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 full-order 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.Figure 7
Expectations and standard deviation of different degrees of reduction: (a) expectations of the normalized SIFsK1/Pπa0.5; (b) expectations of the normalized SIFs K2/Pπa0.5; (c) standard deviation of the normalized SIFs K1/Pπa0.5; (d) standard deviation of the normalized SIFs K2/Pπa0.5.
(a)(b)(c)(d)From Figure7, 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.


## 5.2. Interfacial Central Crack between Isotropic-Orthotropic Materials

The rectangular bimaterial plate composed of isotropic-orthotropic materials is considered in this section, as shown in Figure8(a). The height and the width of the rectangular plate are 2H and H, respectively. The length of the central crack is 2a. Uniform tension P is applied in the vertical direction. The elastic properties are E11=100GPa, E22=30GPa, G12=23.5212GPa, and ν12=0.3 for orthotropic material 1 and E=100GPa and ν=0.3 for isotropic material 2. Plane stress conditions are assumed. The characteristic length is chosen as L=2a.Figure 8
Rectangular bimaterial plate with an interfacial central crack [50]: (a) geometry; (b) mesh of 10th-order elements.
(a)(b)The plate is divided into four subdomains, as shown in Figure8(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 horizontal x-axis by an angle θ1to 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.Figure 9
Different length interpolation results: (a) normalized SIFsK1/Pπa0.5 under different angles θ1; (b) normalized SIFs K2/Pπa0.5 under different angles θ1.
(a)(b)From Figure9, 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.Table 2
Relative error table of approximate value and numerical solution under different interpolation steps.
Angleθ1|RBF-exact|/exactNormalized SIFsK1/Pπa0.5Normalized SIFsK2/Pπa0.5Step=2Step=5Step=10Step=2Step=5Step=107°0.006%0.040%2.753%0.035%0.251%0.800%13°0.004%0.025%1.931%0.044%0.166%0.492%27°0.002%0.004%0.945%0.004%0.005%0.287%37°0.001%0.004%0.604%0.024%0.039%0.212%53°0.004%0.004%0.360%0.057%0.057%0.668%67°0.054%0.066%0.379%0.067%0.131%0.555%73°0.129%0.377%1.217%0.064%0.101%0.601%87°0.051%0.604%2.601%0.071%0.518%3.570%Next, we consider that when the input variableθ1 obeys the Gaussian distribution, the normalized K1 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 K1 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.Figure 10
Expectations and standard deviation of different degrees of reduction: (a) expectations of the normalized SIFsK1/Pπa0.5 under different sample points; (b) standard deviation of the normalized SIFs K1/Pπa0.5 under different sample points.
(a)(b)


## 6. 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 orthotropic-isotropic 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].



---

*Source: 9181415-2021-09-06.xml*