Exponential finite element shape functions for a phase field model of brittle fracture

478 EXPONENTIAL FINITE ELEMENT SHAPE FUNCTIONS FOR A PHASE FIELD MODEL OF BRITTLE FRACTURE CHARLOTTE KUHN∗AND RALF M ¨ ULLER† Institute of Applied Mechanics Technische Universit¨at Kaiserslautern P.O.B. 3049, 67653 Kaiserslautern http://mechanik.mv.uni-kl.de/ ∗e-mail: [email protected] †e-mail: [email protected] Key words: Phase Field Model, Fracture, Finite Elements, Exponential Shape Functions Abstract. In phase field models for fracture a continuous scalar field variable is used to indicate cracks, i.e. the value 1 of the phase field variable is assigned to sound material, while the value 0 indicates fully broken material. The width of the transition zone where the phase field parameter changes between 1 and 0 is controlled by a regularization parameter. As a finite element discretization of the model needs to be fine enough to resolve the crack field and its gradient, the numerical results are sensitive to the choice of the regularization parameter in conjunction with the mesh size. This is the main challenge and the computational limit of the finite element implementation of phase field fracture models. To overcome this limitation a finite element technique using special shape functions is introduced. These special shape functions take into account the exponential character of the crack field as well as its dependence on the regularization length. Numerical examples show that the exponential shape functions allow a coarser discretization than standard linear shape functions without compromise on the accuracy of the results. This is due to the fact, that using exponential shape functions, the approximation of the surface energy of the phase field cracks is impressively precise, even if the regularization length is rather small compared to the mesh size. Thus, these shape functions provide an alternative to a numerically expensive mesh refinement. 1 INTRODUCTION Variational formulations of brittle fracture as suggested by Francfort and Marigo [1] overcome some of the limitations of classical Griffith theory. However, a direct discretization of such fracture models is faced with significant technical difficulties. A regularized approximation by means of Γ-convergence as presented by Bourdin [2] offers a new 1 XI International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XI E. Oñate, D.R.J. Owen, D. Peric and B. Suárez (Eds) 479 Charlotte Kuhn and Ralf M¨uller perspective towards the computational implementation of the model. The core of the regularization is the approximation of the total energy functional, in which a continuous scalar field variable is introduced to indicate cracks, i.e. the value of 1 is assigned to sound material and a value of 0 indicates fracture. With this crack field the regularized model resembles a phase field model for fracture, where additionally a Ginzburg–Landau type equation is used to describe the evolution of the crack field and cracking is addressed as a phase transition problem. Similar phase field fracture models have been introduced e.g. in [3, 4, 5, 6, 7, 8]. Differing in technical details all of these models introduce a regularization length which controls the width of the transition zone where the crack field interpolates between broken and unbroken material; i.e. the smaller the regularization parameter, the smaller the transition zone and the higher the gradients of the crack field in the vicinity of the cracks. Numerical implementations are faced with the difficulty that the spacial discretization has to be fine enough to resolve these high gradients of the crack field, which leads to high computational costs for small values of the regularization parameter. On the other hand the regularization length needs to be chosen sufficiently small in conjunction with the global geometric dimension of the sample in order to get reasonable results. The most common approach to meet the requirement for a sufficiently fine resolution on the one hand and to keep the computation time within bounds on the other hand are adaptive mesh refinement techniques as used e.g. in [9], where the mesh is only refined where it is necessary, i.e. in the vicinity of a crack. Another approach to increase the efficiency of the computations was introduced in [5], where Fourier transforms are used to solve the linear part of the problem. However, this technique restricts the simulations to problems with periodic boundary conditions. In this work we follow a different approach which is inspired by [10], where exponential finite element (FE) shape functions are introduced as an alternative to an extensive mesh refinement in the simulation of extrusion processes. These special shape functions qualitatively capture the shape of the solution and thus allow a much coarser discretization than the standard discretization using linear shape functions. In contrast to the simulation of the extrusion process, where the exponential shape functions are used in one distinct direction only, the discretization of the crack field in a two dimensional setting requires an extension of the concept to the full 2d case. 2 A PHASE FIELD MODEL FOR FRACTURE 2.1 Governing Equations The present phase field model of fracture is based on a regularized version of the variational formulation of brittle fracture by [1] which was introduced in [11]. The core of the regularization is the approximation of the total energy of a cracked linear elastic 2 Charlotte Kuhn and Ralf M¨uller body Ω, with the stiffness tensor Cand the cracking resistance Gc, by the functional E(ε, s) = ZΩ ψ(ε, s)dV =ZΩ 1 2(s2+η)ε·(Cε) |{z } =ψe(ε,s) +Gc1 4ǫ(1 −s)2+ǫ|∇s|2 |{z } =ψs(s) dV . (1) This energy as well as the energy density ψare functions of the linearized strain tensor ε=1 2(∇u+ (∇u)T), i.e. the symmetric part of the gradient of the displacements uand the continuous scalar crack field s, which takes the value 1, if the material is undamaged, and 0 if there is a crack. The degradation of the elastic energy in the bulk Ee=RΩψedV upon cracking is modeled by the factor (s2+η), where the small positive parameter ηis introduced to obtain an artificial rest stiffness ηCat fully broken state (s= 0) in order to circumvent numerical difficulties. The parameter ǫ, appearing twice in the surface energy Es=RΩψsdV , has the dimension of length and controls the width of the transition zone between broken and unbroken material, where sinterpolates between 0 and 1. If body forces and inertia terms are neglected, the mechanical part of the problem is described by the local balance law for the Cauchy stress tensor σ divσ=0,(2) plus the according boundary conditions σn =t∗ non ∂Ωt, where nis the outer normal vector, and the material law (3) derived from the energy density ψ σ=∂ψ ∂ε= (s2+η)Cε.(3) Interpreting sas order parameter of a phase field model, its evolution in time is assumed to follow a Ginzburg–Landau type evolution equation, where ˙sis proportional to the variational derivative of the energy density ψwith respect to s. ˙s=−M·δψ δs =−Msε·(Cε)−Gc2ǫ∆s+1−s 2ǫ (4) The mobility factor Mis a positive constant, which controls the dissipation in the process zone. For sufficiently large values of Mthe solution of the evolution equation can be considered as stationary. In order to take into consideration the irreversible character of cracking, s(x, t) is fixed to 0 for all future times t > t∗if it becomes 0 at any time t∗. 2.2 Evolution Equation in 1d In a 1d setting, the evolution equation for a stationary ( ˙s= 0) crack field reduces to s′′ −s 4ǫ2=−1 4ǫ2,(5) 480 Charlotte Kuhn and Ralf M¨uller −1 −0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 x/L s ǫ→0 Figure 1: 1d stationary crack field if elastic contributions are neglected. With boundary conditions s(0) = 0 and s′(±∞) = 0 the analytic solution of Eq. (5) is given by s(x) = 1 −exp −|x| 2ǫ(6) Figure 1 illustrates the impact of the regularization length ǫon the crack field s(x). The smaller ǫgets, the higher gradients and curvatures of the solution s(x) appear in the vicinity of the crack at x= 0. The limit ǫ→0 yields a discontinuous function, which is 0 at x= 0 and 1 elsewhere. 3 NUMERICAL IMPLEMENTATION 3.1 Weak Forms Starting point for the FE implementation of the coupled problem of mechanical balance equation (2) and evolution equation (4) are the weak forms of these field equations. With virtual displacements δuand δs, they read ZΩ∇δu·σdV =Z∂Ωt δu·t∗ ndA (7) with prescribed surface traction t∗ non part ∂Ωtof the boundary and ZΩδs ˙s M−∇δs ·q+δs sε: [Cε] + Gc 2ǫ(s−1)dV = 0 (8) with q=−2Gcǫ∇s. The normal flux q·nis assumed to vanish on the boundary ∂Ω. 3.2 Finite Element Discretization In a 2d setting the weak forms of the field equations (7) and (8) are discretized with 4 node quadrilateral elements with 3 degrees of freedom (ux, uy, s) per node. The displacements u, the crack field s, as well as their virtual counterparts δuand δs are approximated 481 Charlotte Kuhn and Ralf M¨uller by shape functions Nu I,Ns I,Nδu I, and Nδs I, which interpolate the respective nodal values ˆ uI, ˆsI,δˆ uI, and δˆsI. Using Voigt–notation - denoted by an underline in the following - the approximations read u= N X I=1 Nu Iˆ uI, s = N X I=1 Ns IˆsI, δu= N X I=1 Nδu Iδˆ uI,and δs = N X I=1 Nδs IδˆsI.(9) Accordingly the approximations of the gradient expressions yield ε= N X I=1 [Bu I]ˆ uI,∇s= N X I=1 [Bs I]ˆsI, δε= N X I=1 [Bδu I]δˆ uI,and ∇δs = N X I=1 [Bδs I]δˆsI,(10) where the derivative matrices [Bu I] =   Nu I,x 0 0Nu I,y Nu I,y Nu I,x  ,[Bs I] = Ns I,x Ns I,y,[Bδu I] =   Nδu I,x 0 0Nδu I,y Nδu I,y Nδu I,x  ,and [Bδs I] = Nδs I,x Nδs I,y(11) are obtained from the derivatives of the shape functions. By a standard argument for finite element approximations the nodal values δˆ uIand δˆsI of the virtual quantities δuand δs drop out of the system of equations, leading to the nodal residuals [RI] = Ru I Rs I=ZΩ [Bδu I]Tσ Nδs I ˙s M−[Bδs I]Tq+Nδs IsεT·(Cε) + Gc 2ǫ(s−1) dV . (12) The time integration of the transient terms is performed with the implicit Euler method. Together with the nonlinear character of the phase field model this yields a nonlinear system of equations, which has to be solved in every time step ∆t. This is done with a Newton–Raphson algorithm, which requires the derivation of the consistent tangent matrix SIJ which has the following structure: SIJ =KIJ +1 ∆tDIJ .(13) The stiffness matrix KIJ and the damping matrix DIJare obtained by derivation of the nodal residuals RIwith respect to the nodal values (ˆ uJ,ˆsJ) and (ˆ ˙uJ,ˆ ˙sJ), respectively. KIJ =ZΩ [Bδu I]T(s2+η)C[Bu J] [Bδu I]T2sCεNs J Nδs I2s(Cε)T[Bu J] 2Gcǫ[Bδs I]T[Bs J]+Nδs IεT·Cε+Gc 2ǫNs J dV (14) 482 Charlotte Kuhn and Ralf M¨uller x y η 21 4 3 ξ h2 h3 h4 h1 2 3 1 4 Figure 2: Node and edge numbering of the quadrilateral element in global (left) and natural coordinates (right) and DIJ =ZΩ"0 0 01 MNδs INs J#dV . (15) If the same shape functions are chosen for the approximation of actual values and the virtual quantities, i.e. Nu I=Nδu Iand Ns I=Nδs I, the system matrix SIJ becomes symmetric. This is due to the fact, that the constitutive law (3) as well as the evolution equation (4) are derived from a potential. Different shape functions however, render a non-symmetric system matrix SIJ . 4 EXPONENTIAL SHAPE FUNCTIONS The standard implementation with 4 node quadrilateral elements makes use of the linear Lagrangian shape functions Nlin I(ξ, η) = 1 4(1 + ξIξ)(1 + ηIη), I = 1, ..., 4 (16) with (ξI, ηI) according to Fig. 2 for all the shape functions Nu I,Nδu I,Ns I, and Nδs Ias well as for the approximation of the geometry in the isoparametric concept x= N X I=1 Nlin Iˆ xI.(17) In [12] it is shown that triangular elements with linear shape functions overestimate the surface energy by a factor f(h/ǫ) = 1 + h/4ǫ , (18) where his the edge length of the elements. As a sufficiently good approximation of the surface energy is crucial in order to obtain reasonable results, this yields the necessity of a 483 484 485 486 Charlotte Kuhn and Ralf M¨uller 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x/L s exp/exp lin/exp lin/lin 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x/L s exp/exp lin/exp lin/lin 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x/L s exp/exp lin/exp lin/lin Figure 5: Solution of the 1d stationary evolution equation with n= 4 (left), n= 8 (middle) and n= 16 elements (right) 100101102103 0 1 2 3 4 5 6 7 number of elements per direction surface energy [ GcL] quadrature: 02×02 Gauß points exp/exp lin/exp lin/lin 100101102103 10−4 10−3 10−2 10−1 100 101 number of elements per direction relative error of surface energy quadrature: 02×02 Gauß points exp/exp lin/exp lin/lin 100101102103 10−4 10−3 10−2 10−1 100 101 number of elements per direction relative error of surface energy quadrature: 10×10 Gauß points exp/exp lin/exp lin/lin Figure 6: Evaluation of the surface energy in Fig. 5. While there is almost no visible error in the solutions with exponential shape functions, the linear shape functions fail to adequately resolve the transition zone even for the smallest tested element size h=L/16. 5.2 Surface Energy of an Edge Notched Sample For the first numerical assessment of the 2d exponential shape functions, the stationary evolution equation is solved on the domain L×Lunder the constraint s(x, y)=0 if (x, y)∈[0, L/2] ×{0}. Again, the regularization length is set to =0.01L, and no mechanical loads are applied. A regular mesh with square elements is used for the discretization. Figure 6 shows an evaluation of the surface energy Esassociated with the computed crack field. Regular meshes within the range of 2 ×2 to 400 ×400 elements were used for the discretization. The results are compared to the error estimate (18) for the triangular elements with linear shape functions (black dotted line). The reference solution Es=0.51017344300 GcLwas computed with standard linear shape functions and a nonuniform mesh with square elements of edge length h=7.1429 ·10−4Lin the vicinity the crack. The performance of the tested linear shape functions is slightly better than it is to be expected from the error estimate. However, especially for discretizations with only few 9