|Computer Modeling in Engineering & Sciences|
Monte Carlo Simulation of Fractures Using Isogeometric Boundary Element Methods Based on POD-RBF
1Key Laboratory of In-Situ Property-Improving Mining of Ministry of Education, Taiyuan University of Technology, Taiyuan, 030024, China
2School of Architecture and Civil Engineering, Huanghuai University, Zhumadian, 463003, China
3College of Architecture and Civil Engineering, Xinyang Normal University, Xinyang, 464000, China
4Artificial Intelligence Research Center, National Innovation Institute of Defense Technology, Beijing, 100071, China
5Department of Mechanical Engineering, Suzhou University of Science and Technology, Suzhou, 215009, China
*Corresponding Author: Zhongwang Wang. Email: email@example.com
Received: 25 March 2021; Accepted: 22 April 2021
Abstract: This paper presents a novel framework for stochastic analysis of linear elastic fracture problems. Monte Carlo simulation (MCs) is adopted to address the multi-dimensional uncertainties, whose computation cost is reduced by combination of Proper Orthogonal Decomposition (POD) and the Radial Basis Function (RBF). In order to avoid re-meshing and retain the geometric exactness, isogeometric boundary element method (IGABEM) is employed for simulation, in which the Non-Uniform Rational B-splines (NURBS) are employed for representing the crack surfaces and discretizing dual boundary integral equations. The stress intensity factors (SIFs) are extracted by M integral method. The numerical examples simulate several cracked structures with various uncertain parameters such as load effects, materials, geometric dimensions, and the results are verified by comparison with the analytical solutions.
Keywords: Monte Carlo simulation; POD; RBF; isogeometric boundary element method; fracture
Uncertainties are ubiquitous in engineering applications that may arise from different sources such as inherent material randomness, geometric dimensions, manufacturing errors, and dynamic loading. Because deterministic analysis fails to characterize randomness field, stochastic analysis techniques have been extensively studied to strengthen the credibility of computational prediction of uncertainty problems . There are three main variants of stochastic analysis: perturbation based techniques , stochastic spectral approaches [3,4], and Monte Carlo simulation (MCs) [5–7]. Among them, the MCs is regarded as the most versatile and simplest approach, and often used as the reference solution to verify the results of Perturbation method and Spectral method, but the exhaustive sampling in MCs leads to a heavy computational burden arising from both solving physical problems and constructing analysis-suitable geometric models that must be addressed carefully.
Combining Proper Orthogonal Decomposition (POD) and Radial Basis Functions (RBF) [8–11] is an effective technique of model order reduction [12–14]. The POD represents the solution field with the special ordered orthogonal functions in a low-dimensional subspace, which is constructed based on a discrete number of system responses obtained from the evaluation of the full order models (FOM) . POD reduces the degrees of freedom by capturing the dominant components of high-dimensional processes because it offers the optimal basis in the sense that the approximation error is minimal in L2 norm. On the other hand, the RBF builds a surrogate model through interpolating the data in the reduced space, whereby it admits continuous approximation of system responses for any arbitrary combination of input parameters  and thus does not need to solve partial differential equation for each sample.
The pre-processing time of MCs in constructing geometric models can be reduced with isogeometric analysis (IGA) . The key idea of IGA is employing spline functions used to construct geometric models in Computer Aided Design (CAD), for example Non-Uniform Rational B-splines (NURBS), T-splines , PHT-splines [19,20], and subdivision surfaces , as the basis functions to discretize physical fields. Compared to traditional Lagrange polynomial based methods, the main advantage of IGA lies in its ability of integrating numerical analysis and CAD. IGA enables one to perform numerical analysis directly from CAD models without meshing, which is particularly beneficial to uncertainty qualification since it requires fast generation of a large number of models. IGA also offers the benefits of geometric exactness, flexible refinement scheme and high order continuity that are amenable to numerical analysis.
Fracture mechanics is crucial in structural integrity assessment and damage tolerance analyses, but simulation of fracture behaviors poses significant challenges to Finite Element Methods (FEM) for the following reasons: (1) the mesh in proximity to cracks should be several orders of magnitude finer than that used for stress analysis; (2) a remeshing procedure is inevitable when cracks extend; (3) stress singularity or high stress gradients need to be captured. The extended finite element method (XFEM) incorporates enrichment functions to solution space and thus allows for crack propagation in a fixed mesh, but it still relies on a fine and good-quality mesh and necessitates special techniques to represent crack surfaces like Level Set Method. In comparison, Boundary Element Method (BEM) [22–26] has proven a useful tool for fracture simulation. Since BEM only discretizes the boundary of the domain, it not only reduces the degrees of freedom and facilities mesh generation, but more importantly, extends the crack by simply adding new elements at the front of cracks. In addition, as a semi-analytical method, BEM evaluates stress more accurately which is critical in extracting stress intensity factors. Furthermore, isogeometric analysis within the context of boundary element method (IGABEM) inherits the advantage of IGA in integration of CAD and numerical analysis and that of BEM in dimension reduction. IGA with BEM is also natural because both of them are boundary-represented. Since its inception, IGABEM has been successfully applied to potential [27–31], linear elasticity [32–36], acoustics [37–44], electromagnetics , structural optimization [40,46–49], etc., Peng et al. [50,51] applied IGABEM to two dimensional and three dimensional linear elasticity fracture mechanics, and demonstrates the accuracy and efficiency of IGABEM in this area.
This paper presents a novel procedure for solving uncertainty problems of fracture mechanics. In this method, the IGABEM is employed for fracture analysis, and MCs for addressing multiple random parameters. The computational cost is reduced by combination of POD and RBF. The remaining of this paper is structured as follows. Section 2 introduces the fundamentals of MCs in stochastic analysis. Section 3 illustrates how to apply POD and RBF to MCs. Section 4 formulates IGABEM in linear elasticity fracture mechanics. Several numerical examples are given in Section 5 to test the reliability, accuracy and efficiency of the proposed method, followed by conclusions in Section 6.
2 Stochastic Analysis with Monte Carlo Simulation
Monte Carlo simulation (MCs) directly characterizes uncertainties by calculating expectations and variances from a large number of samples. For a random variable X associated with the probability density function p(x), the two probabilistic moments are defined as
where E[X] is the expectation of the random variable X, and represents its variance, whose square root is the standard deviation.
According to the law of large numbers, the average of the results obtained from a number of samples should converge to the expectation as more sampling points are selected, which is the theoretical basis of MCs. Suppose g(X) is an arbitrary function of the random variable X. The expectation and variance of can be approximated by
where N is the sample size, and the order of convergence rate is O(N−1/2).
MCs in structural analysis can be conducted in the following steps : (1) Identify the random variables that are the source of the uncertainties in the system. (2) Determine the probability density functions of the random variables. (3) Use a random number generator to produce a set of samples which are adopted as the input parameters. (4) Employ a numerical method in deterministic analysis to evaluate the solution for each sample. (5) Based on the outputs of numerical analysis for all of the samples, we calculate the expectation and variance of the system using Eqs. (2) and (3) . From the above, it can be seen that MCs is easy to implement because the existing numerical simulation codes can be directly used without modification. In addition, MCs is versatile and is suitable for complex uncertainty problems. However, Step 3 of MCs is very time consuming because the numerical simulation needs to be conducted as many times as the number of sampling points. The large sample size can enhance the accuracy but may lead to higher computational cost.
3 Proper Orthogonal Decomposition (POD) and Radial Basis Functions (RBF)
As mentioned above, MCs is prohibitively expensive because it needs to solve the physical problems at many samples. This procedure can be accelerated by the reduced-order moeling based on POD and RBF. Let be an input random variable that has an influence on the structural responses. At the first step, we solve the FOM for a series of samples of the random variable. The system responses corresponding to the i-th sample point is called snapshots, which can be collected to form the following snapshot matrix :
where , n indicates the number of responses for any input variable, and m the number of samples. is the i-th response of the structure for the j-th sample. Decomposing the matrix through Singular Value Decomposition (SVD), we arrive at
where . and are orthogonal matrices whose entries are denoted by uij and vij, respectively. and are the eigenvectors of and , respectively. The two eigenvectors are also called the left and right singular matrices of matrix , respectively. is a diagonal matrix in which the diagonal elements (also called singular values) are arranged in descending order.
By defining and , Eq. (5) can be rewritten as
where is defined as the orthogonal basis, and is the corresponding amplitude. Using Eq. (6), the system responses at the selected samples can be expressed by a linear combination of and , which constitutes a reduced-order model (ROM) with less degrees of freedom than the FOM.
Eq. (6) only approximates a discrete number of system responses that are already computed using the FOM. To achieve a continuous approximation of system responses for any arbitrary input parameters, the radial basis functions (RBF) are used to interpolate the amplitudes in the reduced subspace
where N is the number of samples, is the i-th RBF and the corresponding coefficient. takes the form of Gaussian kernel function
in which the symbol represents the Euclidean norm, and the coefficient determines the width of the basis functions.
By letting , the values of the coefficient can be determined by solving Eq. (7). Substituting Eq. (7) into Eq. (6) gives rise to
Hence, the system response for any sample of the stochastic variable can be obtained straightforwardly by Eq. (9) without needing to solve the partial differential equations repeatedly. The algorithm mentioned above is illustrated in the Fig. 1. The POD-RBF enable us to conduct MCs without needing to perform FOM simulation for all of the sample points. However, the FOM is still essential for getting snapshots which is solved by IGABEM in our work as detailed in the following section. Therefore, the number of the samples or snapshots should be selected carefully to strike a balance between the accuracy and efficiency. In the future, we will introduce error estimation technique to improve the performance of our method. In addition, it is highlighted that in order to further enhance the computational efficiency, the number of basis functions in Eqs. (5) and (9) can be decreased by selecting the basis functions corresponding to the elements with larger values in .
4 Fracture Modeling with IGABEM
4.1 Boundary Integral Equations in Fracture Mechanics
Consider an arbitrary domain enclosed by boundary , where the portion of boundary is prescribed displacement boundary condition, is imposed traction boundary condition, and is a traction free crack with the upper surface and the lower surface (Fig. 2).
Because the crack upper surfaces and lower surfaces are geometrically overlapping, the source points on the two surfaces will coincide and thus the boundary integral equations (BIE) corresponding to them are identical to each other, which leads to the degradation of the system matrix. An effective approach—dual boundary element method (DBEM)  can overcome this difficulty by using the traction BIE on one of the crack surface (), and the displacement BIE on the other crack surface () and the rest of the boundary () (Fig. 2b) , i.e.,
where and are the source point and the field point, respectively. and are the fundamental solutions of displacement boundary integral equations, and D*ijk and S*ijk are fundamental solutions of traction boundary integral equations (Appendix A). denotes the jump term and equals 0.5 for smooth boundaries.
4.2 IGABEM for Fracture Mechanics
As a powerful geometric modeling technique, NURBS is the industrial standard of CAD and also central to isogeometric analysis. Given a knot vector , which is a non-decreasing set of coordinates in parametric space, a B-spline basis function is defined as follows:
where n denotes the order of basis functions and the parametric coordinate. As an extension of B-splines, NURBS basis function is written as
where is the NURBS basis function and wi the weight. A NURBS curve can be constructed by linear combination of NURBS basis functions and the associated coefficients
where stands for the Cartesian coordinate of a point and is its parametric coordinate. represents the Cartesian coordinates of control points in the physical space.
In IGABEM, NURBS are used not only for building geometries but also discretizing the boundary integral equations. Hence, the displacement and traction around the boundary can be expressed piecewisely as a linear combination of the NURBS basis function and the nodal parameters,
where and are the coefficients for discretizing the displacements and the tractions, respectively. It is noteworthy that they are not the nodal displacements or tractions because the NURBS basis functions lack the Kronecker-delta property.
By substituting Eq. (15) into the boundary integral equation (10), we can get the following discretization formulation of Eq. (10)
where represents the -th NURBS element, and is the local coordinate of the collocation point. In order to perform Gaussian quadrature, we map all the variables from physical space to local coordinate space , and the Jacobian transformation matrix is calculated as follows:
The boundary element method involves weakly singular, strongly singular and hyper singular integrals, which have to be addressed carefully. For this purpose, the subtraction of singularity technique is used, whose implementation details can be seen in . After solving the governing linear equations of IGABEM, we can obtain the displacement and traction field around the boundary. Then, we are able to evaluate the stress intensity factors (SIFs) with M integral  and predict the crack propagation direction with maximum hoop stress criterion . See Appendix B for details.
5 Numerical Examples
In the following examples, the input random variables are supposed to follow the Gaussian distribution with the expectation being 0.5 and the standard deviation being 0.033. The sampling method for Gaussian distribution function used in this paper is to select sample points in the interval . It is noted that only the problems with a single random input variable are considered in the present paper.The application of the method for the problems with multidimensional random input variables will be investigated in the future.
5.1 Inclined Center Crack Problem
A plate model with an inclined center crack under remote biaxial tension is considered in this section, as shown in Fig. 3. The crack inclination angle is , edge length L = 1, and crack length 2a = 0.02. The load is applied in the X2-direction and is applied in the X1-direction, in which is the load ratio that varies from 0 to 1 and . The Young’s modulus is E = 1 and Poisson’s ratio . The M integral method is used for evaluating SIFs in this example. The analytical solution of the SIFs for the inclined center crack problem is available, as follows :
The number of NURBS elements on the initial crack surface is set as 3, and the total DOFs is 86. We fix the crack angle to be , and set the load ratio as the input parameter. The response vector consists of 48 nodal displacements in X2-direction and two SIFs (K1 and K2). Before MCs is carried out we firstly test the accuracy of ROM in evaluating SIFs. We select 21 sample points of that are evenly distributed in the interval to get the reduced bases of POD and interpolate the system responses with RBF. It is found that the values of K1 and K2 computed by the numerical solution with ROM have a good agreement with the analytical solutions, as shown in Fig. 4.
Tab. 1 lists the normalized values of SIFs (the values divided by the exact solution) calculated by the FOM with IGABEM and that by the ROM with POD-RBF. The normalized value of K1 changes slightly with the increase of , while the normalized value of K2 increases steadily. Generally, both the results of the FOM and ROM are in good agreement with the analytical solutions. It is worth noting that the width of Gaussian kernel function has an important influence on the accuracy of RBF interpolation. According to the literature , the parameter can be determined via , where N represents the number of samples and m represents the dimension of the random variable.
To investigate the influence of the sample size on the accuracy of ROM, we construct a vector, , for different values of . From the vector, 20, 40, 60 and 80 samples are selected to form the reduced space, respectively, and the remaining points of the vector are used as the prediction points whose SIFs are computed with the ROM based on POD-RBF. We introduce R-squared method to quantify the accuracy of numerical results, as follows:
where n is the total number of the predication points, the subscript i the i-th prediction point, yi the numerical value of the SIF, the analytical solution of the SIF, and the average of the analytical solutions of the SIFs over these prediction points. The value of evaluation coefficient R2 is between 0 and 1, which indicates the deviation of numerical solutions and the exact solutions. The closer its value is to 1, the higher the accuracy of the numerical result is. As shown in Tab. 2, the evaluation coefficient R2 gradually increases with the sample number. When the sample number is 80, the evaluation coefficient R2 is 1.000, suggesting that the result of ROM is consistent with the analytical solution. It is also noted that R2 = 0.859 when the sample size is 20, so a good accuracy can still be reached with moderate sample size.
Next, we conduct MCs to evaluate the expectation and standard deviation of stress intensity factors. The following three different schemes are adopted for comparison.
1. Stochastic fracture analysis with MCs using FOM is conducted based on 51 samples of that are evenly distributed between 0 and 1. Correspondingly the system equations need to be solved with IGABEM 51 times.
2. We build the ROM with 10, 20, 30, and 40 samples, respectively, and then perform MCs using the 51 samples. In this case, the dimension of the matrix in Eq. (5) is selected as .
3. Similar to Case 2, the MCs is performed with the ROM that is established with 10, 20, 30 and 40 samples, respectively. However, we will further decrease the number of the reduced bases of the ROM by truncating the matrix by .
As can be seen from Figs. 5 and 6, when the number of samples for reduced-order modeling is small, the results of Cases 2 and 3 have large deviation from that of Case 1. However, with the increase of samples, the solutions of Cases 2 and 3 approach that of Case 1 rapidly. With 30 samples, the solutions of Cases 2 agree with Case 1 well, which demonstrates that the combination of POD and RBF can evaluate the expectation and standard deviation in stochastic analysis accurately. In addition, because the number of reduced bases in Case 3 is decreased, a larger error occurs compared to Case 2 although the computational time is further accelerated. Therefore, it is very important to select an appropriate order for the ROM based on POD-RBF to strike a balance between its accuracy and efficiency.
Now we fix the load ratio at 0.5, and study the influence of the crack inclination angle on the fracture analysis results. 50 samples of the input parameter are used for Monte Carlo simulation, in which 25 samples are used for ROM. The SIFs in terms of different input variables are evaluated by POD-RBF and FOM, respectively (Fig. 7). Tab. 3 lists the values of the expectation and standard deviation of SIFs computed by MCs with FOM and ROM, respectively. It can be found from this table that the results of ROM and the FOM have good agreement to each other, which verifies the correctness and effectiveness of the algorithm.
5.2 Cracks at Rivet Holes
In this section, we further test the performance of this algorithm using an example of a plate with two rivet holes, as shown in Fig. 8. The Young’s modulus is E = 1000 and Poisson’s ratio . The load is and . The crack inclination angle , initial crack length a = 0.1, and the number of NURBS elements on the crack surface is 3.
Suppose the radius of the right circular hole is a random variable and follows a Gaussian distribution. Based on the FOM simulation with IGABEM, Fig. 9 shows that the K1 of the right crack increases significantly with r. We select 20 samples to construct ROM, and conduct MCs with 40 sampling points. Tab. 4 shows the expectation, standard deviation and calculation time of MCs with FOM and ROM. It can be seen that the ROM leads to the results that are consistent with the FOM, and meanwhile reduces computational time.
Now we consider the influence of material parameter on displacement of cracks. The elastic modulus is set as the input variable, which follows Gaussian distribution. 51 samples are used for MCs, and 25 samples are used to construct the ROM based on POD and RBF. We choose the displacements of four collocation points on the crack surface as the response functions, whose coordinates are (2.91598244, 2.49410744), (2.85147368, 2.42959868), (2.08401756, 2.50589256), and (2.14852632, 2.57040132), respectively. Tab. 5 shows the expectation, relative standard deviation (RSD) of the displacement in y direction and calculation time corresponding to reduced-order and full-order MCs, respectively. It can be found from Tab. 5 that the stochastic analysis results of ROM are consistent with that of FOM.
We also study how the crack propagation path is influenced by the input parameter of . The crack propagation distance is fixed at 0.05 for each iteration step. A total of 16 steps are solved here. The path of crack propagation for is depicted in Fig. 10. Supposing the input parameter follows Gaussian distribution, we calculate the expectation and standard deviation of crack propagation path, which are represented by the coordinates of a series of crack tip positions. 10 samples are employed for reduced-order modeling and and 20 samples for MCs. It can be seen from Tab. 6 that the expectation and standard deviation of the crack extension path obtained from the ROM and that from the FOM have good consistency. It demonstrates the correctness and effectiveness of the stochastic analysis method proposed in this work.
The paper presents a novel framework for Monte Carlo simulation of multi-dimensional uncertainties of two-dimensional linear fracture mechanics. We use the NURBS to build the geometric model and discretize boundary integral equations. The IGABEM eliminates the repeatedly meshing procedure in uncertainty quantification, and retains the geometric exactness. The combination of POD and RBF accelerates the stochastic analysis and maintains good accuracy at low frequencies. In this work, the POD-RBF is used to approximate the structural response of a single random input variable. The method can also be extended to the problems with multidimensional random input variables, which will be investigated in the future. In addition, the present method will be applied to three dimensional problems and multi-physics coupling problems.
Funding Statement: The authors thank the financial support of National Natural Science Foundation of China (NSFC) under Grant (Nos. 51904202, 11902212, 11901578).
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
Appendix A: Fundamental Solution
The fundamental solutions of displacement integral equation and are
The fundamental solutions of traction integral equation D*ijk and S*ijk are
Appendix B: Evaluation of Stress Intensity Factors
The M integral based on the J integral is an effective method to extract SIFs. The Jk integral is defined as follows:
where Pkj is the Eshelby tensor, strain energy density . The center point of polar coordinate system is located at the crack tip. The limited integral boundary is an open circle with a center at the crack tip, as shown in Fig. A.1.
Now, consider two independent equilibrium states of an elastically deformable object, the actual state (superscript ‘1’) and the auxiliary state (superscript ‘2’). Superimpose these two equilibrium states into one equilibrium state. It is assumed that the crack surface is flat in a sufficiently small radius of the contour circle. The conservation law can be reduced to the path-independent J integral along an arbitrary path enclosing the crack tip located at the origin of the coordinates. Thus, the J integral of the superposition state can be expressed as
Rearranging Eq. (B.2) into the following formulas:
with W(1, 2) being the mutual potential energy density of the elastic body.
The J integral under two states can be expressed as a function associated with K, as follows:
where for plane strain problem, otherwise for plane stress problem.
Upon substitution of Eq. (B.3) into (B.5), we can obtain new expression of the M integral, as follows:
The M integral shown in Eqs. (B.4) and (B.6) only involves interaction terms, which can be directly used to solve the mixed mode crack problem of linear elastic solids.
After SIFs are calculated, the direction of crack propagation can be determined by the maximum hoop stress criterion, the discriminant equation is as follows:
where is the direction of crack propagation, in which the hoop stress is maximum. It can be seen that the accuracy of crack propagation simulation is mainly determined by K2/K1.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|