- Email: [email protected]
Implementation details
CHAPTER ELEVEN
Implementation details 11.1. Computer implementation of enriched methods It is well known that the standard finite element method based on Lagrange polynomials have been developed for more than 50 years and have a rich history of both theoretical and practical improvements. By contrast, the enriched or partition-of-unity based methods, particularly non-standard basis functions are much newer and although have been intensively studied in recent years they have not reached the same level of acceptance and usage in computational mechanics. However, we mention that XFEM in particular has been natively implemented in some commercial software such as Abaqus and ANSYS. Abaqus also has built-in modules for the discrete element method (DEM) and of smoothed-particle hydrodynamics (SPH). It also allows the implementation of user-defined elements, which has been used to develop plug-ins for various other methods, including isogeometric analysis (IGA) [17]. IGA has also been subject to intensive research and it has been incorporated in some commercially available software such as LS-DYNA [3,4], FEAP [27], and Altair Radioss [20]. Among open-source software, we mention OOFEM [23] for XFEM and IGA, and various research codes such as OpenXFEM++ [6] and ElemFreGalerkin [19] for XFEM and meshfree methods, as well as GeoPDEs [29], IGAFEM [18], PetIGA [10], igatools [21], and IGAPack [2] for isogeometric analysis implementations. While established commercial and community-based projects have the advantage of better reliability, robustness and user support, they are not always well suited for learning the underlying methods. This is particularly true for large open-source and commercial projects which already contain many thousands of lines of code since one has to understand the already existing structure. Even for modular or object-oriented software that allows the development of customized user elements, extending the existing code to new capabilities is not always straightforward. From this point of view, shorter implementations written in a high-level language such as MATLAB® can be better suited for the understanding the basics of the implemented method and for extending it to different applications. Regardless of the details of the method and application considered, a finite element program typically consists of three stages: pre-processing, Extended Finite Element and Meshfree Methods Copyright © 2020 Elsevier Inc. https://doi.org/10.1016/B978-0-12-814106-9.00017-6 All rights reserved.
581
582
Extended Finite Element and Meshfree Methods
processing, and post-processing. In the pre-processing stage, the problem data such domain geometry, material properties, boundary conditions and other model parameters are initialized. In a mesh-based discretization, the mesh that describes the domain geometry is generated along with the node-element connectivities. For mesh-free methods, only a set of nodes (particles) is generated inside the domain, however some mesh-free methods use a background mesh for integration. In the case of strong or weak discontinuities which are tracked by enriched methods, typically a level set representation of the discontinuity must be generated as detailed in Section 1.4. In the processing stage of the program, the stiffness matrices and load vectors are assembled, the boundary conditions are imposed and a linear system is solved by either direct or iterative methods. In the case of nonlinear finite elements, typically a Newton-Raphson iteration is employed to obtain the solution. Finally, in the post-processing stage the quantities of interest such as displacements, strains, stresses are computed and plotted based on the solution vector. In the case of crack-propagation studies, the stress intensity factors are usually computed. Other post-processing tasks may involve the computation of error estimators or indicators (as detailed in Section 6.5 for hierarchical splines) and marking of elements or domain regions for refinement. In this section, we will give a general overview of the three stages mentioned above, highlighting some of differences and commonalities between the different methods considered. Later in the chapter, more in-depth details for XFEM, enriched meshfree and XIGA with numerical examples will be given. We will focus on the implementation details for the MATLAB/Octave language, although the same ideas can be implemented in other programming languages as well.
11.1.1 Pre-processing The most important part of this stage of a finite element program (and often the most time-consuming) is the mesh generation. Fortunately, mesh generators have improved quite significantly over the years and for most standard geometries the meshing tasks can be automated to a great extent. Still, manual intervention can be required which makes the meshing process very time-intensive in particular for problems involving moving boundaries or design optimization and has motivated the development of alternative discretizations as detailed in this book. Among the meshing methods, we differentiate between simplex meshes (which use triangles in 2D and tetrahedra in 3D) and quadrilateral/hexahe-
Implementation details
583
dral meshes. The former are easier to generate for arbitrary geometries and fairly robust implementations exist even in open-source programs such as Gmsh [14] or Tetgen [24]. In MATLAB, simplex meshes in 2D and 3D can be generated with the generateMesh command, which is part of the PDE toolbox. A drawback of triangular and tetrahedral meshes is that the linear elements tend to be inaccurate (too stiff, particularly in bending-dominated problems). Therefore, it is recommended to use quadratic or higher order elements, although this results in a significant increase in the number of degrees of freedom. Alternatively, hexahedral meshes can be used although they are typically more difficult to generate and the results can be quite sensitive to element distortion. In the case of meshfree methods, the nodal points can be generated from some structured grid, using some randomized algorithm such as latinhypercube sampling (LHS), or even by considering the element vertices resulting from a mesh-based approach as in [28]. As stated before, the benefits of meshfree approaches include more flexibility in generating the nodal points and insensitivity to mesh distortion. However, there are some computational costs associated with calculating the neighbor lists and numerical integration which may require a background mesh and an increased number of quadrature points. When considering the standard (tensor product) isogeometric analysis, the mesh requirements are in some ways similar to those for structured of quadrilateral and hexahedral meshes. However, the knot-vector and control-point structures which are typically used in IGA to define the domain geometry are different from those used in finite elements. A bridge that can allow easier integration of IGA concepts into existing finite element codes has been realized through the concept of Bézier extraction [7] as described in Section 6.1.2. This concept is also particularly useful in the case of hierarchical or non-tensor product discretizations, such as the PHT splines. We mention that there has also been interest the integration of unstructured mesh discretizations based on Bézier triangles [13], tetrahedra [31] and quadrilaterals [5,16]. For isogeometric analysis, a useful tool for generating parametrized geometries is the Octave NURBS toolbox [11,25], which is also compatible with MATLAB. An example for generating the plate with a hole geometry is given in Listing 11.1. In this file, we generate the model for a given radius of the hole R and a plate size L by creating the boundary segments, and parametrizing the interior with the Coon patch method [9] via the nrbcoons command. Degree elevation (p-refinement) to quadratic
584
Extended Finite Element and Meshfree Methods
NURBS and knot-insertion (h-refinement) are handled by the nrbdegelev and nrbkntins commands respectively. We refer to the documentation of the NURBS toolbox for more details on these commands. The plots showing the control points and the elements for a mesh with 10 × 10 knot-spans are given in Fig. 11.1. % s c r i p t p l o t P l a t e W i t h H o l e .m %g e n e r a t e s the geometry f o r a q u a r t e r p l a t e with a hole % u s e s t h e NURBS t o o l b o x ( h t t p s : / / o c t a v e . s o u r c e f o r g e . i o / n u r b s / i n d e x . h t m l ) addpath ( ' . / nurbs −1.3.13/ i n s t ' ) R = 1; %r a d i u s of the hole L = 4 ; % s i z e o f t h e modeled p l a t e % c r e a t e t h e c u r v e d b o u n d a r y and t h e o t h e r f o u r l i n e a r b o u n d a r y s e g m e n t s a r c = n r b c i r c (R, [ 0 , 0 ] , p i / 2 , p i ) ; s i d e B o t t o m = n r b l i n e ([ −R, 0 ] , [ − L , 0 ] ) ; s i d e R i g h t = n r b l i n e ( [ 0 , R] , [ 0 , L ] ) ; pointUpRight = [ 0 ; L ; 0 ] ; p o i n t U p L e f t = [ −L ; L ; 0 ] ; p o i n t D o w n L e f t = [ −L ; 0 ; 0 ] ; knotUpLeft = [ 0 , 0 , 1/2 , 1 , 1 ] ; s i d e U p L e f t = nrbmak ( [ p o i n t U p R i g h t , p o i n t U p L e f t , p o i n t D o w n L e f t ] , k n o t U p L e f t ) ; %d e g r e e e l e v a t e t h e l i n e a r s i d e and c r e a t e a 2− e l e m e n t v o l u m e t r i c model sideBottom = nrbdegelev ( sideBottom , 1) ; plateModel = nrbcoons ( arc , sideUpLeft , sideRight , sideBottom ) ; %p e r f o r m k n o t i n s e r t i o n s t o r e f i n e t h e model newKnotsU = [ 0 . 1 : 0 . 1 : 0 . 4 , 0 . 6 : 0 . 1 : 0 . 9 ] ; newKnotsV = [ 0 . 1 : 0 . 1 : 0 . 9 ] ; p l a t e M o d e l = n r b k n t i n s ( p l a t e M o d e l , { newKnotsU , newKnotsV } ) ; % p l o t t h e model and t h e c o n t r o l p o i n t s figure subplot (1 ,2 ,1) , n r b c t r l p l o t ( plateModel ) view ( 2 ) subplot (1 ,2 ,2) , nrbkntplot ( plateModel ) view ( 2 )
Listing 11.1: Matlab code for generating the plate with a hole model. In the framework of enriched methods, additional structures are needed to represent the sets of enriched nodes and the type of enrichment. A unified way to implement this concept is to define an array corresponding to each node that determines its set of enrichment functions. For the functions that are not enriched, this array will consist of just the unity function {ξ1 (x) = 1}. For other functions, for example the Heavisideenriched basis, the array will also include corresponding enrichments, e.g. {ξ1 (x) = 1, ξ2 (x) = H (x)}. Then the derivatives of the basis functions of the form φi ξji can be computed by the product rule from the derivatives of φi , i = 1, . . . , N, and ξji , j = 1, . . . , ni where N is the number of (unenriched)
Implementation details
585
Figure 11.1 Plate with a hole model generated by Listing 11.1, showing the control points on the left and the refined mesh on the right.
global shape functions and ni is the number of enrichments (including the unity) for the ith basis function. Finally, in the case level-sets are used to track the strong or weak discontinuities, a discretization for the level-sets is needed, usually in terms of the global basis functions. We refer to the included codes and Sections 1.4 and 8.3 for additional implementation details.
11.1.2 Processing The processing stage of a finite element program consists of constructing and solving a linear system which is obtained by assembling the stiffness matrices and load vector and by imposing the boundary conditions. The assembly stage is typically conducted as an outer loop over the elements (or integration cells in the case of meshfree methods) and an inner loop over the quadrature Gauss points. In the case of a MATLAB implementation, it is preferable that the values of all the shape functions with support at a particular Gauss point are output as a vector by the basis function subroutine. Then the process of assembling the local stiffness matrices can be vectorized, resulting in a significant performance gain. Another MATLAB-specific speedup is the use of the so-called “triplet arrays” for the assembly of the sparse-format stiffness matrices. The idea is to replace matrix update statements of the form K(rowIndex,colIndex)=K(rowIndex,colIndex)+localK with statements like the ones in Listing 11.2. The reason is that the sparse matrices are stored in the compressed sparse column (CSC) format [15], for which adding additional non-zero
586
Extended Finite Element and Meshfree Methods
entries requires expensive resorting operations. Therefore, it is much more efficient to store the non-zero entries in triplet arrays and build the sparse matrix at the end. A drawback of this procedure is that the triplet arrays consume more memory than the sparse matrix, particularly when there are many repeated row and column indices whose corresponding value needs to be summed. For large matrices (with more than 100,000 degrees of freedom) or when there is a lot of overlap between the basis functions, it may be required to break up the assembly into subdomains and sum up the assembled stiffness matrices. This procedure can also be helpful in the parallelization of the assembly subroutine. iCounter = 0 % A l l o c a t e s p a c e f o r t h e a r r a y s rowSet , c o l S e t , v a l S e t %. . . f o r i n d e x E l e m = 1 : numElements % s e t c o l I n d e x and r o w I n d e x from t h e c o n n e c t i v i t y a r r a y s and %a s s e m b l e t h e l o c a l s t i f f n e s s m a t r i x by l o o p i n g o v e r G a u s s p o i n t s %. . . numEnt = l e n g t h ( r o w I n d e x ) ; r o w S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e p m a t ( rowIndex , 1 , numEnt ) ; c o l I n d e x S e t = r e p m a t ( c o l I n d e x , numEnt , 1 ) ; c o l S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e s h a p e ( c o l I n d e x S e t , 1 , numEnt ^ 2 ) ; v a l S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e s h a p e ( l o c a l K , 1 , nument ^ 2 ) ; i C o u n t e r = i C o u n t e r + numEnt ^ 2 ; end K = s p a r s e ( rowSet , c o l S e t , v a l S e t , dimK , dimK ) ;
Listing 11.2: Snippet of Matlab code for triplet matrix assembly. For the imposition of the boundary conditions, we distinguish between Dirichlet (e.g. prescribed displacement) and Neumann (e.g. prescribed traction) boundary conditions. The Dirichlet boundary conditions are imposed by changing the stiffness matrix and right-hand side vector, while for Neumann boundary conditions only the right-hand side vector is changed. In the case of finite element method with Lagrange polynomials, a common procedure is to the zero-out the rows of the stiffness matrix corresponding to the degrees of freedom on the Dirichlet boundaries, set the diagonal entries in these rows to one and set the corresponding entries in the righthand side vector to the prescribed value. This has the drawback that the resulting stiffness matrix is no longer symmetric, which makes the process of solving the linear system less efficient. Symmetry can be restored by condensing out the prescribed degrees of freedom. To give a simple example, suppose we have a matrix with 4 degrees of freedom, of which the second one is given by x2 = u¯ , with u¯ known. The initial linear symmetric system
587
Implementation details
can be written as: ⎡
k11 ⎢k ⎢ 12 ⎢ ⎣k13 k14
k12 k22 k23 k24
k13 k23 k33 k34
⎤⎡ ⎤
⎡ ⎤
k14 x1 f1 ⎢ ⎥ ⎢ ⎥ k24 ⎥ ⎥ ⎢x2 ⎥ ⎢f2 ⎥ ⎥⎢ ⎥ = ⎢ ⎥ k34 ⎦ ⎣x3 ⎦ ⎣f3 ⎦ k44 x4 f4
(11.1)
After setting x2 = u¯ , we now have: ⎡
⎤⎡ ⎤
⎡ ⎤
k11 k12 k13 k14 x1 f1 ⎢0 ⎥ ⎢x ⎥ ⎢ u¯ ⎥ 1 0 0 ⎢ ⎥ ⎢ 2⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ = ⎢ ⎥ ⎣k13 k23 k33 k34 ⎦ ⎣x3 ⎦ ⎣f3 ⎦ k14 k24 k34 k44 x4 f4
(11.2)
To restore the symmetry, we substitute x2 = u¯ in the equations corresponding to the rows 1, 3 and 4 and move the known quantities to the right hand side: ⎡
k11 ⎢0 ⎢ ⎢ ⎣k13 k14
⎤⎡ ⎤
⎡
⎤
0 k13 k14 x1 f1 − k12 u¯ ⎢ ⎥ ⎢ ⎥ ⎥ 1 0 0 ⎥ ⎢x2 ⎥ ⎢ u¯ ⎥ ⎥⎢ ⎥ = ⎢ ⎥ 0 k33 k34 ⎦ ⎣x3 ⎦ ⎣f3 − k23 u¯ ⎦ f4 − k24 u¯ 0 k34 k44 x4
(11.3)
This procedure can naturally be applied for the case of multiple prescribed degrees of freedom. A MATLAB implementation of this algorithm is given in Listing 11.3. In this code, additionally the diagonal entries on the stiffness matrix corresponding to the Dirichlet degrees of freedom and the entries in the load vector are multiplied by the mean value of the diagonal of the original matrix to improve conditioning. We note that in the case of IGA or mesh methods, the shape functions on the boundary do not, in general, satisfy the Kronecker-Delta property. However, the B-Spline and NURBS functions, as well as some meshfree basis such as the local maximum entropy (LME) basis functions on convex domains or generalized meshfree (GMF) approximations [30] do satisfy a weak Kronecker-delta property. In this case, the shape functions corresponding to the interior nodes (or control points) do not have support on the boundary, however the value of the basis functions on the boundary may not reach 1 at the corresponding nodes as in the strong Kroneckerdelta case. This requires an interpolation or projection of the boundary date on the approximation space spanned by the boundary basis functions. In the case where interior shape functions also have support on the bound-
588
Extended Finite Element and Meshfree Methods
f u n c t i o n [ kk , f f ]= a p p l y D i r i c h l e t S y m ( kk , f f , b c d o f , b c v a l ) % A p p l y t h e D i r i c h l e t c o n s t r a i n t s t o m a t r i x e q u a t i o n [ kk ] { x }={ f f } % w h i l e o u t p u t t i n g a s y m m e t r i c m a t r i x kk % Input : % kk − s y s t e m m a t r i x b e f o r e a p p l y i n g c o n s t r a i n t s % f f − system vector before applying c o n s t r a i n t s % bcdof − a vector containging c o n s t r a i n e d d . o . f % bcval − a vector containing constrained value % Output : % kk − t h e new s y s t e m m a t r i x a f t e r a p p l y i n g t h e c o n s t r a i n t s % f f − t h e new l o a d v e c t o r a f t e r a p p l y i n g t h e c o n s t r a i n t s bcwt=mean ( d i a g ( kk ) ) ; % a m e a s u r e o f t h e a v e r a g e s i z e o f an e l e m e n t i n K f f = f f −kk ( : , b c d o f ) ∗ b c v a l ' ; % m o d i f y t h e force vector f f ( b c d o f ) = bcwt ∗ b c v a l ; kk ( b c d o f , : ) =0; % z e r o o u t t h e r o w s and c o l u m n s o f t h e K m a t r i x kk ( : , b c d o f ) = 0 ; kk ( b c d o f , b c d o f ) =bcwt ∗ s p e y e ( l e n g t h ( b c d o f ) ) ; % p u t o n e s ∗ bcwt on t h e d i a g o n a l
Listing 11.3: Matlab code for imposing Dirichlet boundary conditions. ary (such as in the case of MLS or RKP meshfree methods), then other techniques such as Lagrange multipliers are required. After the linear system is changed to take into account the boundary conditions, the solution vector can be obtained by either a direct or iterative solver. Direct solvers (such as the MATLAB backslash “\” operator) tend to be more robust for ill-conditioned linear systems, but the memory requirements can become very high, particularly for 3D problems with a lot of degrees of freedom. The reason is due to the “fill-in” which occurs when the zero entries within the bandwidth of a matrix become non-zero when the matrix is factorized. Iterative (Krylov subspace) solvers rely on matrix vector multiplications and can be used to solve very large linear systems efficiently. Examples of iterative solvers are pre-conditioned conjugate gradient (pcg command in MATLAB) which is efficient for symmetric positive definite matrices or the generalized minimal residual method (GMRES, gmres command in MATLAB) for general matrices. These solvers work better when a pre-conditioner is applied to improve the convergence. Standard preconditioners are incomplete Cholesky factorizations (ichol command in MATLAB) for symmetric positive definite matrices and the incomplete LU factorization (ilu command in MATLAB) for general sparse matrices. The idea of incomplete factorization is that only the non-zero entries in the original matrix can be non-zero in the factorization, the other entries are set to zero. Other preconditioners, such as Jacobi preconditioners which
Implementation details
589
perform a diagonal scaling of the matrices, or multi-grid preconditioners can be used as well. The solution of linear systems also appears in nonlinear analysis, where usually Newton-Raphson iterations are performed to solve a linearized model. In this case, on the left-hand side we have the “tangent-stiffness” matrix, as described e.g. in Section 3.6.1.3 for hydro-mechanical coupled problems, and on the left-hand side we have a residual vector. Solving linear systems is also required for time-dependent problems with implicit timeintegration (where a system involving the stiffness matrix is solved) or with explicit time-integration where the consistent mass matrix is used. Therefore, for large simulations, the use of efficient solvers or of accurate lumped mass methods is of great importance in keeping the computational times within acceptable limits.
11.1.3 Post-processing In the post-processing stage, various quantities of interest can be computed and displayed based on the solution vector. A correct and detailed visualization of the output date is very important, particularly for identifying any errors in the model such as incorrect imposition of boundary conditions. In the case of XFEM or basis functions which do not satisfy the (strong) Kronecker-delta property, the displacements cannot be simply inferred from the values of the solution vector. Instead, the basis functions (including the enrichments) need to be evaluated at the chosen points (e.g. element vertices) and the values multiplied by the corresponding coefficients of the solution. The strain and stress values are calculated in the same way, by considering the derivatives of the basis functions and of the enrichments. For meshless methods, typically a visualization mesh needs to be constructed for example from a Delaunay triangulation of the domain. Moreover, to precisely show the discontinuities in a fracture simulation, the visualization mesh should be conforming to the crack path. For many examples, we use output files in the VTK (Visualization Toolkit) format. These files can be opened with open-source viewers such as Paraview, and are particularly convenient for 3D visualization as one can easily manipulate the graphs and interactively plot the deformations, stresses and other quantities of interest on structured or unstructured meshes. Highlevel output subroutines such as vtkwrite1 are available to make the output of scalar or vector quantities relatively straightforward. 1 https://github.com/joe-of-all-trades/vtkwrite.
590
Extended Finite Element and Meshfree Methods
For fracture propagation problems, a crack increment is computed according to some failure as well as propagation magnitude and angle criteria. In XFEM (and other enriched methods), stress intensity factors are often used for the computation of these criteria. They can be calculated in 2D by the interaction integral methods [26], which are derived from the J-integral [22], as discussed in Sections 5.9.1.4 and 7.2.2.5. This method considers a contour around the crack tip and is converted to a domain integral for improved integration accuracy. Since the resulting code is quite lengthy, we refer to the included examples and the documentation for more details. Finally, even though the idea of enriched methods is that the mesh is independent of the crack, often times the results can be improved by considering a finer mesh around the crack tips. In the framework of hierarchical splines, this can be conveniently implemented with a recovery-type error estimator, as detailed in [2]. This type of estimator is quite effective at refining in regions with high stress concentrations. Another aspect to be considered is the selection of a marking scheme and the amount of elements being marked at each refinement step. When the mesh is coarse, the error can be quite large even relatively far from the regions with non-smooth solution such as the crack tip. As the mesh is refined and the singular regions are better resolved, the areas with large error can shrink considerably. Therefore, it is more economical in terms of degrees of freedom to refine as few elements as possible at each step, however this leads to more refinement steps which increases the computational cost. A strategy that works well and has been implemented in [2] is the Dörfler marking scheme [12], which refines elements are in the 100θ percentile of the largest error. For problems with singularity, we typically choose θ = 0.5, while for problems with smoother solutions (but which can still benefit from adaptivity) θ = 0.75 is more efficient.
11.2. Numerical examples In this section we provide a few numerical examples of benchmark problems in XFEM, including that of a coupled hydro-mechanical model. Further examples and documentation, including 3D, plate and shell models can be found in the open-source repository IGA(X)FEM at https:// sourceforge.net/projects/cmcodes.
591
Implementation details
Figure 11.2 Schematic illustration of mechanical loading.
Table 11.1 Material parameters for model under pure mechanical loading. Young’s modulus E = 25.85 GPa Poisson’s ratio ν = 0.18 uniform stress field σ1 = 1.0e7 Pa crack length 2a = 2.15 m
11.2.1 Crack propagation angle This section presents a square domain with an inclined center crack under pure mechanical loading (Fig. 11.2). Hereafter, the initial propagation angles of the induced crack with different pre-defined angle β and loading factor α are studied. The material parameters are list in Table 11.1. Fig. 11.3 shows the stress profiles of domains with pre-defined angle β ranging from 10 to 80 degree and loading factor α is set to 0. The comparison between numerical calculation and theoretical model [1] in terms of the propagation angle of the induced crack with different loading ratio α is shown in Fig. 11.4A. Furthermore, the convergence study (case β = 10) of stiffness intensity factor K1 of mode I respect to total number of elements is presented in Fig. 11.4B.
11.2.2 Hydro-mechanical model with center cracks This section presents the domain with inner fractures under fluid flux loading. The geometry and boundary conditions are shown in Fig. 11.5. The material parameters are listed in Table 11.2. Domain with single and five
592
Extended Finite Element and Meshfree Methods
Figure 11.3 Stress profile in x direction with different pre-defined angle β s.
randomly distributed inner cracks are studied (Fig. 11.6). Fig. 11.6 clearly indicates that the existing crack has great impact on deformation, pressure profiles of the whole domain.
11.2.3 Hydro-mechanical model with edge crack This section simulates the rupture of a plate under plane-strain conditions. As shown in Fig. 11.7A, a pre-existing fracture is located at the symmetry
Implementation details
593
Figure 11.4 (A) Comparison results of the propagation angle of the induced crack between numerical calculation and theoretical model with different loading ratio α , solid lines represent the result from theoretical model. (B) Convergence study of stiffness intensity factor K1 of mode I respect to total element number.
Figure 11.5 Schematic illustration of geometry and boundary conditions. p represents pressure, q represents fluid flux. Table 11.2 Material parameters for hydro-mechanical model with inner cracks. Young’s modulus E = 25.85 GPa Poisson’s ratio ν = 0.18 Initial porosity n = 0.2 Intrinsic permeability k = 2.78e-21 m2 Biot’s coefficient α=1 Solid bulk modulus Ks = 13.46 GPa Water bulk modulus Kf = 0.2 GPa ρ = 1000 kg/m3 Water density Water dynamic viscosity μw = 5e-4 N/m2 s Irreducible saturation Sirr = 0.0 Reference pressure Rref = 18.6 MPa Porous index l = 0.4396
594
Extended Finite Element and Meshfree Methods
Figure 11.6 Displacement, pressure and stress profile for domain with single and five cracks.
axis of the plate. A fixed vertical velocity v = 2.35 × 10−2 µm/s is prescribed in an opposite direction at the bottom and at the top of the plate (tensile loading). All boundaries of the plate are assumed to be impervious. The induced pressure profile is presented in Fig. 11.7B.
11.2.3.1 Hydro-mechanical model with edge crack under fluid flux loading This subsection follows an experimental setup [8] (Fig. 11.8). The setup of this experiment can be summarized as follows: a square sample (Fig. 11.8A) with two ports, high pore pressure is introduced to the lower port region by injecting pressurized fluid and maintain the pressure, hereafter a wedge
Implementation details
595
Figure 11.7 (A) Schematic illustration of geometry and boundary conditions, (B) pressure profile.
Figure 11.8 (A) Schematic illustration of geometry and boundary conditions, (B) experimental results.
is driven into the pre-shaped concave edge to induce fracture and fracture propagation. The crack propagation path is then identified. Further more, the impact of the permeability of the solid on the crack propagation path is also studied (Fig. 11.9). The fracture development path, displacement of the domain and stress profiles are shown in Fig. 11.9 for porous medias (under identical loading) with different intrinsic permeabilities (k1, k2, k3). Fig. 11.9D-F clearly demonstrates that the fracture paths are attracted to the fluid pressurized port. Furthermore, we found that under identical injection rate, porous media with lower permeability (e.g. k1) will induce
596
Extended Finite Element and Meshfree Methods
Figure 11.9 Displacement (m) plot of porous media with the intrinsic permeabilities (A) k1 = 2.78e-18 m2 ; (B) k2 = 1.52e-17 m2 ; (C) k3 = 2.78e-17 m2 and (D)-(F) are the stress (Pa) plots in cylinder coordinate accordingly.
Implementation details
597
higher stress/pressure state around the injection port comparing to porous media with higher permeability (e.g. k3), which, in return influences the propagation direction of the fracture.
References [1] T. Anderson, Fracture Mechanics: Fundamentals and Applications, fourth edition, CRC Press, 2017. [2] C. Anitescu, M.N. Hossain, T. Rabczuk, Recovery-based error estimation and adaptivity using high-order splines over hierarchical t-meshes, Computer Methods in Applied Mechanics and Engineering 328 (2018) 638–662. [3] D. Benson, Y. Bazilevs, M. Hsu, T. Hughes, Isogeometric shell analysis: the Reissner– Mindlin shell, Computer Methods in Applied Mechanics and Engineering 199 (5) (2010) 276–289. [4] D. Benson, Y. Bazilevs, M.-C. Hsu, T. Hughes, A large deformation, rotation-free, isogeometric shell, Computer Methods in Applied Mechanics and Engineering 200 (13) (2011) 1367–1378. [5] M. Bercovier, T. Matskewich, Smooth Bézier Surfaces Over Unstructured Quadrilateral Meshes, Springer International Publishing, 2017. [6] S. Bordas, V. Nguyen, C. Dunant, H. Nguyen-Dang, A. Guidoum, An extended finite element library, International Journal for Numerical and Analytical Methods in Engineering 71 (6) (2007) 703–732. [7] M.J. Borden, M.A. Scott, J.A. Evans, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of NURBS, International Journal for Numerical Methods in Engineering 87 (1–5) (2011) 15–47. [8] M. Bruno, F. Nakagawa, Pore pressure influence on tensile fracture propagation in sedimentary rock, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28 (4) (1991) 261–273. [9] S.A. Coons, Surfaces for Computer-Aided Design of Space Forms, Technical report, Cambridge, MA, USA, 1967. [10] L. Dalcin, N. Collier, P. Vignal, A. Côrtes, V. Calo, PetIGA: a framework for highperformance isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 308 (2016) 151–181. [11] C. de Falco, A. Reali, R. Vázquez, GeoPDEs: research tool for isogeometric analysis of PDEs, Advances in Engineering Software 42 (12) (2011) 1020–1034. [12] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM Journal on Numerical Analysis 33 (3) (1996) 1106–1124. [13] L. Engvall, J.A. Evans, Isogeometric triangular Bernstein–Bézier discretizations: automatic mesh generation and geometrically exact finite element analysis, Computer Methods in Applied Mechanics and Engineering 304 (2016) 378–407. [14] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331. [15] J. Gilbert, C. Moler, R. Schreiber, Sparse matrices in MATLAB: design and Implementation, SIAM Journal on Matrix Analysis and Applications 13 (1) (1992) 333–356. [16] M. Kapl, G. Sangalli, T. Takacs, An isogeometric C1 subspace on unstructured multipatch planar domains, Computer Aided Geometric Design 69 (2019) 55–75.
598
Extended Finite Element and Meshfree Methods
[17] Y. Lai, Y.J. Zhang, L. Liu, X. Wei, E. Fang, J. Lua, Integrating CAD with Abaqus: a practical isogeometric analysis software platform for industrial applications, in: HighOrder Finite Element and Isogeometric Methods 2016, Computers & Mathematics with Applications 74 (7) (2017) 1648–1660. [18] V.P. Nguyen, C. Anitescu, S.P. Bordas, T. Rabczuk, Isogeometric analysis: an overview and computer implementation aspects, Mathematics and Computers in Simulation 117 (2015) 89–116. [19] V.P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshfree methods: a review and key computer implementation aspects, Mathematics and Computers in Simulation 79 (3) (2008) 763–813. [20] M. Occelli, T. Elguedj, S. Bouabdallah, L. Morançay, LR B-Splines implementation in the Altair RadiossTM solver for explicit dynamics IsoGeometric Analysis, Advances in Engineering Software 131 (2019) 166–185. [21] M. Pauletti, M. Martinelli, N. Cavallini, P. Antolin, Igatools: an isogeometric analysis library, SIAM Journal on Scientific Computing 37 (4) (2015) C465–C496. [22] J.R. Rice, A path independent integral and the approximate analysis of stress concentration by notches and cracks, Journal of Applied Mechanics 35 (1968) 379–386. [23] D. Rypl, B. Patzák, From the finite element analysis to the isogeometric analysis in an object oriented computing environment, in: CIVIL-COMP, Advances in Engineering Software 44 (1) (2012) 116–125. [24] H.Si. TetGen, A Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathematical Software 41 (2) (2015) 1–36. [25] M. Spink, D. Claxton, C. de Falco, R. Vázquez, The NURBS toolbox, http://octave. sourceforge.net/nurbs/index.html. [26] M. Stern, E. Becker, R. Dunham, A contour integral computation of mixed-mode stress intensity factors, International Journal of Fracture 12 (3) (1976) 359–368. [27] R.L. Taylor, FEAP - Finite Element Analysis Program, 2014. [28] N. Valizadeh, Y. Bazilevs, J. Chen, T. Rabczuk, A coupled IGA–meshfree discretization of arbitrary order of accuracy and without global geometry parameterization, Computer Methods in Applied Mechanics and Engineering 293 (2015) 20–37. [29] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Computer Mathematics and Its Applications 72 (3) (2016) 523–554. [30] C. Wu, C. Park, J. Chen, A generalized approximation for the meshfree analysis of solids, International Journal for Numerical Methods in Engineering 85 (6) (2011) 693–722. [31] S. Xia, X. Qian, Isogeometric analysis with Bézier tetrahedra, in: Special Issue on Isogeometric Analysis: Progress and Challenges, Computer Methods in Applied Mechanics and Engineering 316 (2017) 782–816.
Implementation details 11.1. Computer implementation of enriched methods It is well known that the standard finite element method based on Lagrange polynomials have been developed for more than 50 years and have a rich history of both theoretical and practical improvements. By contrast, the enriched or partition-of-unity based methods, particularly non-standard basis functions are much newer and although have been intensively studied in recent years they have not reached the same level of acceptance and usage in computational mechanics. However, we mention that XFEM in particular has been natively implemented in some commercial software such as Abaqus and ANSYS. Abaqus also has built-in modules for the discrete element method (DEM) and of smoothed-particle hydrodynamics (SPH). It also allows the implementation of user-defined elements, which has been used to develop plug-ins for various other methods, including isogeometric analysis (IGA) [17]. IGA has also been subject to intensive research and it has been incorporated in some commercially available software such as LS-DYNA [3,4], FEAP [27], and Altair Radioss [20]. Among open-source software, we mention OOFEM [23] for XFEM and IGA, and various research codes such as OpenXFEM++ [6] and ElemFreGalerkin [19] for XFEM and meshfree methods, as well as GeoPDEs [29], IGAFEM [18], PetIGA [10], igatools [21], and IGAPack [2] for isogeometric analysis implementations. While established commercial and community-based projects have the advantage of better reliability, robustness and user support, they are not always well suited for learning the underlying methods. This is particularly true for large open-source and commercial projects which already contain many thousands of lines of code since one has to understand the already existing structure. Even for modular or object-oriented software that allows the development of customized user elements, extending the existing code to new capabilities is not always straightforward. From this point of view, shorter implementations written in a high-level language such as MATLAB® can be better suited for the understanding the basics of the implemented method and for extending it to different applications. Regardless of the details of the method and application considered, a finite element program typically consists of three stages: pre-processing, Extended Finite Element and Meshfree Methods Copyright © 2020 Elsevier Inc. https://doi.org/10.1016/B978-0-12-814106-9.00017-6 All rights reserved.
581
582
Extended Finite Element and Meshfree Methods
processing, and post-processing. In the pre-processing stage, the problem data such domain geometry, material properties, boundary conditions and other model parameters are initialized. In a mesh-based discretization, the mesh that describes the domain geometry is generated along with the node-element connectivities. For mesh-free methods, only a set of nodes (particles) is generated inside the domain, however some mesh-free methods use a background mesh for integration. In the case of strong or weak discontinuities which are tracked by enriched methods, typically a level set representation of the discontinuity must be generated as detailed in Section 1.4. In the processing stage of the program, the stiffness matrices and load vectors are assembled, the boundary conditions are imposed and a linear system is solved by either direct or iterative methods. In the case of nonlinear finite elements, typically a Newton-Raphson iteration is employed to obtain the solution. Finally, in the post-processing stage the quantities of interest such as displacements, strains, stresses are computed and plotted based on the solution vector. In the case of crack-propagation studies, the stress intensity factors are usually computed. Other post-processing tasks may involve the computation of error estimators or indicators (as detailed in Section 6.5 for hierarchical splines) and marking of elements or domain regions for refinement. In this section, we will give a general overview of the three stages mentioned above, highlighting some of differences and commonalities between the different methods considered. Later in the chapter, more in-depth details for XFEM, enriched meshfree and XIGA with numerical examples will be given. We will focus on the implementation details for the MATLAB/Octave language, although the same ideas can be implemented in other programming languages as well.
11.1.1 Pre-processing The most important part of this stage of a finite element program (and often the most time-consuming) is the mesh generation. Fortunately, mesh generators have improved quite significantly over the years and for most standard geometries the meshing tasks can be automated to a great extent. Still, manual intervention can be required which makes the meshing process very time-intensive in particular for problems involving moving boundaries or design optimization and has motivated the development of alternative discretizations as detailed in this book. Among the meshing methods, we differentiate between simplex meshes (which use triangles in 2D and tetrahedra in 3D) and quadrilateral/hexahe-
Implementation details
583
dral meshes. The former are easier to generate for arbitrary geometries and fairly robust implementations exist even in open-source programs such as Gmsh [14] or Tetgen [24]. In MATLAB, simplex meshes in 2D and 3D can be generated with the generateMesh command, which is part of the PDE toolbox. A drawback of triangular and tetrahedral meshes is that the linear elements tend to be inaccurate (too stiff, particularly in bending-dominated problems). Therefore, it is recommended to use quadratic or higher order elements, although this results in a significant increase in the number of degrees of freedom. Alternatively, hexahedral meshes can be used although they are typically more difficult to generate and the results can be quite sensitive to element distortion. In the case of meshfree methods, the nodal points can be generated from some structured grid, using some randomized algorithm such as latinhypercube sampling (LHS), or even by considering the element vertices resulting from a mesh-based approach as in [28]. As stated before, the benefits of meshfree approaches include more flexibility in generating the nodal points and insensitivity to mesh distortion. However, there are some computational costs associated with calculating the neighbor lists and numerical integration which may require a background mesh and an increased number of quadrature points. When considering the standard (tensor product) isogeometric analysis, the mesh requirements are in some ways similar to those for structured of quadrilateral and hexahedral meshes. However, the knot-vector and control-point structures which are typically used in IGA to define the domain geometry are different from those used in finite elements. A bridge that can allow easier integration of IGA concepts into existing finite element codes has been realized through the concept of Bézier extraction [7] as described in Section 6.1.2. This concept is also particularly useful in the case of hierarchical or non-tensor product discretizations, such as the PHT splines. We mention that there has also been interest the integration of unstructured mesh discretizations based on Bézier triangles [13], tetrahedra [31] and quadrilaterals [5,16]. For isogeometric analysis, a useful tool for generating parametrized geometries is the Octave NURBS toolbox [11,25], which is also compatible with MATLAB. An example for generating the plate with a hole geometry is given in Listing 11.1. In this file, we generate the model for a given radius of the hole R and a plate size L by creating the boundary segments, and parametrizing the interior with the Coon patch method [9] via the nrbcoons command. Degree elevation (p-refinement) to quadratic
584
Extended Finite Element and Meshfree Methods
NURBS and knot-insertion (h-refinement) are handled by the nrbdegelev and nrbkntins commands respectively. We refer to the documentation of the NURBS toolbox for more details on these commands. The plots showing the control points and the elements for a mesh with 10 × 10 knot-spans are given in Fig. 11.1. % s c r i p t p l o t P l a t e W i t h H o l e .m %g e n e r a t e s the geometry f o r a q u a r t e r p l a t e with a hole % u s e s t h e NURBS t o o l b o x ( h t t p s : / / o c t a v e . s o u r c e f o r g e . i o / n u r b s / i n d e x . h t m l ) addpath ( ' . / nurbs −1.3.13/ i n s t ' ) R = 1; %r a d i u s of the hole L = 4 ; % s i z e o f t h e modeled p l a t e % c r e a t e t h e c u r v e d b o u n d a r y and t h e o t h e r f o u r l i n e a r b o u n d a r y s e g m e n t s a r c = n r b c i r c (R, [ 0 , 0 ] , p i / 2 , p i ) ; s i d e B o t t o m = n r b l i n e ([ −R, 0 ] , [ − L , 0 ] ) ; s i d e R i g h t = n r b l i n e ( [ 0 , R] , [ 0 , L ] ) ; pointUpRight = [ 0 ; L ; 0 ] ; p o i n t U p L e f t = [ −L ; L ; 0 ] ; p o i n t D o w n L e f t = [ −L ; 0 ; 0 ] ; knotUpLeft = [ 0 , 0 , 1/2 , 1 , 1 ] ; s i d e U p L e f t = nrbmak ( [ p o i n t U p R i g h t , p o i n t U p L e f t , p o i n t D o w n L e f t ] , k n o t U p L e f t ) ; %d e g r e e e l e v a t e t h e l i n e a r s i d e and c r e a t e a 2− e l e m e n t v o l u m e t r i c model sideBottom = nrbdegelev ( sideBottom , 1) ; plateModel = nrbcoons ( arc , sideUpLeft , sideRight , sideBottom ) ; %p e r f o r m k n o t i n s e r t i o n s t o r e f i n e t h e model newKnotsU = [ 0 . 1 : 0 . 1 : 0 . 4 , 0 . 6 : 0 . 1 : 0 . 9 ] ; newKnotsV = [ 0 . 1 : 0 . 1 : 0 . 9 ] ; p l a t e M o d e l = n r b k n t i n s ( p l a t e M o d e l , { newKnotsU , newKnotsV } ) ; % p l o t t h e model and t h e c o n t r o l p o i n t s figure subplot (1 ,2 ,1) , n r b c t r l p l o t ( plateModel ) view ( 2 ) subplot (1 ,2 ,2) , nrbkntplot ( plateModel ) view ( 2 )
Listing 11.1: Matlab code for generating the plate with a hole model. In the framework of enriched methods, additional structures are needed to represent the sets of enriched nodes and the type of enrichment. A unified way to implement this concept is to define an array corresponding to each node that determines its set of enrichment functions. For the functions that are not enriched, this array will consist of just the unity function {ξ1 (x) = 1}. For other functions, for example the Heavisideenriched basis, the array will also include corresponding enrichments, e.g. {ξ1 (x) = 1, ξ2 (x) = H (x)}. Then the derivatives of the basis functions of the form φi ξji can be computed by the product rule from the derivatives of φi , i = 1, . . . , N, and ξji , j = 1, . . . , ni where N is the number of (unenriched)
Implementation details
585
Figure 11.1 Plate with a hole model generated by Listing 11.1, showing the control points on the left and the refined mesh on the right.
global shape functions and ni is the number of enrichments (including the unity) for the ith basis function. Finally, in the case level-sets are used to track the strong or weak discontinuities, a discretization for the level-sets is needed, usually in terms of the global basis functions. We refer to the included codes and Sections 1.4 and 8.3 for additional implementation details.
11.1.2 Processing The processing stage of a finite element program consists of constructing and solving a linear system which is obtained by assembling the stiffness matrices and load vector and by imposing the boundary conditions. The assembly stage is typically conducted as an outer loop over the elements (or integration cells in the case of meshfree methods) and an inner loop over the quadrature Gauss points. In the case of a MATLAB implementation, it is preferable that the values of all the shape functions with support at a particular Gauss point are output as a vector by the basis function subroutine. Then the process of assembling the local stiffness matrices can be vectorized, resulting in a significant performance gain. Another MATLAB-specific speedup is the use of the so-called “triplet arrays” for the assembly of the sparse-format stiffness matrices. The idea is to replace matrix update statements of the form K(rowIndex,colIndex)=K(rowIndex,colIndex)+localK with statements like the ones in Listing 11.2. The reason is that the sparse matrices are stored in the compressed sparse column (CSC) format [15], for which adding additional non-zero
586
Extended Finite Element and Meshfree Methods
entries requires expensive resorting operations. Therefore, it is much more efficient to store the non-zero entries in triplet arrays and build the sparse matrix at the end. A drawback of this procedure is that the triplet arrays consume more memory than the sparse matrix, particularly when there are many repeated row and column indices whose corresponding value needs to be summed. For large matrices (with more than 100,000 degrees of freedom) or when there is a lot of overlap between the basis functions, it may be required to break up the assembly into subdomains and sum up the assembled stiffness matrices. This procedure can also be helpful in the parallelization of the assembly subroutine. iCounter = 0 % A l l o c a t e s p a c e f o r t h e a r r a y s rowSet , c o l S e t , v a l S e t %. . . f o r i n d e x E l e m = 1 : numElements % s e t c o l I n d e x and r o w I n d e x from t h e c o n n e c t i v i t y a r r a y s and %a s s e m b l e t h e l o c a l s t i f f n e s s m a t r i x by l o o p i n g o v e r G a u s s p o i n t s %. . . numEnt = l e n g t h ( r o w I n d e x ) ; r o w S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e p m a t ( rowIndex , 1 , numEnt ) ; c o l I n d e x S e t = r e p m a t ( c o l I n d e x , numEnt , 1 ) ; c o l S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e s h a p e ( c o l I n d e x S e t , 1 , numEnt ^ 2 ) ; v a l S e t ( i C o u n t e r +1: i C o u n t e r+numEnt ^ 2 ) = r e s h a p e ( l o c a l K , 1 , nument ^ 2 ) ; i C o u n t e r = i C o u n t e r + numEnt ^ 2 ; end K = s p a r s e ( rowSet , c o l S e t , v a l S e t , dimK , dimK ) ;
Listing 11.2: Snippet of Matlab code for triplet matrix assembly. For the imposition of the boundary conditions, we distinguish between Dirichlet (e.g. prescribed displacement) and Neumann (e.g. prescribed traction) boundary conditions. The Dirichlet boundary conditions are imposed by changing the stiffness matrix and right-hand side vector, while for Neumann boundary conditions only the right-hand side vector is changed. In the case of finite element method with Lagrange polynomials, a common procedure is to the zero-out the rows of the stiffness matrix corresponding to the degrees of freedom on the Dirichlet boundaries, set the diagonal entries in these rows to one and set the corresponding entries in the righthand side vector to the prescribed value. This has the drawback that the resulting stiffness matrix is no longer symmetric, which makes the process of solving the linear system less efficient. Symmetry can be restored by condensing out the prescribed degrees of freedom. To give a simple example, suppose we have a matrix with 4 degrees of freedom, of which the second one is given by x2 = u¯ , with u¯ known. The initial linear symmetric system
587
Implementation details
can be written as: ⎡
k11 ⎢k ⎢ 12 ⎢ ⎣k13 k14
k12 k22 k23 k24
k13 k23 k33 k34
⎤⎡ ⎤
⎡ ⎤
k14 x1 f1 ⎢ ⎥ ⎢ ⎥ k24 ⎥ ⎥ ⎢x2 ⎥ ⎢f2 ⎥ ⎥⎢ ⎥ = ⎢ ⎥ k34 ⎦ ⎣x3 ⎦ ⎣f3 ⎦ k44 x4 f4
(11.1)
After setting x2 = u¯ , we now have: ⎡
⎤⎡ ⎤
⎡ ⎤
k11 k12 k13 k14 x1 f1 ⎢0 ⎥ ⎢x ⎥ ⎢ u¯ ⎥ 1 0 0 ⎢ ⎥ ⎢ 2⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ = ⎢ ⎥ ⎣k13 k23 k33 k34 ⎦ ⎣x3 ⎦ ⎣f3 ⎦ k14 k24 k34 k44 x4 f4
(11.2)
To restore the symmetry, we substitute x2 = u¯ in the equations corresponding to the rows 1, 3 and 4 and move the known quantities to the right hand side: ⎡
k11 ⎢0 ⎢ ⎢ ⎣k13 k14
⎤⎡ ⎤
⎡
⎤
0 k13 k14 x1 f1 − k12 u¯ ⎢ ⎥ ⎢ ⎥ ⎥ 1 0 0 ⎥ ⎢x2 ⎥ ⎢ u¯ ⎥ ⎥⎢ ⎥ = ⎢ ⎥ 0 k33 k34 ⎦ ⎣x3 ⎦ ⎣f3 − k23 u¯ ⎦ f4 − k24 u¯ 0 k34 k44 x4
(11.3)
This procedure can naturally be applied for the case of multiple prescribed degrees of freedom. A MATLAB implementation of this algorithm is given in Listing 11.3. In this code, additionally the diagonal entries on the stiffness matrix corresponding to the Dirichlet degrees of freedom and the entries in the load vector are multiplied by the mean value of the diagonal of the original matrix to improve conditioning. We note that in the case of IGA or mesh methods, the shape functions on the boundary do not, in general, satisfy the Kronecker-Delta property. However, the B-Spline and NURBS functions, as well as some meshfree basis such as the local maximum entropy (LME) basis functions on convex domains or generalized meshfree (GMF) approximations [30] do satisfy a weak Kronecker-delta property. In this case, the shape functions corresponding to the interior nodes (or control points) do not have support on the boundary, however the value of the basis functions on the boundary may not reach 1 at the corresponding nodes as in the strong Kroneckerdelta case. This requires an interpolation or projection of the boundary date on the approximation space spanned by the boundary basis functions. In the case where interior shape functions also have support on the bound-
588
Extended Finite Element and Meshfree Methods
f u n c t i o n [ kk , f f ]= a p p l y D i r i c h l e t S y m ( kk , f f , b c d o f , b c v a l ) % A p p l y t h e D i r i c h l e t c o n s t r a i n t s t o m a t r i x e q u a t i o n [ kk ] { x }={ f f } % w h i l e o u t p u t t i n g a s y m m e t r i c m a t r i x kk % Input : % kk − s y s t e m m a t r i x b e f o r e a p p l y i n g c o n s t r a i n t s % f f − system vector before applying c o n s t r a i n t s % bcdof − a vector containging c o n s t r a i n e d d . o . f % bcval − a vector containing constrained value % Output : % kk − t h e new s y s t e m m a t r i x a f t e r a p p l y i n g t h e c o n s t r a i n t s % f f − t h e new l o a d v e c t o r a f t e r a p p l y i n g t h e c o n s t r a i n t s bcwt=mean ( d i a g ( kk ) ) ; % a m e a s u r e o f t h e a v e r a g e s i z e o f an e l e m e n t i n K f f = f f −kk ( : , b c d o f ) ∗ b c v a l ' ; % m o d i f y t h e force vector f f ( b c d o f ) = bcwt ∗ b c v a l ; kk ( b c d o f , : ) =0; % z e r o o u t t h e r o w s and c o l u m n s o f t h e K m a t r i x kk ( : , b c d o f ) = 0 ; kk ( b c d o f , b c d o f ) =bcwt ∗ s p e y e ( l e n g t h ( b c d o f ) ) ; % p u t o n e s ∗ bcwt on t h e d i a g o n a l
Listing 11.3: Matlab code for imposing Dirichlet boundary conditions. ary (such as in the case of MLS or RKP meshfree methods), then other techniques such as Lagrange multipliers are required. After the linear system is changed to take into account the boundary conditions, the solution vector can be obtained by either a direct or iterative solver. Direct solvers (such as the MATLAB backslash “\” operator) tend to be more robust for ill-conditioned linear systems, but the memory requirements can become very high, particularly for 3D problems with a lot of degrees of freedom. The reason is due to the “fill-in” which occurs when the zero entries within the bandwidth of a matrix become non-zero when the matrix is factorized. Iterative (Krylov subspace) solvers rely on matrix vector multiplications and can be used to solve very large linear systems efficiently. Examples of iterative solvers are pre-conditioned conjugate gradient (pcg command in MATLAB) which is efficient for symmetric positive definite matrices or the generalized minimal residual method (GMRES, gmres command in MATLAB) for general matrices. These solvers work better when a pre-conditioner is applied to improve the convergence. Standard preconditioners are incomplete Cholesky factorizations (ichol command in MATLAB) for symmetric positive definite matrices and the incomplete LU factorization (ilu command in MATLAB) for general sparse matrices. The idea of incomplete factorization is that only the non-zero entries in the original matrix can be non-zero in the factorization, the other entries are set to zero. Other preconditioners, such as Jacobi preconditioners which
Implementation details
589
perform a diagonal scaling of the matrices, or multi-grid preconditioners can be used as well. The solution of linear systems also appears in nonlinear analysis, where usually Newton-Raphson iterations are performed to solve a linearized model. In this case, on the left-hand side we have the “tangent-stiffness” matrix, as described e.g. in Section 3.6.1.3 for hydro-mechanical coupled problems, and on the left-hand side we have a residual vector. Solving linear systems is also required for time-dependent problems with implicit timeintegration (where a system involving the stiffness matrix is solved) or with explicit time-integration where the consistent mass matrix is used. Therefore, for large simulations, the use of efficient solvers or of accurate lumped mass methods is of great importance in keeping the computational times within acceptable limits.
11.1.3 Post-processing In the post-processing stage, various quantities of interest can be computed and displayed based on the solution vector. A correct and detailed visualization of the output date is very important, particularly for identifying any errors in the model such as incorrect imposition of boundary conditions. In the case of XFEM or basis functions which do not satisfy the (strong) Kronecker-delta property, the displacements cannot be simply inferred from the values of the solution vector. Instead, the basis functions (including the enrichments) need to be evaluated at the chosen points (e.g. element vertices) and the values multiplied by the corresponding coefficients of the solution. The strain and stress values are calculated in the same way, by considering the derivatives of the basis functions and of the enrichments. For meshless methods, typically a visualization mesh needs to be constructed for example from a Delaunay triangulation of the domain. Moreover, to precisely show the discontinuities in a fracture simulation, the visualization mesh should be conforming to the crack path. For many examples, we use output files in the VTK (Visualization Toolkit) format. These files can be opened with open-source viewers such as Paraview, and are particularly convenient for 3D visualization as one can easily manipulate the graphs and interactively plot the deformations, stresses and other quantities of interest on structured or unstructured meshes. Highlevel output subroutines such as vtkwrite1 are available to make the output of scalar or vector quantities relatively straightforward. 1 https://github.com/joe-of-all-trades/vtkwrite.
590
Extended Finite Element and Meshfree Methods
For fracture propagation problems, a crack increment is computed according to some failure as well as propagation magnitude and angle criteria. In XFEM (and other enriched methods), stress intensity factors are often used for the computation of these criteria. They can be calculated in 2D by the interaction integral methods [26], which are derived from the J-integral [22], as discussed in Sections 5.9.1.4 and 7.2.2.5. This method considers a contour around the crack tip and is converted to a domain integral for improved integration accuracy. Since the resulting code is quite lengthy, we refer to the included examples and the documentation for more details. Finally, even though the idea of enriched methods is that the mesh is independent of the crack, often times the results can be improved by considering a finer mesh around the crack tips. In the framework of hierarchical splines, this can be conveniently implemented with a recovery-type error estimator, as detailed in [2]. This type of estimator is quite effective at refining in regions with high stress concentrations. Another aspect to be considered is the selection of a marking scheme and the amount of elements being marked at each refinement step. When the mesh is coarse, the error can be quite large even relatively far from the regions with non-smooth solution such as the crack tip. As the mesh is refined and the singular regions are better resolved, the areas with large error can shrink considerably. Therefore, it is more economical in terms of degrees of freedom to refine as few elements as possible at each step, however this leads to more refinement steps which increases the computational cost. A strategy that works well and has been implemented in [2] is the Dörfler marking scheme [12], which refines elements are in the 100θ percentile of the largest error. For problems with singularity, we typically choose θ = 0.5, while for problems with smoother solutions (but which can still benefit from adaptivity) θ = 0.75 is more efficient.
11.2. Numerical examples In this section we provide a few numerical examples of benchmark problems in XFEM, including that of a coupled hydro-mechanical model. Further examples and documentation, including 3D, plate and shell models can be found in the open-source repository IGA(X)FEM at https:// sourceforge.net/projects/cmcodes.
591
Implementation details
Figure 11.2 Schematic illustration of mechanical loading.
Table 11.1 Material parameters for model under pure mechanical loading. Young’s modulus E = 25.85 GPa Poisson’s ratio ν = 0.18 uniform stress field σ1 = 1.0e7 Pa crack length 2a = 2.15 m
11.2.1 Crack propagation angle This section presents a square domain with an inclined center crack under pure mechanical loading (Fig. 11.2). Hereafter, the initial propagation angles of the induced crack with different pre-defined angle β and loading factor α are studied. The material parameters are list in Table 11.1. Fig. 11.3 shows the stress profiles of domains with pre-defined angle β ranging from 10 to 80 degree and loading factor α is set to 0. The comparison between numerical calculation and theoretical model [1] in terms of the propagation angle of the induced crack with different loading ratio α is shown in Fig. 11.4A. Furthermore, the convergence study (case β = 10) of stiffness intensity factor K1 of mode I respect to total number of elements is presented in Fig. 11.4B.
11.2.2 Hydro-mechanical model with center cracks This section presents the domain with inner fractures under fluid flux loading. The geometry and boundary conditions are shown in Fig. 11.5. The material parameters are listed in Table 11.2. Domain with single and five
592
Extended Finite Element and Meshfree Methods
Figure 11.3 Stress profile in x direction with different pre-defined angle β s.
randomly distributed inner cracks are studied (Fig. 11.6). Fig. 11.6 clearly indicates that the existing crack has great impact on deformation, pressure profiles of the whole domain.
11.2.3 Hydro-mechanical model with edge crack This section simulates the rupture of a plate under plane-strain conditions. As shown in Fig. 11.7A, a pre-existing fracture is located at the symmetry
Implementation details
593
Figure 11.4 (A) Comparison results of the propagation angle of the induced crack between numerical calculation and theoretical model with different loading ratio α , solid lines represent the result from theoretical model. (B) Convergence study of stiffness intensity factor K1 of mode I respect to total element number.
Figure 11.5 Schematic illustration of geometry and boundary conditions. p represents pressure, q represents fluid flux. Table 11.2 Material parameters for hydro-mechanical model with inner cracks. Young’s modulus E = 25.85 GPa Poisson’s ratio ν = 0.18 Initial porosity n = 0.2 Intrinsic permeability k = 2.78e-21 m2 Biot’s coefficient α=1 Solid bulk modulus Ks = 13.46 GPa Water bulk modulus Kf = 0.2 GPa ρ = 1000 kg/m3 Water density Water dynamic viscosity μw = 5e-4 N/m2 s Irreducible saturation Sirr = 0.0 Reference pressure Rref = 18.6 MPa Porous index l = 0.4396
594
Extended Finite Element and Meshfree Methods
Figure 11.6 Displacement, pressure and stress profile for domain with single and five cracks.
axis of the plate. A fixed vertical velocity v = 2.35 × 10−2 µm/s is prescribed in an opposite direction at the bottom and at the top of the plate (tensile loading). All boundaries of the plate are assumed to be impervious. The induced pressure profile is presented in Fig. 11.7B.
11.2.3.1 Hydro-mechanical model with edge crack under fluid flux loading This subsection follows an experimental setup [8] (Fig. 11.8). The setup of this experiment can be summarized as follows: a square sample (Fig. 11.8A) with two ports, high pore pressure is introduced to the lower port region by injecting pressurized fluid and maintain the pressure, hereafter a wedge
Implementation details
595
Figure 11.7 (A) Schematic illustration of geometry and boundary conditions, (B) pressure profile.
Figure 11.8 (A) Schematic illustration of geometry and boundary conditions, (B) experimental results.
is driven into the pre-shaped concave edge to induce fracture and fracture propagation. The crack propagation path is then identified. Further more, the impact of the permeability of the solid on the crack propagation path is also studied (Fig. 11.9). The fracture development path, displacement of the domain and stress profiles are shown in Fig. 11.9 for porous medias (under identical loading) with different intrinsic permeabilities (k1, k2, k3). Fig. 11.9D-F clearly demonstrates that the fracture paths are attracted to the fluid pressurized port. Furthermore, we found that under identical injection rate, porous media with lower permeability (e.g. k1) will induce
596
Extended Finite Element and Meshfree Methods
Figure 11.9 Displacement (m) plot of porous media with the intrinsic permeabilities (A) k1 = 2.78e-18 m2 ; (B) k2 = 1.52e-17 m2 ; (C) k3 = 2.78e-17 m2 and (D)-(F) are the stress (Pa) plots in cylinder coordinate accordingly.
Implementation details
597
higher stress/pressure state around the injection port comparing to porous media with higher permeability (e.g. k3), which, in return influences the propagation direction of the fracture.
References [1] T. Anderson, Fracture Mechanics: Fundamentals and Applications, fourth edition, CRC Press, 2017. [2] C. Anitescu, M.N. Hossain, T. Rabczuk, Recovery-based error estimation and adaptivity using high-order splines over hierarchical t-meshes, Computer Methods in Applied Mechanics and Engineering 328 (2018) 638–662. [3] D. Benson, Y. Bazilevs, M. Hsu, T. Hughes, Isogeometric shell analysis: the Reissner– Mindlin shell, Computer Methods in Applied Mechanics and Engineering 199 (5) (2010) 276–289. [4] D. Benson, Y. Bazilevs, M.-C. Hsu, T. Hughes, A large deformation, rotation-free, isogeometric shell, Computer Methods in Applied Mechanics and Engineering 200 (13) (2011) 1367–1378. [5] M. Bercovier, T. Matskewich, Smooth Bézier Surfaces Over Unstructured Quadrilateral Meshes, Springer International Publishing, 2017. [6] S. Bordas, V. Nguyen, C. Dunant, H. Nguyen-Dang, A. Guidoum, An extended finite element library, International Journal for Numerical and Analytical Methods in Engineering 71 (6) (2007) 703–732. [7] M.J. Borden, M.A. Scott, J.A. Evans, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of NURBS, International Journal for Numerical Methods in Engineering 87 (1–5) (2011) 15–47. [8] M. Bruno, F. Nakagawa, Pore pressure influence on tensile fracture propagation in sedimentary rock, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28 (4) (1991) 261–273. [9] S.A. Coons, Surfaces for Computer-Aided Design of Space Forms, Technical report, Cambridge, MA, USA, 1967. [10] L. Dalcin, N. Collier, P. Vignal, A. Côrtes, V. Calo, PetIGA: a framework for highperformance isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 308 (2016) 151–181. [11] C. de Falco, A. Reali, R. Vázquez, GeoPDEs: research tool for isogeometric analysis of PDEs, Advances in Engineering Software 42 (12) (2011) 1020–1034. [12] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM Journal on Numerical Analysis 33 (3) (1996) 1106–1124. [13] L. Engvall, J.A. Evans, Isogeometric triangular Bernstein–Bézier discretizations: automatic mesh generation and geometrically exact finite element analysis, Computer Methods in Applied Mechanics and Engineering 304 (2016) 378–407. [14] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331. [15] J. Gilbert, C. Moler, R. Schreiber, Sparse matrices in MATLAB: design and Implementation, SIAM Journal on Matrix Analysis and Applications 13 (1) (1992) 333–356. [16] M. Kapl, G. Sangalli, T. Takacs, An isogeometric C1 subspace on unstructured multipatch planar domains, Computer Aided Geometric Design 69 (2019) 55–75.
598
Extended Finite Element and Meshfree Methods
[17] Y. Lai, Y.J. Zhang, L. Liu, X. Wei, E. Fang, J. Lua, Integrating CAD with Abaqus: a practical isogeometric analysis software platform for industrial applications, in: HighOrder Finite Element and Isogeometric Methods 2016, Computers & Mathematics with Applications 74 (7) (2017) 1648–1660. [18] V.P. Nguyen, C. Anitescu, S.P. Bordas, T. Rabczuk, Isogeometric analysis: an overview and computer implementation aspects, Mathematics and Computers in Simulation 117 (2015) 89–116. [19] V.P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshfree methods: a review and key computer implementation aspects, Mathematics and Computers in Simulation 79 (3) (2008) 763–813. [20] M. Occelli, T. Elguedj, S. Bouabdallah, L. Morançay, LR B-Splines implementation in the Altair RadiossTM solver for explicit dynamics IsoGeometric Analysis, Advances in Engineering Software 131 (2019) 166–185. [21] M. Pauletti, M. Martinelli, N. Cavallini, P. Antolin, Igatools: an isogeometric analysis library, SIAM Journal on Scientific Computing 37 (4) (2015) C465–C496. [22] J.R. Rice, A path independent integral and the approximate analysis of stress concentration by notches and cracks, Journal of Applied Mechanics 35 (1968) 379–386. [23] D. Rypl, B. Patzák, From the finite element analysis to the isogeometric analysis in an object oriented computing environment, in: CIVIL-COMP, Advances in Engineering Software 44 (1) (2012) 116–125. [24] H.Si. TetGen, A Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathematical Software 41 (2) (2015) 1–36. [25] M. Spink, D. Claxton, C. de Falco, R. Vázquez, The NURBS toolbox, http://octave. sourceforge.net/nurbs/index.html. [26] M. Stern, E. Becker, R. Dunham, A contour integral computation of mixed-mode stress intensity factors, International Journal of Fracture 12 (3) (1976) 359–368. [27] R.L. Taylor, FEAP - Finite Element Analysis Program, 2014. [28] N. Valizadeh, Y. Bazilevs, J. Chen, T. Rabczuk, A coupled IGA–meshfree discretization of arbitrary order of accuracy and without global geometry parameterization, Computer Methods in Applied Mechanics and Engineering 293 (2015) 20–37. [29] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Computer Mathematics and Its Applications 72 (3) (2016) 523–554. [30] C. Wu, C. Park, J. Chen, A generalized approximation for the meshfree analysis of solids, International Journal for Numerical Methods in Engineering 85 (6) (2011) 693–722. [31] S. Xia, X. Qian, Isogeometric analysis with Bézier tetrahedra, in: Special Issue on Isogeometric Analysis: Progress and Challenges, Computer Methods in Applied Mechanics and Engineering 316 (2017) 782–816.