Abstract
The extend finite element method (XFEM) is popular in structural mechanics when dealing with the problem of the cracked domains. XFEM ends up with a linear system. However, XFEM usually leads to nonsymmetric and ill-conditioned stiff matrix. In this paper, we take the linear elastostatics governing equations as the model problem. We propose a new iterative method to solve the linear equations. Here we separate two variables and Enr, so that we change the problems into solving the smaller scale equations iteratively. The new program can be easily applied. Finally, numerical examples show that the proposed method is more efficient than common methods; we compare the -error and the CPU time in whole process. Furthermore, the new XFEM can be applied and optimized in many other problems.
1. Introduction
Many physical phenomena can be modeled by partial differential equations with singularities and interfaces. Much attention has been focused on the development of interface problem in fluid dynamics and material science. The standard finite difference and finite element methods may not be successful in giving satisfactory numerical results for such problems. Hence, many new methods have been developed. Some of them are developed with the modifications in the standard methods, so that they can deal with the discontinuities and the singularities. Peskin developed immersed boundary method (IBM) [1] basically to model blood flow in the heart. Immersed interface method (IIM) [2–6] was basically developed because of the need for a second-order accurate finite difference method to solve linear elliptic, parabolic partial differential equations in which an integral equation may not be available. IIM will give a nonsymmetric system of equations. The extended finite element method (XFEM) [7, 8] enables the accurate approximation of solutions that involve jumps, kinks, singularities, and other locally nonsmooth features within elements. This is achieved by enriching the polynomial approximation space of the classical finite element method. The generalized finite element method (GEFM)/XFEM [9–11] has shown its potential in a variety of applications that involve nonsmooth solutions near interfaces. Among them are the simulation of cracks, shear bands, dislocations, solidification, and multifield problems. A large number of examples can be found in the real world where field quantities change rapidly over length scales that are small with respect to the observed domain. For the modeling of such phenomena, the resulting solutions typically involve discontinuities, singularities, high gradients, or other nonsmooth properties. For the approximation of nonsmooth solutions, one of the approaches is to enrich a polynomial approximation space such that the nonsmooth solutions can be modeled independently of the mesh. Methods that extend polynomial approximation spaces are called “enriched methods” in this work. The enrichment can be achieved by adding special shape functions (which are customized so as to capture jumps, singularities, etc.) to the polynomial approximation space. The total stiffness matrix consists of four parts: the FEM part, the enrich nodes part, and another two parts are the connecting parts between the real nodes and enrich nodes. Compared with FEM, the matrix size increases , while the FEM stiff matrix size is . But engineering calculations in XFEM always encountered the problem of solving ill-conditioned matrix. It made the XFEM widely.
In this paper, according to the characteristics of the stiffness matrix, we use an iterative method to avoid dealing with complicated stiff matrix. First we can use the result of linear equation from standard FEM. The result can be used to be the initial value for iteration. Then we decompose the unknown values to two parts; each part is governed by two linear equations, whose dimension is smaller than the whole stiff matrix. The new programm can be easily applied. Finally, numerical examples show that the proposed method is more efficient than common methods. Furthermore, the new method works well for some problems while the common method (such as LU) fails. The new iterative method can be applied to many other problems.
The outline of this paper is as follows. Section 2 introduces the XFEM. Subsequently, the main part of this paper is Section 3 that we describe a new efficient and feasible method. We also analyze the time complexity of the algorithm. In Section 4, we briefly mention details on numerical examples. A final conclusion is drawn in Section 5.
2. The Extended Finite Element Method
The basic mathematical foundation of the partition of unity finite element method (PUFEM) was discussed in [12, 13]. They illustrated that PUFEM can be used to employ the structure of the differential equation under consideration to construct effective and robust methods. The global solution of PUFEM has been the theoretical basis of the local partition of unity finite element method, to be called later the extended finite element method.
The first effort for developing the extended finite element method can be traced back to 1999, [7] presented a minimal enmeshing finite element method for crack growth. They added discontinuous enrichment functions to the finite element approximation to account for the presence of the crack. The method allowed the crack to be arbitrarily aligned within the mesh, though it required enmeshing for severely curved cracks. The method was improved and called extend Finite Element Method (XFEM) in [8, 9, 11, 14–16].
In mathematics, relevant theoretical parts of doing are almost same in [17–19]. Of course, XFEM is still evolving currently and keeping the new extension about (such as fluid-structure interaction): mathematics stiff matrix remained prone to the sick.
Consider a domain in which is discredited by a set of elements such that and a set of nodes
The total ordering of the element basis is given by and the total ordering of the nodes is given by . The standard finite element basis is defined as follows: where is the finite element basis or shape function of node ; is assumed to be of compact support and piecewise continuously differentiable. Typically, spans the space of piecewise continuous polynomials of a specific order. The XFEM aims to alleviate the burden associated with mesh generation for problems with voids and interfaces.. It does not require the finite element mesh to conform to internal boundaries. The essence of the XFEM lies in subdividing a model into two distinct parts: mesh generation for the domain (excluding internal boundaries) and enrichment of the finite-element approximation by additional functions. A partitioning of a typical domain into its enriched and unenriched subdomains is shown in Figure 2(a). We defined the subdomain where the enrichment is applied as This region’s feature is that the enrichment is dominant. The elements covering are the enriched elements, so where is the ordering of the subset of elements which are to be enriched. Any node for which the support of overlaps or is contained in is enriched node, and the set is defined by the ordering , . For an enrichment function , the enriched basis for the local partition of unity method is The approximation enriched with a local partition of unity the field is given by the sum of a standard finite element approximation and a linear combination of the enriched basis functions where and are the shape functions for the standard part and the partition of unity, respectively; different order interpolates may be used for the standard and enriched shape functions.
For example, we consider the linear elastostatics governing equations where is cauchy stress and is the body force per unit volume. The first function is called equilibrium equation. The second function is the equilibrium equation that describes the relation of strain and displacement. The third equation is for Hooke’s law. The parameter implies two parameters and shear elasticity, (Lamé constants are , and Poisson’s ratio is ). Consider The boundary is composed of the sets , , , and such that where is the unit outward normal to and and are prescribed displacements and traction, respectively. The weak formulation of (12) is well known and reads as follows: find such that Here and .
For the finite element discretization, let be a regular rectangular element of domain and define the mesh parameter .
Consider the Bubnov-Galerkin implementation for the XFEM in two-dimensional linear elasticity. In the XFEM, finite-dimensional subspaces and are used as the approximating trial and test spaces. The finite dimensional space is given by The weak form of the discrete problem can be stated as follows: find such that where and . In a Bubnov-Galkerkin procedure, the trial function and the test function are represented as linear combinations of the same shape functions. The trial and test functions are where denotes the set of all nodes in the mesh and denotes the set of nodes near the interface. are the finite-element shape functions, is the level set function, and is the enrichment function for material interface.
There are several kinds of enrichment functions, such as abs-enrichment function, step-enrichment function, and Ramp-enrichment function. The choice of function is based on the behavior of the solution near the interface.
The Ramp function is defined as This enrichment function yields only continuous solutions. The advantage is that it automatically satisfies the continuity condition and does not require the use of Lagrange multipliers.
The step-enrichment function is defined as This enrichment function can yield a continuous or discontinuous solution across the interface but requires Lagrange multipliers to apply the Dirichlet jump condition.
On substituting the trial and test functions form (17), and using the arbitrariness of nodal variations, the following discrete system of linear equations is obtained: where for a finite-element displacement degree of freedom, and for an enriched degree of freedom. In the above equations, is the constitutive matrix for an isotropic linear elastic material: and the matrix is given by
We can get the integral terms, and the matrices computed from the integral terms are block matrices of the form
The next step is solving the linear equations system.
3. The Program for Solving the Equation
We need to solve the following equations: Here and , the block matrix is same in FEM and are sparse matrix.
We must point out that the matrix has up to nonzero elements. We can use the modified program, in (22) can be converted to two sub-problems: to solve and solve by -order linear equations. Now, we list the detailed steps.
Algorithm 1. (i) We use the ordinary method to solve
the results for the FEM equation and set the result as , . The matrix is block tridiagonal symmetric, sparse matrix. There are a lot of methods to solve this equation efficiently.
(ii) Solve the small scale equation
The matrix is a sparse tridiagonal, where the dimension is much smaller than . This equation can be solved quicker than (i).
(iii) Then we use (25) to solve ;
the computation time complexity is the same with the first one.
(iv) If the norm or (giving the maximum iterative times), output the result , . Else, return to step (ii).
In this part, we can use the common FEM solver such as LU or other useful solvers. We do not change the stiffness matrix and , these two matrixes are symmetric, positive, definite and banded. We can compute easily for (23), (24), (25). The matrix conditions is not changed; [17]. First of all, the store can choose a one-dimensional. We can choose a one dimensional storage when we use triangular decomposition for the matrix .
Remark 2. If the matrix and are not tridiagonal, we often use LU or GMRES method [20] to solve (23), (24), and (25).
The iteration has the following advantages. (1)Other than the half bandwidth of zero elements do not have storage, nor involved in the calculation, the amount of storage crunch, the impact of rounding error is relatively small, even if it is ill-conditioned matrix, we often obtain a higher accuracy.(2)The operator only uses FEM linear equation solver; it is a simple calculation.(3)The original matrix deposit triangular matrix and diagonal matrix deposit area can overlap, thus avoiding to take up a lot of the middle unit of work, saving memory.(4)This solution can be extended to block direct solution.(5)With the iterative method, it is easy to build and test the general program and you can more accurately estimate the computation time; the user is relatively easy to grasp the algorithm process.
Also if we using in [21] to solve the equation, equations use LU solver which need Multiplication, here if the number of iterations is , we use up to multiplications.
4. Numerical Test
In this section, we report two experiments to verify the efficacy and accuracy of the new iterative method in XFEM. We use Matlab and modify software [22] to demonstrate Algorithm 1. For the convenience of presentation, we introduce the following notations: SFEM means the standard finite element method; XFEM means the extended finite element method; means the mesh size; GMRES means the generalized minimal residual iterative method; -error means ; means the -error when the mesh size equals ; iterations mean the iterative times in the algorithm; orders mean the numerical error order by XFEM:
Here . First, we list Example 1 for one-dimensional problems.
Example 1. The exact solution is
Here and the coefficient , if, , if. And is point force from the righ; . We use Figure 1 to describe the error change when we use different . We can easily see that, when the step is smaller than , the whole computing time has economized one order of magnitude in Example 1.
(a)
(b)
(a)
(b)
Example 2. For the two-dimensional solid test, we use the same enrichment function and shape function to compute (4). We define with the material interface
Let Young’s modulus and Poisson’s ratio in
be , and let those in
be , , but area force in domain .
We denote and , and the exact solution is
Plate with a circular hole. Figure 2(b) shows the consistent rectangular subdivision and labels the nodes which need to be extended. The CPU time we counted is the total time for the entire programm running. Table 1 shows the linear equation scale and the condition numbers of FEM and XFEM, the density of sparse matrix show the matrix is sparse. It can be easily find the condition number of XFEM is bigger than that of FEM. The CPU times are also compared with LU solvers and GMRES in Table 2. The cheapest is LU by FEM and the new iterative method is the second cheapest; however, the error for LU by FEM is the biggest in all methods. Table 3 shows an -error comparison between the new iterative method and other solvers. And the numerical -error is nearly . -error, order, and iterations are given in Table 4.
5. Conclusion
In this paper, we discuss how to use XFEM for the stiff matrix ill-conditioned and came to the following conclusions. Firstly, we found why some of problems cannot be computed by XFEM; the stiff matrix is not then symmetric and ill-conditioned. So, a lot of storage space is wasted. The time of computing XFEM equation is longer, which cannot be computed in certain subdivision. Secondly, we use XFEM for solving two-dimensional solid test; when solving the XFEM equations, the iterative method and the standard finite element equations solver are adopted, so that we change the problems into solving the smaller scale equations and standard finite element stiffness matrix need not be changed. Finally, numerical examples show that the proposed method is more efficient than the common method. We can use this programm in many other interface problems later. Besides, we also consider combining the stable GFEM with our iterative method in the future. And it is believed that the XFEM can be widely used in more problems in fluid dynamics and material science. The calculation module in some software can be improved for users.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author also thanks the anonymous reviewers for their comments in improving the paper. This research was supported by the special funds for NSFC (1127313, 11171269, 11201369, and 61163027) and the Ph.D. Programs Foundation of the Ministry of Education of China (Grant no. 20110201110027), and it was partially subsidized by the Fundamental Research Funds for the Central Universities (Grant nos. 08142003 and 08143007).