Isogeometric collocation for phase-field fracture models Dominik Schillingera , Michael J. Bordenb , Henryk K. Stolarskia a

b

Department of Civil, Environmental and Geo-Engineering, University of Minnesota, USA Institute for Computational Engineering and Sciences, The University of Texas at Austin, USA

Abstract Phase-field models based on the variational formulation for brittle fracture have recently been shown capable of accurately and robustly predicting complex crack behavior. Their numerical implementation requires costly operations at the quadrature point level, which may include finding eigenvalues and forming tensor projection operators. We explore the application of isogeometric collocation methods for the discretization of second-order and fourthorder phase-field fracture models. We show that a switch from isogeometric Galerkin to isogeometric collocation methods has the potential to significantly speed up phase-field fracture computations due to a reduction of point evaluations. We advocate a hybrid collocationGalerkin formulation that provides a consistent way of weakly enforcing Neumann boundary conditions and multi-patch interface constraints, is able to handle the multiple boundary integral terms that arise from the weighted residual formulation, and offers the flexibility to adaptively improve the crack resolution in the fracture zone. We present numerical examples in one and two dimensions that illustrate the advantages of our approach. Keywords: Phase-field modeling, Fracture, Isogeometric collocation

Corresponding author; Department of Civil, Environmental, and Geo-Engineering, University of Minnesota, 500 Pillsbury Drive S.E., Minneapolis, MN 55455-0116, USA; Phone: +1 612 624-0063; Fax: +1 612 626-7750; E-mail: [email protected]

Preprint submitted to Computer Methods in Applied Mechanics and Engineering

July 8, 2014

Contents 1 Introduction

3

2 Variational formulation of phase-field fracture models 2.1 Variational description of fracture . . . . . . . . . . . . . 2.2 Phase-field approximation of the fracture surface . . . . . 2.3 Additive split of the strain tensor . . . . . . . . . . . . . 2.4 Energy degradation using a phase-field . . . . . . . . . . 2.5 Variational formulation . . . . . . . . . . . . . . . . . . . 2.6 Irreversibility and history-based formulation . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

3 Numerical treatment with isogeometric collocation 3.1 The weighted residual form . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Hybrid collocation-Galerkin discretization in space . . . . . . . . . . . . 3.2.1 Collocation at the Greville points . . . . . . . . . . . . . . . . . 3.2.2 Local switch between Galerkin and collocation . . . . . . . . . . 3.2.3 Weak enforcement of Neumann boundary conditions . . . . . . 3.2.4 Handling more complex Neumann boundary conditions . . . . . 3.3 Discrete forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Phase-field part of the second-order theory . . . . . . . . . . . . 3.3.2 Phase-field part of the fourth-order theory . . . . . . . . . . . . 3.3.3 Elasticity part . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Weak enforcement of interface constraints in multi-patch discretizations 3.5 Staggered solution algorithm . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

4 Analysis in one dimension 4.1 One-dimensional model problem . . . . . . . . . . . . . . . . . . . . . . . 4.2 Undamaged elastic bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Elastic bar with a crack . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Adaptive crack resolution using the hybrid collocation-Galerkin approach 5 Numerical examples in two dimensions 5.1 Single edge notched tension test . . . . . . . . . . . . . . . . . . 5.1.1 Cost estimate based on the number of point evaluations . 5.1.2 Crack resolution vs. mesh size . . . . . . . . . . . . . . . 5.2 Single edge notched shear test . . . . . . . . . . . . . . . . . . . 6 Summary, conclusions and outlook

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

5 5 6 7 8 9 10

. . . . . . . . . . . .

11 11 12 13 14 14 14 15 15 16 17 18 19

. . . .

20 20 21 22 25

. . . .

26 26 27 29 31 34

2

1. Introduction Phase-field fracture models introduced by Bourdin et al. [1] are based on the variational approach to brittle fracture due to Francfort and Marigo [2] where the solution to the fracture problem is found as the minimizer of a global energy functional. The basic idea of phase-field fracture is to represent cracks by a continuous scalar field that has a value of one away from the crack and is zero at the crack location. The phase-field serves as a multiplication factor to tensile energy components such that it locally penalizes the capability of the material to carry tensile stress at the crack location. In this sense, the phase-field idea is conceptually very similar to the fictitious domain approach applied e.g. in the finite cell method [3–8]. The diffusiveness of the crack approximation, i.e. the local slope of the phase-field at the crack location, is controlled by a length-scale parameter. From a mathematical point of view, the concept of Γ-convergence can be used to show that the phase-field model converges to the correct minimizing solution as the length-scale parameter goes to zero [9–11]. From a numerical point of view, phase-field fracture models open the door for computational analysis of complex crack patterns in three dimensions. While many established approaches that rely on the introduction of discontinuities into the displacement field by remeshing or by enriching basis functions (see for example [12–18]) have shown much success in modeling two-dimensional fracture problems, these methods have proven difficult to be extended to three-dimensional problems. The diffusive approximation of the crack by a continuous phase-field does not require the introduction of discontinuities. Cracks can be represented independently from the mesh and its topology by the solution of an additional differential equation. Moreover, the reformulation of the fracture problem as a system of partial differential equations that completely determine the evolution of cracks eliminates the need for any adhoc rules or conditions to determine crack nucleation, propagation, or bifurcation. The phase-field fracture approach has proven to accurately and robustly capture crack behavior in two and three dimensions for quasi-static fracture [19–23], dynamic crack propagation [21, 24–28], at finite strains [29], for fracture in piezo- and ferroelectric materials [30–32] and for cohesive fracture [33, 34]. Isogeometric analysis (IGA) was introduced by Hughes and coworkers [35, 36] to bridge the gap between computer aided geometric design (CAD) and finite element analysis. Its core idea is to use the same smooth and higher-order basis functions, e.g. non-uniform rational B-splines (NURBS), for the representation of both the geometry in CAD and the approximation of solution fields in analysis. The initial motivation for IGA was to simplify the cost-intensive mesh generation process required for standard finite elements and to support a more tightly connected interaction between CAD and finite element tools. Since then IGA has developed into an innovative computational mechanics technology, offering a range of new perspectives and opportunities that go far beyond the geometric point of view of analysis. One such opportunity is isogeometric collocation [37, 38], which emerged from the search for more efficient quadrature rules for spline discretizations [39–41]. In contrast to Galerkin-type formulations, collocation is based on the discretization of the strong form of the underlying partial differential equations, which requires basis functions of sufficiently high order and smoothness. Consequently, the use of IGA for collocation suggests itself, 3

since spline functions can be readily adjusted to any order of polynomial degree and continuity required by the differential operators at hand. The major advantage of isogeometric collocation is the minimization of the computational effort with respect to quadrature, since for each degree of freedom only one point evaluation at a so-called collocation point is required. This constitutes a significant advantage for applications where the efficiency and success of an analysis technology is directly related to the cost of quadrature. Isogeometric collocation has been applied successfully in elasticity and explicit elastodynamics [42], for the development of locking-free structural elements [43, 44], in contact and plasticity problems [45], and for the Cahn-Hilliard phase-field model [46]. In this paper, we start to explore the application of isogeometric collocation to phase-field models of brittle fracture, considering the second-order phase-field formulations developed by Miehe, Welschinger and Hofacker [23, 26] and Borden et al. [21], and the fourth-order phase-field formulation recently introduced by Borden et al. [24]. In addition to common quadrature point operations such as evaluating basis functions or computing matrix-matrix products, their numerical implementation requires a range of costly operations at the quadrature point level, such as finding eigenvalues, taking derivatives of tensor-valued functions, and forming tensor projection operators. Therefore, the formation and assembly of stiffness matrices constitutes a significant portion of the total computational cost. So far Galerkin methods have been the dominant tool for the discretization of phase-field fracture models, which require more than one quadrature point per basis function. In particular for higherorder isogeometric Galerkin discretizations, the number of quadrature point evaluations can be overwhelming (see e.g. [38]). Isogeometric collocation has the potential to significantly speed up phase-field fracture computations by reducing the number of point evaluations. At the same time, it maintains the favorable approximation properties of higher-order smooth basis functions, which on a per-degree-of-freedom basis exhibit increased accuracy and robustness in comparison to standard C 0 basis functions [5, 47, 48] and superior spectral accuracy in the higher modes [49–51]. Moreover, the higher-order continuity of spline basis functions allows for the direct numerical treatment of higher-order phase-field fracture models. The fourth-order phasefield model of Borden et al. [24] leads to higher regularity in the exact phase-field solution, which has been shown to improve the convergence rate of the numerical solution. The present work is organized as follows: In Section 2, we provide a brief review of the second- and fourth-order phase-field models of brittle fracture and their variational formulation based on an incremental energy balance principle. We then use integration by parts to transform these expressions into the weighted residual form. In Section 3, we describe a hybrid collocation-Galerkin method that discretizes the weighted residual form using smooth spline basis functions. Away from the boundary and the fracture zone we use Dirac δ distributions at the Greville collocation points as test functions, which leads to discrete collocation equations. Following De Lorenzis et al. [45], we use splines as test functions in the Galerkin sense adaptively along patch boundaries and crack paths. The resulting thin layers of discrete Galerkin equations provide a consistent and stable way to weakly enforce Neumann boundary conditions and multi-patch interface constraints. Galerkin equations at the boundary also provide the flexibility to account for any boundary terms that appear in 4

the weighted residual form of higher-order phase-field models. In addition, it turns out that for a given mesh size Galerkin equations are able to accurately resolve smaller length-scale parameters than collocation equations. Using thin layers of discrete Galerkin equations in the fracture zone thus enables us to achieve a given accuracy with a significantly coarser mesh as compared to the resolution of the crack with collocation equations. In Section 4, we consider a simple one-dimensional model problem to assess the solution behavior of different isogeometric Galerkin and collocation variants, demonstrating the advantages of the hybrid collocation-Galerkin approach. In Section 5, we present numerical examples in two dimensions. We show that the computational cost for sufficiently fine meshes is still dominated by the collocation equations, since their number increases quadratically in 2D, but the number of local Galerkin equations increases only linearly. Section 6 summarizes our most important points and motivates future research. 2. Variational formulation of phase-field fracture models We start by briefly reviewing the variational formulation of the second-order and fourthorder phase-field fracture models based on an incremental variational principle that balances internal power with the power due to external loading. We follow the presentations given by Miehe, Welschinger and Hofacker [23, 26] and Borden et al. [21, 24]. In the scope of this work we focus on quasi-static fracture problems, assuming that loads are applied infinitely slowly. We therefore neglect inertia effects and use a pseudo-time t, which indexes each deformation state during the loading process. 2.1. Variational description of fracture We let Ω ⊂ Rd (with d ∈ {1; 2; 3}) be an arbitrary body with external boundary ∂Ω and internal displacement discontinuity boundary Γ. We assume small strain isotropic elasticity and recall the following key equations. The symmetric strain tensor at each material point is given as  1 ε= (1) ∇u + ∇uT 2 where u(x, t) denotes the displacement at each material point x and time instant t. The elastic energy stored in the bulk of the solid per unit volume is described at each material point by the energy density function 1 ψ0 (ε) = λ tr (ε)2 + µ ε : ε 2

(2)

where λ and µ are the Lam´e constants [52]. The stresses at each material point can then be obtained as σ := ∂ε ψ0 (ε) = λ tr(ε) I + 2µ ε (3) The evolving internal discontinuity boundary Γ(t) represents a set of discrete cracks, often simply referred to as the fracture topology. Following the classical treatments of 5

∂Ωg

∂Ωg

c( x , t)



c=1



Γ

l0 x2

x2

c=0 x3

x1

∂Ωh

x3

(a) A solid body Ω with internal discontiuity boundary Γ that represents a crack.

x1

∂Ωh

(b) Phase-field approximation of Γ. Parameter l0 controls the diffusive fracture width.

Figure 1: The concept of phase-field fracture.

Griffith [53] and Irwin [54], the work required to create a unit area of fracture surface is equal to the critical energy release rate Gc . The total internal potential of the body is then Z Z Ψ(ε, Γ) = ψ0 (ε) dΩ + Gc dΓ (4) Ω

Γ

The first term represents the energy stored in the bulk of the elastic solid, and the second term represents the work necessary to create the current fracture topology Γ(t). The variational approach to fracture proposed by Francfort and Marigo [2] predicts the nucleation, propagation and interaction of cracks by finding a global minimizer of (4) for a given load. 2.2. Phase-field approximation of the fracture surface Solving this variational problem numerically for discrete cracks can be difficult because of the evolution of the crack path Γ in time. In order to circumvent this difficulty, Bourdin et al. [1] introduced a volumetric approximation to the surface integral Z Z Gc dΓ ≈ Gc Γc,n dΩ (5) This approximation uses a smooth scalar-valued phase-field, c ∈ [0; 1], to represent the crack. This phase-field has a value of one away from the crack and a value of zero at the crack (see Fig. 1). The phase-field approximation introduces a crack density function Γc,n that depends on a length-scale parameter l0 and the phase-field c with derivatives up to order n. In the scope of this work, we consider the following two crack density functions. The function originally proposed by Bourdin et al. [1] reads Γc,2 =

 1  (c − 1)2 + 4l02 |∇c|2 4l0

(6)

and has been utilized in several contributions since then (see e.g. [21, 23, 26]). In a recent 6

1

1 Γ

c(x)

Γ

x

−2ℓ0 2ℓ0 (a) Second-order phase-field theory.

c(x)

−2ℓ0 2ℓ0

x

(b) Fourth-order phase-field theory.

Figure 2: One-dimensional phase-field approximation to the crack surface Γ = {0} (see [24]).

paper, Borden et al. [24] presented a higher-order functional, which reads Γc,4 =

 1  (c − 1)2 + 2l02 |∇c|2 + l04 (∆c)2 4l0

(7)

Since the strong form problems corresponding to (6) and (7) feature second-order derivatives and fourth-order derivatives of c, respectively, they are referred to as the second-order and fourth-order phase-field fracture theories [24]. We note that the definition (6) is well-posed variationally for all c ∈ H 1 (Ω), but its solutions will have in general no greater regularity. The definition (7) is well-posed for all c ∈ H 2 (Ω). Figure 2 illustrates the difference in regularity in the second-order and fourthorder phase-field models by plotting the exact solution to a crack surface Γ = {0} in a one-dimensional bar of infinite length. Borden et al. [24] showed that the greater regularity in the exact solution of the fourth-order theory provides better accuracy and convergence rates for the corresponding numerical solutions. 2.3. Additive split of the strain tensor Tensile stresses cannot be carried across a crack (crack opening), but compressive stresses can (crack closing). In order to be able to account for stress degradation in tension only, we require a split of the strain tensor into a positive and a negative part ε = ε+ + ε−

(8)

where the former describes the tensile mode and the latter the compressive mode contained in ε. The split is defined based on the spectral decomposition +

ε = ε− =

d X

i=1 d X

hεi i+ ni ⊗ ni

(9)

hεi i− ni ⊗ ni

(10)

i=1

7

where εi (with i = 1, . . . , d) are principal strains and ni the corresponding principal directions of the strain tensor. The bracket operators h·i+ and h·i− describe ramp functions and are defined (for an arbitrary input value x) as follows  0 if x < 0 + hxi = (11) x if x ≥ 0  x if x < 0 − hxi = (12) 0 if x ≥ 0 The derivatives of (9) and (10) with respect to the total strains   P+ := ∂ε ε+ (ε)   P− := ∂ε ε− (ε) = I − P+

(13) (14)

define two projection operators, which are isotropic tensor-valued functions of ε. These fourth-order tensors project the total strains onto its positive and negative parts, i.e. ε+ = P+ : ε and ε− = P− : ε. Algorithms to compute these objects are described by De Souza Neto et al. [55] (see Appendix A.5) and Miehe [56]. Similar projection tensors are very common in damage mechanics [57–59]. 2.4. Energy degradation using a phase-field We introduce the phase-field as a multiplication factor to tensile components of the reference energy density function (2) associated with the undamaged elastic solid. From a physical point of view, we use the phase-field to locally penalize the capability of the material to carry tension across cracks. This constitutive assumption constitutes the key component of phase-field fracture modeling and can be expressed in the form   ˜ c) = (1 − κ)c2 + κ ψ + (ε) + ψ − (ε) ψ(ε, (15) 0 0

The reference energy density function ψ0 is additively decomposed into a positive part due to tension and a negative part due to compression 2 1 ψ0+ (ε) = λ htr (ε)i+ + µ ε+ : ε+ 2 2 1 ψ0− (ε) = λ htr (ε)i− + µ ε− : ε− 2

(16) (17)

using the additive split of the strain tensor (9) and (10) and the bracket operators (11) and (12). The model parameter κ ≪ 1 introduced by Ambrosio and Tortorelli [60] prevents the full degradation of the stored energy by maintaining a small artificial tensile energy density κψ0+ (ε) at the fully-broken state c = 0. Although Braides [10] from a theoretical point of view and Borden et al. [24] from a computational point of view found κ to be unnecessary, we maintain this parameter as a way to control the algebraic conditioning number of the discrete system of equations in the case of partly or fully broken systems. 8

2.5. Variational formulation The variational form of the phase-field fracture problem can be derived from an incremental variational principle that balances the rate of different energy terms. In the quasi-static case considered here, we have the rate of the stored energy, the rate of dissipated energy due to the work done by fracture, and the rate of the energy due to the work of external forces. Their balance leads to the following rate of energy functional E˙int + F˙ frac,n − P˙ ext = 0

(18)

The rate of the stored energy can be computed from the stored energy function (15) as Z Z Z  d + − ˜ ˙ Eint = σ + σ : ε˙ dΩ + 2(1 − κ)ψ0+ c c˙ dΩ (19) ψ dΩ = dt where the tensile and compressive parts of the stress tensor follow from (3) and (15) as   σ + := (1 − κ)c2 + κ ∂ε ψ(ε)+ 0    2 (20) = (1 − κ)c + κ λ htr(ε)i+ I + 2µ ε+ − − σ − := ∂ε ψ(ε)− 0 = λ htr(ε)i I + 2µ ε

(21)

The second term of (19) contains the positive part of the reference energy density function ψ0+ and describes the local intensity of the tensile part of the deformation. This term can be interpreted as an energetic force that drives the crack evolution [23]. The rate of the dissipated energy due to work done by fracture follows with the phase-field approximation of the fracture topology (5) as Z d ˙ Ffrac,n = Gc Γc,n dΩ (22) dt and can be computed Z ˙ Ffrac,2 = Z ˙ Ffrac,4 =

for the second-order and fourth-order phase-field models as Z Gc ˙ dΩ (c − 1) c˙ dΩ + 2Gc l0 ∇c (∇c) 2l0 Z Z Gc Gc l02 ˙ ˙ dΩ (c − 1) c˙ dΩ + 4Gc l0 ∇c (∇c) dΩ + ∆c (∆c) 2l0 2

The rate of energy due to the work done by external forces is simply Z Z Pext = b · u˙ dΩ + tˆ · u˙ d∂Ω

(23) (24)

(25)

where b and tˆ denote body forces and boundary tractions, respectively. We can bring all terms of the rate of energy functional (18) in a form where they depend only on the rate of the displacements u˙ or the rate of the phase-field c. ˙ With the argument 9

that the balance (18) must hold for arbitrary u˙ and c˙ we can separate the rate of energy functional into a phase-field part and an elasticity part. The phase-field part reads  Z  Z Z 4l0 (1 − κ)ψ0+ 2 + 1 c c˙ dΩ + 4l0 ∇c ∇c˙ dΩ = c˙ dΩ (26) Gc  Z  Z Z Z 4l0 (1 − κ)ψ0+ 2 4 + 1 c c˙ dΩ + 2l0 ∇c ∇c˙ dΩ + l0 ∆c ∆c˙ dΩ = c˙ dΩ (27) Gc for the second-order and fourth-order model, respectively. The elasticity part reads Z Z Z  + − σ + σ : ∇u˙ dΩ = b · u˙ dΩ + t · u˙ d∂Ω

(28)

They constitute the weak form of the phase-field fracture problem, where we identify u˙ and c˙ as test functions. In the following, we will refer to these test functions as w and q. Applying integration by parts to (26) or (27) and (28) we can derive the strong form of the corresponding phase-field fracture problem by identifying a system of coupled partial differential equations and boundary conditions [24]. We note that the strong and weak forms of the phase-field fracture problem can be alternatively derived using calculus of variations, where we search for stationary points of the rate of energy functional (18) [22, 23, 26]. 2.6. Irreversibility and history-based formulation The evolution of the crack topology is driven by the positive part of the reference energy density function ψ0+ that emanates from (19). This can be best illustrated by the strong form of the phase-field equation that reads for the second-order model   4l0 (1 − κ)ψ0+ 1 − c + 4l0 ∆c = c (29) Gc Considering the limit cases, we can easily observe that the phase-field solution c is one if ψ0+ is zero, and goes to zero if ψ0+ becomes very large. However, there is nothing in (29) that prevents cracks from healing if loads are removed and ψ0+ decreases. An intuitive, simple and effective idea to incorporate the irreversibility of the crack topology is the introduction of a local history field H(x, t) = max ψ0+ (ε(x, s)) s∈[0;t]

(30)

that records the maximum positive reference energy density that has occurred during loading at a specific location x up to the current time instant t. Replacing ψ0+ by H in (29) prevents cracks from healing once they have opened up, since the corresponding value of ψ0+ is maintained and cannot be reversed when the solid is unloaded. The constitutive assumption (30) accommodates the irreversibility condition F˙ frac,2 (c, c) ˙ ≥0 10

(31)

preventing a decrease of the work done by fracture. Enforcing (31) guarantees that the corresponding energy is dissipated from the system and cannot be returned. The history field based procedure directly transfers to the fourth-order phase-field model [24]. 3. Numerical treatment with isogeometric collocation In this section, we extend the isogeometric hybrid collocation-Galerkin method outlined by De Lorenzis et al. [45] for the discretization of phase-field fracture models. In particular, we will show that the hybrid approach offers a natural and consistent way of imposing weak Neumann boundary conditions and weak multi-patch interface constraints and is able to handle the multiple boundary integral terms that appear in the weighted residual form in the fourth-order phase-field model. We assume that the reader has basic knowledge of B-spline and NURBS basis functions. For an introduction in the context of isogeoemtric analysis, we refer to Hughes et al. [36] and Cottrell et al. [35]. More details on different aspects can be found e.g. in [61–64]. We note that when we talk about NURBS in the following, we also include B-splines, as every B-spline is a special case of a NURBS. 3.1. The weighted residual form The variational background for the derivation of isogeometric collocation methods is the weighted residual form [65, 66]. We can derive weighted residual statements for the phasefield and elasticity parts from the variational forms (26), (27) and (28) using integration by parts, shifting the derivative operators from the test functions w and q back onto the the primary variables u and c. We note that the integration by parts procedure generates surface flux terms that introduce the Neumann boundary conditions explicitly. Applying integration by parts to the second term in (26), we find the weighted residual form of the phase-field part for the second-order phase-field theory as  Z  Z Z Z 4l0 H 2 2 + 1 c q dΩ − 4l0 ∆c q dΩ + 4l0 (∇c · n) q d∂Ω = q dΩ (32) Gc Ω Ω ∂Ω Ω where n denotes the outward unit normal vector. Applying integration by parts to the second term and integration by parts twice to the third term in (27), we find the weighted residual form of the phase-field part for the fourth-order phase-field theory as  Z  Z Z 4l0 H 2 + 1 c q dΩ − 2l0 ∆c q dΩ + l04 ∆(∆c) q dΩ + G c Ω Ω Z Z ZΩ   2 4 4 ∇ 2l0 c − l0 ∆c · n q d∂Ω + l0 ∆c n · ∇q d∂Ω = q dΩ (33) ∂Ω

∂Ω



Integrating the first term in (28) by parts, we find the weighted residual form of elasticity Z Z    +    + − (34) σ + σ − n − tˆ · w d∂Ω = 0 div σ + σ + b · w dΩ + − Ω

∂Ω

11

1 Collocation points

0.5

0 −1

−0.875

−0.75

−0.625

−0.5

−0.375

−0.25

−0.125

0

0.125

0.25

0.375

0.5

0.625

0.75

0.875

1

Figure 3: 1D patch of quadratic NURBS basis functions and corresponding Greville points (red dots).

Before discretization the variational statements (26) to (28) and (32) to (34) are identical. However, their discretization requires basis functions with different minimum continuity requirements, since the variational index, that is, the maximum order of derivatives, changes as a consequence of the integration by parts procedure. In the cases considered here the weighted residual form requires the following function spaces on each patch: Sc,2 = {c ∈ C 1 (Ω)} Vc,2 = {q ∈ L1loc (Ω)}

(35)

for the phase-field part for the second-order theory, Sc,4 = {c ∈ C 3 (Ω)} Vc,4 = {q ∈ L1loc | q ∈ C 0 (∂Ω)}

(36)

for the phase-field part for the fourth-order theory, and ¯ Su = {u ∈ C 1 (Ω) | u(∂Ωg ) = u} 1 Vu = {w ∈ Lloc | w(∂Ωg ) = 0}

(37)

for the elasticity part. C m (m≥0) denotes the class of continuous and differentiable functions whose derivatives up to order m are continuous, and L1loc denotes the class of locally integrable functions without any continuity requirement. ∂Ωg denotes the Dirichlet part of ¯ are prescribed. the domain boundary where displacements u 3.2. Hybrid collocation-Galerkin discretization in space We consider approximations u∗ and c∗ to the exact displacements u and phase-field c ∗

u (x, t) = c∗ (x, t) =

n X

i=1 n X

Ni (x) ui (t)

(38)

Ni (x) ci (t)

(39)

i=1

The displacements and the phase-field are discretized by the same set of spline basis functions. The basis Ni (x) (i = 1 . . . n) consist of a set of n NURBS functions that satisfy 12

corresponding smoothness requirements (35) or (36), and (37). Each basis function is multiplied by coefficients ui (t) and ci (t), where we neglect the time dependence until Section 3.5. Each vector ui contains d scalar coefficients for the displacement field of each spatial direction, where d denotes the dimension. NURBS are defined over patches, where higher-order smoothness is a global property. Each patch is constructed with open knot vectors, so that boundary functions are interpolatory. Figure 3 shows an example patch of quadratic NURBS basis functions in one dimension. Multiple patches can be easily connected using their interpolatory property at the patch interfaces [35, 36]. For the discretization of the test functions w and q we consider two different sets of functions. The first set contains the same NURBS basis functions, Ni (x), used in the approximation of the solution fields (38) and (39). The second set consists of Dirac δ functions, constructed as the limit of a sequence of smooth functions with compact support that converge to a distribution [37, 42]. Depending on which set we choose, we obtain the standard Galerkin method or a point-collocation scheme [37, 38, 67, 68]. 3.2.1. Collocation at the Greville points Point-collocation methods make use of the sifting property of Dirac δ functions. This property can be summarized as Z g(x) δ(x − x ˆi ) dΩ = g(ˆ xi ) (40) Ω

provided that g is a continuous function about the point x ˆi . The Dirac δ test functions are defined at specific points with coordinates x ˆi , the so-called collocation points. In isogeometric collocation based on NURBS, we choose the Greville points as collocation points [38, 42]. For a one-dimensional NURBS patch of polynomial degree p, Greville points can be constructed with the Greville abscissa that reads τˆiξ

i+p−1 1 X = ξj p j=i

(41)

where ξj denote consecutive entries of the corresponding knot vector. Greville points for multi-variate NURBS discretizations can be constructed by taking the tensor product of one-dimensional Greville points, set up individually for each parametric direction. The advantages of the Greville points in the context of isogeometric collocation are twofold: First their number exactly corresponds to the number of NURBS basis functions, since both are constructed from the same knot vector structure. Hence a one-to-one correspondence can be established between a specific point and a specific NURBS function. Second Greville points have been shown to lead to stable and robust collocation schemes in previous studies [37, 38, 42, 46]. Figure 3 also illustrates the Greville points for the one-dimensional example patch.

13

3.2.2. Local switch between Galerkin and collocation For each NURBS basis function used in the approximation of the solution fields (38) and (39), we can make a choice: we can either use the function itself to weight the residual in the Galerkin sense or we can collocate at the corresponding collocation point. Based on our specific choice the discretization w∗ and q ∗ of the exact test functions w and q read ∗

w (x) = ∗

q (x) =

k X

Ni (x) wi +

n X

δ(x − x ˆ i ) wi

i=1

i=k+1

k X

n X

Ni (x) qi +

i=1

δ(x − x ˆ i ) qi

(42) (43)

i=k+1

where we gather the k NURBS functions in the first term and the n-k Dirac δ functions at the Greville points in the second term. Each vector wi contains again a scalar coefficient for each spatial direction. 3.2.3. Weak enforcement of Neumann boundary conditions The original concept of isogeometric collocation envisaged the use of Dirac δ test functions at the patch boundary in the sense of a basic point-collocation method [37, 38, 42]. However, De Lorenzis et al. [45] discovered that the point-wise exact enforcement of Neumann boundary conditions at boundary collocation points can lead to spurious oscillations. The flexibility of the local switch between Galerkin and collocation based on (42) and (43) allows us to use the outer 1-layer of NURBS basis functions at patch boundaries as test functions. The Galerkin method enforces a weighted average between the domain and boundary operators, where the weighting factors naturally occur as a result of the integration of the volume and boundary terms [69]. De Lorenzis et al. [45] showed that this approach effectively prevents oscillations. The same idea has also been successfully used in the context of spectral collocation [67, 70] and collocated C 0 finite elements [69]. Figure 4a illustrates Galerkin and collocated NURBS basis functions for the one-dimensional example patch. 3.2.4. Handling more complex Neumann boundary conditions A further advantage of the local switch between Galerkin and collocation is the possibility to consistently handle more complex Neumann boundary conditions as they appear in the fourth-order phase-field theory. Looking at the corresponding weighted residual form (33), we observe that there are two Neumann boundary condition terms, one of which contains first-order derivatives on the phase-field test function. Note that first-order derivatives in boundary terms cannot be removed by an integration by parts procedure. The hybrid collocation-Galerkin method can naturally deal with this situation by simply using more NURBS basis functions at the boundary in the test function discretization (43). The rule is as follows: If there is no derivative on the test function along the boundary, only the outer 1-layer of NURBS basis functions is required. If there are derivatives on the test function along the boundary, the number of outer layers of NURBS functions needs to be increased, so that all NURBS basis functions with derivatives of that order appear in the 14

1 Collocation point

0.5

0 −1

−0.875

−0.75

−0.625

−0.5

−0.375

−0.25

−0.125

0

0.125

0.25

0.375

0.5

0.625

0.75

0.875

1

(a) Second-order phase-field problem (quadratic NURBS, 1-layer of Galerkin test functions). 1 Collocation points

0.5

0 −1

−0.875

−0.75

−0.625

−0.5

−0.375

−0.25

−0.125

0

0.125

0.25

0.375

0.5

0.625

0.75

0.875

1

(b) Fourth-order phase-field problem (quartic NURBS, 2-layer of Galerkin test functions). Figure 4: Example patch in 1D: NURBS basis functions used in a Galerkin sense (dashed green lines), collocated basis functions (solid blue lines), and Greville collocation points (red dots). Note that all B´ezier elements supporting at least one NURBS test function require Gauss quadrature.

test function discretization (43) in a Galerkin sense. For the fourth-order phase-field model, this means we use the outer 2-layer of NURBS functions at the patch boundary, which is illustrated in Fig. 4b for the 1D example patch. 3.3. Discrete forms We describe the stiffness matrices and force vectors that arise when we discretize the weighted residual forms of the phase-field fracture models using (38), (39), (42), and (43). 3.3.1. Phase-field part of the second-order theory For the second-order theory, substitution of the phase-field approximation (39) and the test function (43) into the corresponding weighted residual form (32) yields k X i=1

k X

i=1 n X i=k

Z  n h X i  qi Dc,2 Nj (x) cj − 1 Ni (x) dΩ + Ω

qi

Z

j=1

Bc,2 ∂Ω

n h X

Nj (x) cj

j=1

i

Ni (x) d∂Ω +

n  h X i  qi Dc,2 Nj (ˆ xi ) cj − 1 = 0 j=1

15

(44)

where we use a compact notation based on the domain and boundary operators    ∗ 4l0 H Dc,2 c = + 1 c∗ − 4 l02 ∆c∗ Gc   Bc,2 c∗ = 4 l02 ∇c∗ · n

(45) (46)

We observe that for all collocated basis functions the integrals are naturally eliminated due to the sifting property (40) of the Dirac δ functions. We assume in (44) that the value of the history field H(u∗ ) is known at each location x. Since coefficients qi are arbitrary, equation (44) yields a system of linear algebraic equations with unknowns cj . Collecting terms and factoring out cj we obtain the stiffness matrix  Z hP i   Dc,2 Nj (x) Ni (x) dΩ +     Z for 1 ≤ i ≤ k h i  P c,2 Bc,2 Nj (x) Ni (x) d∂Ω Kij = (47)    hP i     Dc,2 Nj (ˆ xi ) for k + 1 ≤ i ≤ n

with k NURBS basis functions that have support at the patch boundary and are used in the discretization of the test function in a Galerkin sense, and n-k collocated NURBS basis functions in the interior. The corresponding force vector reads  Z  Ni (x) dΩ for 1 ≤ i ≤ k (48) Fic,2 =  1 for k + 1 ≤ i ≤ n 3.3.2. Phase-field part of the fourth-order theory When we use the fourth-order theory, substitution of the phase-field approximation (39) and the test function (43) into the weighted residual form (33) yields k X i=1

k X

Z  n h X i  qi Dc,4 Nj (x) cj − 1 Ni (x) dΩ + Ω

qi

i=1

k X

i=1 n X i=k

qi

Z Z

∂Ω

∂Ω

j=1

1 Bc,4

n h X

Nj (x) cj

i

Ni (x) d∂Ω +

n h X

Nj (x) cj

i

· ∇Ni (x) d∂Ω +

j=1

2 Bc,4

j=1

n  h X i  qi Dc,4 Nj (ˆ xi ) cj − 1 = 0 j=1

16

(49)

with the following domain and boundary operators    ∗ 4l0 H Dc,4 c = + 1 c∗ − 2l02 ∆c∗ + l04 ∆(∆c∗ ) Gc    ∗ 1 Bc,4 c = ∇ 2l02 c∗ − l04 ∆c∗ · n  ∗ 2 Bc,4 c = l04 ∆c∗ n

(50) (51) (52)

We assume again that the value of the history field H(u∗ (x)) is known at each location x. We obtain the following stiffness matrix  Z hP i  D N (x) Ni (x) dΩ +  c,4 j    Z  hP i   1  Nj (x) Ni (x) d∂Ω for 1 ≤ i ≤ k Bc,4  c,4 Z Kij = (53) hP i  2  B N (x) · ∇N (x) d∂Ω  j i c,4     hP i    D Nj (ˆ xi ) for k + 1 ≤ i ≤ n c,4

with k NURBS basis functions whose first derivatives have support at the patch boundary and n-k collocated NURBS basis functions in the interior. We note that the force vector for the fourth-order theory is the same as the force vector (48) for the second-order theory. 3.3.3. Elasticity part The discrete forms of the elasticity part are based on the weighted residual form (34) and have the same form for both the second-order and fourth-order phase-field theories. Substitution of the displacement approximation (38) and the test function (42) yields −

k X

wi

i=1 n X i=k

Du



i=1

k X

Z 

wi

Z

∂Ω



n h X

Nj (x) uj

j=1

Bu

n h X

i

Nj (x) uj

j=1

 + b Ni (x) dΩ + i

 ˆ − t Ni (x) d∂Ω +

n  h X i  wi D u Nj (ˆ xi ) uj + b(ˆ xi ) = 0

(54)

j=1

where we use a compact notation based on the domain and boundary operators     Du u∗ = div σ + (u∗ ) + σ − (u∗ )     Bu u∗ = σ + (u∗ ) + σ − (u∗ ) n

(55) (56)

For all collocated basis functions the integrals are naturally eliminated due to the sifting property (40) of the Dirac δ test functions. We assume that in (54) the value of the phase17

field c∗ of the current fracture state is known at each location x. Since coefficients wi are arbitrary, equation (54) yields a system of linear algebraic equations with unknowns uj . Collecting terms and factoring out uj we obtain the stiffness matrix of the elasticity part  Z hP i   Du Nj (x) Ni (x) dΩ +     Z for 1 ≤ i ≤ k h i  P u Bu Nj (x) Ni (x) d∂Ω Kij = (57)    hP i     Du Nj (ˆ xi ) for k + 1 ≤ i ≤ n

and the corresponding force vector  Z Z  b · Ni (x) dΩ + tˆ · Ni (x) d∂Ω for 1 ≤ i ≤ k Fiu =  b(ˆ x) for k + 1 ≤ i ≤ n

(58)

i

with k NURBS basis functions that have support at the patch boundary and n-k collocated NURBS basis functions in the interior. 3.4. Weak enforcement of interface constraints in multi-patch discretizations NURBS discretizations usually consist of multiple NURBS patches that are connected along patch interfaces. At these interfaces, NURBS basis functions are still conforming, but their continuity degrades to C 0 [35, 36]. As a consequence, weighted residual forms cannot be derived by using integration by parts over the complete domain, since the higher-order continuity requirements (35), (36), and (37) for the approximation spaces S are not met at the patch interfaces. The hybrid collocation-Galerkin approach suggests a treatment of multi-patch discretizations along the lines of collocated finite element methods that use standard nodal C 0 basis functions [67, 69, 70]. In these methods, the weighted residual form is derived individually for each element by using the smoothness over each element domain. The compatibility of the solution fields across different elements is enforced weakly by the additional flux terms that arise due to integration by parts at the element boundaries. They ensure that at each element interface the residual of the original differential equation and the element boundary flux terms from all neighboring elements are in equilibrium in a weighted average sense. The corresponding discrete forms can be simply established by evaluation of the weighted residual form for each element with subsequent assembly. This procedure holds equivalently for multi-patch discretizations with C 0 patch interfaces, where we treat each patch as a collocated finite element. We simply derive the weighted residual form (32), (33), and (34) for each patch, with the only difference that there is no external traction tˆ along patch interfaces. We then evaluate the discrete forms for each patch separately. Interface conditions that tie adjacent patches together are established by the flux terms that occur at patch interfaces. They are automatically taken 18

The displacements un the phase-field cn , and the history field Hn are known at the current time tn . Update prescribed body forces bn+1 and boundary tractions tˆn+1 to the new time tn+1 . 1. Update phase-field (phase-field part):     c,2 Kij cj = Fic,2 second-order is to be employed if the Solve linear system fourth-order Kijc,4 cj = Fic,4 using (47), (48), and (53) at frozen displacements u. 2. Update displacements (elasticity part): Solve linear system Kiju uj = Fiu using (57) and (58) at frozen phase-field c. 3. Update history field: Determine maximum density during loading history  + reference energy + ψ0 (un+1 ) if ψ0 (un+1 ) > Hn Hn+1 (un+1 ) = Hn otherwise Update time variable tn+1 to tn and proceed to next step by restarting this procedure. Table 1: One time step [tn , tn+1 ] of the staggered solution algorithm for phase-field fracture [23, 26]

care of by evaluating the boundary terms in each patch and assembling the results into the global discrete forms. In particular, due to the conforming NURBS basis functions along patch interfaces, there is no need to exchange additional information across patches during formation and assembly of the discrete forms. In analogy to the single-patch case hybrid collocation-Galerkin uses outer layers of NURBS functions along each patch boundary for the discretization of test functions.This ensures that the weighting between domain and boundary terms along patch interfaces is consistent with the integrals of the weighted residual form, which has been shown to be important for accuracy and stability in collocated finite element methods [67, 69, 70]. 3.5. Staggered solution algorithm For the discretization in pseudo-time t, we apply the staggered solution algorithm based on operator splits that solves the elasticity and phase-field parts independently. For phasefield fracture, it was suggested by Miehe et al. [23] and used successfully in several other studies [21, 24, 26, 27]. In our case, we solve the phase-field part first to update the phasefield variable and then solve the elasticity part that updates the displacements. Finally we update the history field that records the maximum reference energy density. The algorithm is summarized in Table 1. 19

sin(πx) x -1.0

0.0

1.0

Figure 5: One-dimensional model problem of a bar with a centered crack with a sinusoidal load.

4. Analysis in one dimension In this section, we will illustrate the accuracy of isogeometric collocation and compare it to the standard Galerkin method. In particular, we will show that we can use the hybrid collocation-Galerkin approach to adaptively resolve the fracture zone with Galerkin test functions. Using the same mesh size, Galerkin equations allow for the accurate resolution of phase-field approximations with significantly smaller length-scale parameters as compared to collocation equations based on Dirac δ test functions. In turn, for a given length-scale parameter, we can use a coarser mesh and still achieve an accurate crack resolution, if we adaptively replace Dirac δ by Galerkin test functions in the fracture zone. In the following examples, we will illustrate this advantage by using the length-scale as a numerical parameter that can be varied arbitrarily. However, the reader should always keep in mind that from a physical point of view the length-scale parameter and the associated diffusiveness are not arbitrary, but need to be perceived of as a specific material property that determines the critical stress at which crack nucleation occurs [21]. For a more general discussion on the role and importance of finite length-scale parameters in computational fracture mechanics, we refer e.g. to [71–73]. 4.1. One-dimensional model problem We consider the 1D model problem shown in Fig. 5, which consists of a bar that is fixed at both ends and loaded along its axis by a sinusoidal load. For this simple case, we assume the parameter κ to be equal to zero, Young’s modulus E and cross-section area A to be one, and the strain ε to be non-negative in the fracture zone. Hence, we can use the stress-strain relationship σ = c2 ε

(59)

without the additive split (9). We impose a crack at the center of the bar by an initial history field  1000.0 if d(x) ≤ l0 H0 (x) = (60) 0.0 if d(x) > l0 where d(x) denotes the distance to the center of the bar at x=0. The solution is obtained by one pass through the staggered algorithm of Table 1, where we solve first for the phase-field 20

−2

10

10 10

−2

−3

−4

Error in H semi−norm

2 1 −6

10

3 1 −8

4

10

−10

10

10 10

−4

2 1

−5

1

Error in H1 semi−norm

10

10

Quadratics (p=2) Cubics (p=3) Quartics (p=4) Quintics (p=5) 20

10 10

1 10

5 1

40 # degrees of freedom

80

10

160

−6

−7

Quadratics (p=2) Cubics (p=3) Quartics (p=4) Quintics (p=5)

−8

−9

10

20

(a) Galerkin.

4 1

40 # degrees of freedom

80

160

(b) Collocation.

Figure 6: Convergence of the error in stress (H1 semi-norm) for the undamaged bar problem.

and then for the displacements. We asses the quality of the solution by the error in stress e=

Z

2

σex − c ε

2

dx

1/2

(61)

which contains both the approximation of the phase-field c and approximation of the axial strain ε = ∂u/∂x. It can be interpreted as an extension of the error in the H1 semi-norm of the elasticity part to the phase-field fracture case. 4.2. Undamaged elastic bar First we illustrate the accuracy of isogeometric Galerkin and collocation methods for smooth problems. To this end, we consider the undamaged bar without the crack that we discretize with a uniform mesh of NURBS basis functions of polynomial degree p and continuity C p−1 . The corresponding exact stress solution is σex =

1 cos(πx) π

(62)

Figures 6a and 6b plot the convergence behavior of the isogeometric Galerkin and collocation methods under mesh refinement for the undamaged case. The Galerkin method achieves optimal rates of convergence with a rate of O(p) for the error in the H1 semi-norm [74, 75]. Isogeometric collocation converges with a rate of O(p), if the polynomial degree p is even, and with a rate of O(p − 1), if the polynomial degree p is odd. Based on numerical evidence, these convergence rates are the best possible rates for the error in the H 1 semi-norm that can be expected from an isogeometric collocation method [37, 38].

21

0

0

10

10

Galerkin Collocation Adaptive

Galerkin Collocation Adaptive −1

−1

10 Error in stress

Error in stress

10

−2

−3

10

0.5

−2

10

10

−3

1

1.5

2 2.5 3 3.5 Length ratio k (l0 = k h)

4

4.5

10

5

(a) Second-order phase-field model.

0.5

1

1.5

2 2.5 3 3.5 Length ratio k (l0 = k h)

4

4.5

5

(b) Fourth-order phase-field model.

Figure 7: Accuracy for varying length-scale l0 : We plot the error in stress vs. the ratio k = l0 /h obtained with 256 B´ezier elements (constant h). We use quadratic and quartic NURBS for the secondand fourth-order models, respectively, and assume the discontinuous solution as a reference.

4.3. Elastic bar with a crack We now consider the 1D model problem with a crack at the center (see Fig. 5). The discontinuous solution fields in the fully fractured case read  1 1+x   if x < 0  2 sin(πx) − π π uex = (63)    1 sin(πx) + 1 − x if x ≥ 0 π2 π σex =

1 1 cos(πx) − π π

(64)

From a physical point of view, we obtain two single bars that are not connected at x=0. Equations (63) and (64) represent the solutions that would be obtained from a phase-field model in the limit h → 0 and l0 → 0, where h denotes the characteristic length of the B´ezier elements and l0 is the length-scale parameter. We solve the test problem using the second- and fourth-order phase-field fracture models on a mesh of 256 B´ezier elements with quadratic and quartic NURBS, respectively. In particular, we are interested in the sensitivity of the accuracy of the isogeometric collocation and Galerkin methods with respect to the length-scale parameter l0 . To this end, we vary the length-scale parameter in the form l0 = k h

(65)

where k denotes the ratio between l0 and constant element length h. Figure 7a and 7b plot 22

0.4

0.6

0.3 0.2

0.2

Displacement u

Displacement u

0.4

Analytical Galerkin Collocation Adaptive

0

Analytical Galerkin Collocation Adaptive

0.1 0 −0.1

−0.2 −0.2 −0.4

−0.6 −1

−0.3

−0.5

0 0.5 Axial coordinate x of the bar

−0.4 −1

1

(a) Displacements.

0.4

Analytical Galerkin Collocation Adaptive

0.2 Stress c2 du/dx

2

Stress c du/dx

1

0.6 Analytical Galerkin Collocation Adaptive

0.2 0 −0.2

0 −0.2

−0.4

−0.4

−0.6

−0.6

−1

0 0.5 Axial coordinate x of the bar

(a) Displacements.

0.6 0.4

−0.5

−0.5

0 0.5 Axial coordinate x of the bar

1

−1

−0.5

0 0.5 Axial coordinate x of the bar

1

(b) Axial stress.

(b) Axial stress.

Figure 8: Solution fields obtained for the second-order phase-field model with 256 quadratic B´ezier elements and l0 = 2 h.

Figure 9: Solution fields obtained for the fourth-order phase-field model with 256 quartic B´ezier elements and l0 = 2 h.

the error in stress (61) with respect to different ratios k, where k varies from one half to five. We assume that the physical length-scale of the material is much smaller than l0 = 0.5 h, so that the discontinuous solution (64) can be used as a reference. We observe that the accuracy improves with decreasing ratio k, since the phase-field approximation approaches the reference as l0 gets smaller. However, for both methods, there exists a distinct critical ratio kc , below which the solution suddenly starts to diverge, although the length-scale parameter is further decreased. We note that this critical ratio is significantly smaller for the Galerkin method (kc ≈ 1.5) than for isogeometric collocation (kc ≈ 3.5). Figures 8 and 9 plot the displacement and stress solutions obtained with the Galerkin and collocation methods on a mesh of 256 B´ezier elements, using the second-order and 23

1

0.8

0.8 Phase field c

Phase field c

1

0.6

0.4

0.2

0 −0.15

−0.05 0 0.05 Axial coordinate x of the bar

0.1

0.4

0.2

Galerkin Collocation Adaptive −0.1

0.6

0 −0.15

0.15

Galerkin Collocation Adaptive −0.1

(a) Phase-field c.

0.1

0.15

(a) Phase-field c.

150

60 Galerkin Collocation Adaptive

100

Galerkin Collocation Adaptive 40 Strain du/dx

50 Strain du/dx

−0.05 0 0.05 Axial coordinate x of the bar

0 −50 −100

20

0

−150 −200 −0.15

−0.1

−0.05 0 0.05 Axial coordinate x of the bar

0.1

−20 −0.15

0.15

−0.1

−0.05 0 0.05 Axial coordinate x of the bar

0.1

0.15

(b) Axial strain ε.

(b) Axial strain ε.

Figure 10: Localization of c and ε in the fracture zone x ∈ [−0.15, 0.15] for the secondorder model (256 quadratic ele’s, l0 = 2 h).

Figure 11: Localization of c and ε in the fracture zone x ∈ [−0.15, 0.15] for the fourthorder model (256 quartic el’s, l0 = 2 h).

fourth-order theories, respectively. All computations use a length-scale of l0 = 2 h. We note that for isogeometric Galerkin, k is above the critical ratio, while for isogeometric collocation k is below the critical ratio (see Figs. 7a and 7b). We observe that the Galerkin method can accurately resolve the crack, yielding good results that do not exhibit any oscillations. Isogeometric collocation, however, leads to spurious oscillations in the solution fields around the crack. We can also see that this phenomenon has an influence on the global accuracy of the solution fields away from the fracture zone. We conclude that there exists a lower bound for the length-scale parameter l0 that each numerical scheme is able to resolve, which can be expressed as a critical ratio kc with respect to the mesh size. In particular, this lower bound is significantly smaller for Galerkin than for collocation. 24

1

0.5

0 −1

−0.875

−0.75

−0.625

−0.5

−0.375

−0.25

−0.125

0

0.125

0.25

0.375

0.5

0.625

0.75

0.875

1

Figure 12: Adaptive resolution of the crack at x = 0 by a small band of NURBS test functions (dashed green lines) that are used in a Galerkin sense. Away from fracture zone, we use collocated NURBS basis functions (blue solid lines) and Greville collocation points (red dots).

Figure 10 and 11 offer further insight by separately plotting the phase-field and strain solutions for the second- and fourth-order models, respectively. The multiplication of the two fields yields the stress solution according to (59). Both plots zoom into the fracture zone. We observe that the phase-field solutions obtained with the Galerkin and collocation methods are virtually lying on top of each other. This indicates that both methods yield an equivalent accuracy in approximating the phase-field. All strain fields are localized in the fracture zone. However, the strains obtained with the collocation method exhibit a broader band and distinct oscillations. This indicates that the difference in strain localization in the fracture zone is the reason for the different performance of Galerkin and collocation. The phase-field is able to successfully mitigate the influence the corresponding stresses in the Galerkin method. In collocation, it fails to sufficiently mitigate the corresponding stresses in the transition region of the diffusive fracture zone due to a larger strain band. 4.4. Adaptive crack resolution using the hybrid collocation-Galerkin approach The hybrid collocation-Galerkin method offers an opportunity to remedy the insufficient strain localization of collocation at the crack. We can use its flexibility to locally resolve the fracture zone with Galerkin test functions, while collocating further away from the crack. For the 1D model problem at hand, we increase the approximation power of isogeometric collocation by using a small band of NURBS basis functions that are closest to x=0 as Galerkin test functions. This idea is illustrated in Fig. 12. Figures 7a and 7b illustrate the accuracy of the adaptive collocation-Galerkin method for different ratios k, comparing the error in stress to the full Galerkin method and the collocation method that uses Dirac δ test functions in the fracture zone. We observe that the adaptive resolution of the crack with a few Galerkin test functions reduces the critical ratio kc and enables isogeometric collocation to use the same small length-scale parameter l0 as the Galerkin method. Figures 8 through 11 illustrate that the corresponding solution fields of the adaptive approach are accurate and free of oscillations and its strain localization at the crack is equivalent to the Galerkin method. At the same time, the adaptive collocation-Galerkin method has practically the same computational efficiency as a basic collocation method, when sufficiently fine meshes are considered. For the one-dimensional problem at hand, the number of collocated basis functions increases linearly under mesh refinement, while the number of B´ezier elements that need to be fully integrated due to the presence of Galerkin 25

u

0.5mm

0.5mm

0.5mm

0.5mm

Figure 13: Geometry and boundary conditions for the tension test of a single edge notched specimen.

test functions does not change. The adaptive approach thus combines the advantages of the Galerkin method in terms of the ability to guarantee strain localization at the crack, which results in accurate stresses and displacements, and the advantages of collocation in terms of a significant reduction of point evaluations. 5. Numerical examples in two dimensions In this section we examine benchmark problems in two dimensions that demonstrate the advantages of isogeometric collocation for the solution of phase-field fracture models. In particular, we show that collocation with an adaptive Galerkin resolution of the fracture zone enables us to accurately solve phase-field models with significantly smaller length-scale parameters than isogeometric collocation that resolves the crack with Dirac δ test functions. In turn, for a given length-scale parameter, we can use a coarser mesh to achieve the same accuracy. At the same time, the adaptive collocation-Galerkin approach significantly reduces the number of point evaluations in the bulk of the domain away from the fracture zone as compared to a full Galerkin approach. 5.1. Single edge notched tension test We first consider the square plate in plain strain with a horizontal notch described by Miehe et al. [22, 23]. The geometry and boundary conditions of this problem are illustrated in Fig. 13, and we use Young’s modulus E=210 GPa, Poisson’s ratio ν=0.3, and a critical energy release rate of Gc =2700 J/m2 . We discretize the square plate with one patch of C p−1 continuous NURBS basis functions and model the notch by an induced crack in the initial phase-field based on equation (60). The displacement at the top is imposed in increments of ∆u=1×10−5 mm. We use quadratic NURBS for the second-order phase-field model and quartic NURBS for the fourth-order phase-field model. We employ the isogeometric Galerkin method and the two isogeometric collocation methods without and with adaptive Galerkin resolution of the fracture zone. 26

(a)

(b)

(c)

Figure 14: 50×50 quadratic B´ezier elements: (a) Quadrature points over the complete domain. (b) Quadrature (green) only at the boundary, collocation (red) in the interior. (c) Quadrature at the boundary and in the fracture zone, collocation elsewhere.

(a)

(b)

(c)

Figure 15: 250×250 quadratic B´ezier elements: (a) Quadrature points over the complete domain. (b) Quadrature (green) only at the boundary, collocation (red) in the interior. (c) Quadrature at the boundary and in the fracture zone, collocation elsewhere.

5.1.1. Cost estimate based on the number of point evaluations We first estimate the cost for the formation and assembly of the discrete forms that are required for each method. In the scope of this work we correlate the formation and assembly cost of a method with the number of point evaluations. We note that this is a conservative assumption, since taking into account the cost per point evaluation tends to increase the advantage of collocation over Galerkin (see e.g. Schillinger et al. [38]). Figures 14 and 15 give a visual impression of the quadrature and collocation points that need to be evaluated in each method for two different mesh sizes. Numerical quadrature is required in each B´ezier element that supports at least one Galerkin test function. We use reduced quadrature based on p × p Gauss points, which has been shown to guarantee full accuracy and stability in isogeometric analysis [41, 76–78]. Collocated basis functions require the evaluation of only 27

10

10

8

9

10 Galerkin Point collocation only

7

Galerkin Point collocation only

8

10

10

10

10

6

5

# point evaluations

# point evaluations

7

10

Collocation w Galerkin at boundary and in fracture zone

4

Collocation w Galerkin at boundary and in fracture zone

6

10

5

10

4

Collocation w Galerkin in 2−layer at boundary only

10 10

Collocation w Galerkin in 1−layer at boundary only

3

3

10

2

2

10 2 10

3

10

4

5

10 10 # Bézier elements

6

10

10

10 2 10

7

(a) Quadratics (second-order phase-field).

3

10

4

5

10 10 # Bézier elements

6

10

10

7

(b) Quartics (fourth-order phase-field).

Figure 16: Square plate: number of point evaluations with increasing mesh size for different methods.

one collocation point generated from tensor-product Greville abscissae. We plot quadrature points in green and collocation points in red. The Galerkin method requires the use of quadrature points over the complete domain (see Figures 14a and 15a). In the finer mesh, single points cannot be distinguished anymore. Assuming we use the second-order phase-field theory, isogeometric collocation employs the 1-layer of Galerkin test functions at the patch boundary and collocation points in the interior of the patch (see Fig. 14b and 15b). It is easy to see that the cost is increasingly dominated by collocation points when the mesh is refined. Their number increases quadratically, while the number of quadrature points at the boundary increases only linearly and the corresponding elements become very thin boundary layers. The same is true for isogeometric collocation with adaptive Galerkin resolution of the fracture zone that in addition employs a band of NURBS test functions around the expected location of the horizontal crack (see Fig. 14c and 15c). In the coarser mesh, about 30% of the B´ezier elements require quadrature. Since the crack topology consists of a one-dimensional line, the number of quadrature points related to its Galerkin resolution increases linearly with mesh refinement, leading to the same cost reduction for sufficiently fine meshes. For the finer mesh shown in Fig. 14c, the number of elements with quadrature is already less than 5%. Figures 16a and 16b plot theoretical estimates for the number of point evaluations with increasing mesh size. We consider quadratic and quartic NURBS basis functions, as they are required for the collocation of the second-order and fourth-order phase-field fracture models, respectively. The upper bound of point evaluations corresponds to the Galerkin method, which requires quadrature in all elements. The lower bound corresponds to a basic point-collocation method that strictly uses one collocation point per basis function (with all disadvantages involved). We observe that isogeometric collocation using Galerkin test functions at the boundary quickly approaches the lower bound. When we adaptively resolve 28

(a) Initial.

(b) u = 6.75 × 10−3 mm.

(c) u = 7.86 × 10−3 mm.

Figure 17: Crack evolution obtained from collocation without Galerkin resolution of the fracture zone. The minimum length-scale that yields a stable analysis with the mesh of Fig. 14b is l0 = 5 h.

(a) Initial.

(b) u = 6.75 × 10−3 mm.

(c) u = 7.86 × 10−3 mm.

Figure 18: Crack evolution obtained from collocation with adaptive Galerkin resolution of the fracture zone. The minimum length-scale that yields a stable analysis with the mesh of Fig. 14c is l0 = 1.5 h.

the fracture zone by a band of Galerkin test functions, the convergence towards the lower bound is only slightly slower. Assuming sufficiently fine meshes, all versions of isogeometric collocation are about three times faster than Galerkin for quadratic NURBS, and about 15 times faster than Galerkin for quartic NURBS. Our comparison assumes that all schemes use NURBS discretizations with the same polynomial degree p. However, it needs to be noted that the Galerkin method can solve the second-order phase-field model with linear basis functions and the fourth-order phase-field model with quadratic NURBS, since it operates on the weak forms (26), (27) and (28) with reduced continuity requirements. 5.1.2. Crack resolution vs. mesh size As discussed in detail in Section 4 for the 1D test problem, the Galerkin method allows for a much smaller l0 than a collocation method that resolves the fracture zone with collocation 29

0.8 0.7

0.7

0.6

0.6

0.5

0.5

Load F [kN]

Load F [kN]

0.8

Galerkin Adaptive collocation−Galerkin

0.4 0.3

0.4 0.3

0.2

0.2

0.1

0.1

0 0

1

2

3 4 5 Displacement [mm]

6

7

0 0

8 x 10

−3

(a) Second-order phase-field model.

Galerkin Adaptive collocation−Galerkin

1

2

3 4 5 Displacement [mm]

6

7

8 x 10

−3

(b) Fourth-order phase-field model.

Figure 19: Load-displacement curves obtained with the Galerkin method and collocation with adaptive Galerkin resolution of the fracture zone.

points and Dirac δ test functions. Collocation that uses an adaptive resolution of the fracture zone with Galerkin test functions can use the same small l0 as the Galerkin method, since it accurately localizes the strains in the fracture zone. In turn, with the adaptive collocationGalerkin approach we can achieve the same accuracy with a coarser mesh. Figures 17 and 18 illustrate the progression of the crack at different deformation states for the two collocation variants, using the same NURBS mesh of 250×250 quadratic B´ezier elements that leads to quadrature and collocation points shown in Figs. 15b and 15c. Collocating in the fracture zone yields accurate results up to a limiting minimum length-scale of l0 = 5 h. Using a smaller parameter l0 induces spurious oscillations in stresses and strains that lead to the diffusive spread of spurious cracks over the complete domain. The adaptive resolution of the crack with Galerkin test functions yields accurate results up to a limiting minimum length-scale of l0 = 1.5 h, which is significantly smaller. Furthermore, we observed that collocating in the fracture zone relies on a finite parameter κ to stabilize the solution fields at the fully broken state. For the current problem, the adaptive collocation-Galerkin method yields accurate results even for κ = 0.0, while collocation with Dirac δ test functions in the fracture zone requires κ ≥ 0.001. Computations with the fourth-order phase-field model on a quartic NURBS mesh leads to the same observations and qualitatively similar results. Figures 19a and 19b plot the force obtained by integrating the shear stress over the top surface of the domain versus the imposed total displacements for the second-order and fourth-order phase-field computations, respectively. We observe that for the current mesh collocation with adaptive Galerkin resolution of the fracture zone achieves exactly the same load-displacement curve as the full Galerkin method. These results demonstrate that the adaptive collocation-Galerkin method successfully combines the advantages of the Galerkin method in terms of accurate crack resolution and the advantages of collocation in terms of reducing the number of point evaluations, offering 30

the best computational efficiency of all methods. For the current mesh of 250×250 quadratic B´ezier elements and the second-order model, the Galerkin method requires 250,000 point evaluations. The adaptive collocation-Galerkin method reduces the number of points evaluations to 72,000. To achieve the same degree of localization of the crack a collocation method that uses Dirac δ test functions and collocation points in the fracture zone would require a mesh of 830×830 B´ezier elements with approx. 700,000 point evaluations. We expect that the advantage of the adaptive approach will increase significantly, when we look at discretizations of three-dimensional problems with higher-order polynomial degree. 5.2. Single edge notched shear test As a second benchmark we consider the square plate with a horizontal notch and a horizontal shear loading [21–23]. The corresponding boundary conditions are illustrated in Fig. 20, and we use the same material parameters as before. To be consistent with reference calculations from the literature, we model the initial crack as a discrete discontinuity in the geometry. To this end, we discretize the square plate with four equal NURBS patches. This introduces C 0 lines that divide the domain in four equal quadrants, and the connectivity of the mesh is broken along the initial crack. As described in Section 3.4, multi-patch interface constraints are automatically taken care of when Galerkin test functions are used at the patch boundaries. The displacement at the top is again imposed in increments of ∆u=1×10−5 mm. We employ the second-order phase-field model and the isogeometric collocation method with adaptive Galerkin resolution of the fracture zone, and compute the solution on a quadratic NURBS mesh of 400×400 B´ezier elements with a length-scale parameter of l0 = 3 h. Our assumptions and parameters are thus consistent with computations reported by Miehe et al. [23] and Borden et al. [21], and we can use their results as reference solutions. Figures 21a, 21b and 21c plot collocation points (in red) and quadrature points (in green) that result from the resolution of the crack with Galerkin test functions at different loading states. The adaptive set of Galerkin test functions was generated by switching all NURBS functions of a B´ezier element from collocation to Galerkin if the phase-field falls below 0.8 u

0.5mm

0.5mm

0.5mm

0.5mm

Figure 20: Geometry and boundary conditions for the shear test of a single edge notched specimen.

31

(a) u = 9.0 × 10−3 mm.

(b) u = 11.0 × 10−3 mm.

(c) u = 13.6 × 10−3 mm.

Figure 21: Notched shear test: Points to be evaluated for different loading states when the fracture zone is adaptively resolved with Galerkin test functions. Quadrature points are plotted in green, collocation points in red (four patches, each discretized with 200×200 quadratic B´ezier elements).

(a) u = 9.0 × 10−3 mm.

(b) u = 11.0 × 10−3 mm.

(c) u = 13.6 × 10−3 mm.

Figure 22: Notched shear test: Crack progression at several loading states obtained with the meshes and point configurations shown above.

anywhere on the element domain. Figures 22a, 22b and 22c illustrate the phase-field approximation of the crack at the corresponding loading states. Using a reduced quadrature rule with p × p Gauss points the Galerkin method requires the evaluation of 640,000 quadrature points, while collocation with adaptive Galerkin resolution of the fully developed crack requires less than 200,000 point evaluations. Our approach thus reduces the formation and assembly cost by a factor three, while the accuracy can be expected to be comparable with Galerkin. For computations with higher-order splines, this ratio grows substantially. For example, using quartic NURBS on the same mesh the Galerkin method requires the evaluation of more than 2,5 million quadrature points, while the adaptive collocation-Galerkin method requires less than 350,000 point evaluations. 32

(a) u = 9.0 × 10−3 mm.

(b) u = 11.0 × 10−3 mm.

(c) u = 13.6 × 10−3 mm.

Figure 23: Notched shear test: Von Mises stress distribution plotted on the deformed domain.

Figures 23a, 23b and 23c plot the corresponding von Mises stresses on the deformed configuration. We do not show visualization points with a phase-field value below 0.1. For better visibility of the deformation, we multiply the displacements by factor five. We observe that the stresses are free of oscillations and reproduce the stress concentration at the reentrant corner where the crack is initiated. As soon as the crack is fully developed, the lower left part is stress-free and does not show any deformation. For better visibility of stresses away from the stress concentration at the re-entrant corner, we limit the stress scale to a maximum of six. We note that the highest stresses on the current mesh and with the current length-scale parameter are about one order of magnitude larger. 0.6

Load [kN]

0.4

0.2 Adaptive IGA collocation IGA Galerkin (Borden) FEA Galerkin (Miehe) 0.0 0.000

0.005

0.010

0.015

0.020

Displacement u [mm]

Figure 24: Force-displacement curves for the quasi-static shear test: Comparison of the curve obtained with the adaptive collocation method to curves obtained with standard C 0 finite elements and the isogeometric Galerkin method.

33

Figure 23 compares the load-displacement curve obtained with the adaptive collocation method to reference solutions reported by Miehe et al. [23] and Borden et al. [21]. We see that our curve is in good agreement with the reference curves that were obtained with standard C 0 finite elements and the isogeometric Galerkin method. 6. Summary, conclusions and outlook Phase-field fracture models are based on a diffusive representation of cracks, thus removing the requirement to numerically track discontinuities in the displacement field. In this paper, we started to explore the use of isogeometric collocation methods for the numerical solution of phase-field fracture models. The main objective is the speed up of phase-field fracture computations with respect to Galerkin methods by reducing the number of point evaluations during the formation and assembly procedure. To this end, we extended the hybrid collocation-Galerkin approach of De Lorenzis et al. [45] to the phase-field fracture case. We first reviewed the variational formulation of second-order and fourth-order phase-field fracture models, drawing on earlier presentations by Miehe et al. [22, 23, 26] and Borden et al. [21, 24]. We then derived the weighted residual form of the coupled phase-field and elasticity equations that constitutes the basis for isogeometric collocation. In the next step, we showed that the possibility of locally using NURBS test functions in the Galerkin sense instead of Dirac δ test functions supports the flexibility and consistency of isogeometric collocation at patch boundaries. First, using outer layers of NURBS functions at the patch boundary as Galerkin test functions provides a consistent way of weakly enforcing Neumann boundary conditions and multi-patch interface constraints. Second, Galerkin test functions at the patch boundary can naturally handle the more complex boundary integral terms that appear in the weighted residual form of the fourth-order phase-field model. We examined the accuracy of isogeometric collocation for one- and two-dimensional benchmark problems and compared it to the isogeometric Galerkin method. The results revealed that a collocation method that uses Dirac δ test functions and collocation points in the fracture zone requires a significantly finer mesh to accurately resolve a phase-field approximation of a crack with a given length-scale parameter l0 . This observation motivated the use of the hybrid collocation-Galerkin concept to adaptively resolve the crack with Galerkin test functions, while collocation is maintained in the bulk of the domain away from the fracture zone. We showed that isogeometric collocation with adaptive Galerkin resolution of the crack can resolve the same length-scale parameter on the same mesh as the Galerkin method. As the mesh is refined, the number of adaptive Galerkin test functions and the corresponding quadrature points along 1D patch boundaries and cracks increases only linearly, while the number of collocation points in the bulk of the domain increases quadratically. For sufficiently fine meshes the computational cost for the formation and assembly of discrete forms is therefore still dominated by collocation, and the number of point evaluations per basis function is close to the optimum of one. The adaptive approach thus combines the advantages of the Galerkin method in terms of accurate crack resolution on coarser meshes with the advantages of isogeometric collocation in terms of reducing the number of point evaluations. 34

The results of the present paper open up several directions for future work. The presented collocation-Galerkin methods can be explored for more complex fracture problems in three dimensions and in the fully dynamic setting. In explicit dynamics calculations, no equation solver is involved, and the computational efficiency is thus directly related to the number of point evaluations. We therefore anticipate a significant advantage of isogeometric collocation for explicit phase-field calculations. A further improvement of the current method can be achieved by combining its adaptivity based on a change from Galerkin to collocation with locally adaptive splines, e.g., hierarchical NURBS [79–83] or analysis suitable T-splines [84–88], that increase the resolution of the NURBS basis in the fracture zone. In this context the phase-field can be used as the indicator for local refinement at the evolving crack tip. Isogeometric collocation with adaptive Galerkin resolution of cracks profits from their dimensionally reduced geometric structure, that is, cracks are always represented by lines in two dimensions and surfaces in three dimensions. From a more general viewpoint, our approach adaptively resolves dimensionally reduced parts of the domain such that isogeometric collocation maintains its cost advantage for sufficiently fine meshes. We think that this collocation-Galerkin concept can be beneficially applied in many other problems where rapid changes in the solutions fields play a crucial role, such as advection-dominated transport with moving fronts or plate and shell structures with sharp boundary layers. Acknowledgments. We acknowledge the Minnesota Supercomputing Institute (MSI) of the University of Minnesota for providing computing resources that have contributed to the research results reported within this paper (https://www.msi.umn.edu/). M. Borden was supported by Grants from the Office of Naval Research (N00014-08-1-0992), the National Science Foundation (CMMI-1101007), and SINTEF (UTA10-000374) with the University of Texas at Austin. We also thank Jia-Liang Le for very helpful discussions on various aspects of fracture mechanics.

35

References [1] B. Bourdin, G.A. Francfort, and J.-J. Marigo. The variational approach to fracture. Journal of Elasticity, 91(1-3):5–148, 2008. [2] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 46(8):1319–1342, 1998. [3] J. Parvizian, A. D¨ uster, and E. Rank. Finite cell method: h- and p- extension for embedded domain methods in solid mechanics. Computational Mechanics, 41:122–133, 2007. [4] D. Schillinger, A. D¨ uster, and E. Rank. The hp-d adaptive finite cell method for geometrically nonlinear problems of solid mechanics. International Journal for Numerical Methods in Engineering, 89:1171– 1202, 2012. [5] D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. D¨ uster, and E. Rank. Small and large deformation analysis with the p- and B-spline versions of the Finite Cell Method. Computational Mechanics, 50(4):445–478, 2012. [6] M. Ruess, D. Schillinger, Y. Bazilevs, V. Varduhn, and E. Rank. Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method. International Journal for Numerical Methods in Engineering, 95(10):811–846, 2013. ¨ [7] M. Ruess, D. Schillinger, A.I. Ozcan, and E. Rank. Weak coupling for isogeometric analysis of nonmatching and trimmed multi-patch geometries. Computer Methods in Applied Mechanics and Engineering, 269:46–71, 2014. [8] D. Schillinger and M. Ruess. The finite cell method: A review in the context of higher-order structural analysis of cad and image-based geometric models. Archives of Computational Methods in Engineering, doi:10.1007/s11831-014-9115-y, 2014. [9] G. Alberti. Variational models for phase transitions, an approach via γ-convergence. In Calculus of variations and partial differential equations, pages 95–114. Springer, 2000. [10] A. Braides. Approximation of free discontinuity problems. Springer, 1998. [11] A. Chambolle. An approximation result for special functions with bounded deformation. Journal de math´ematiques pures et appliqu´ees, 83(7):929–954, 2004. [12] G.T. Camacho and M. Ortiz. Computational modelling of impact damage in brittle materials. International Journal of solids and structures, 33(20):2899–2938, 1996. [13] C. Daux, N. Mo¨es, J. Dolbow, N. Sukumar, and T. Belytschko. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 48(12):1741–1760, 2000. [14] R. De Borst. Numerical aspects of cohesive-zone models. Engineering Fracture Mechanics, 70(14):1743– 1757, 2003. [15] T.P. Fries and T. Belytschko. The generalized/extended finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering, 84:253–304, 2010. [16] N. Mo¨es, J. Dolbow, and T Belytschko. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1):131–150, 1999. [17] J.J.C. Remmers, R. De Borst, and A. Needleman. A cohesive segments method for the simulation of crack growth. Computational mechanics, 31(1-2):69–77, 2003. [18] X.-P. Xu and A. Needleman. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids, 42(9):1397–1434, 1994. [19] G. Lancioni and G.F. Royer-Carfagni. The variational approach to fracture mechanics. A practical application to the French Panth´eon in Paris. Journal of Elasticity, 95(1-2):1–30, 2009. [20] C. Kuhn and R. M¨ uller. A continuum phase field model for fracture. Engineering Fracture Mechanics, 77(18):3625–3634, 2010. [21] M.J. Borden, C.V. Verhoosel, M.A. Scott, Hughes T.J.R., and C.M. Landis. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217–220:77–95, 2012. [22] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of frac-

36

[23]

[24]

[25] [26] [27]

[28] [29] [30] [31] [32] [33] [34] [35] [36]

[37] [38]

[39]

[40] [41]

[42]

[43]

ture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 83(10):1273–1311, 2010. C. Miehe, M. Hofacker, and F. Welschinger. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199(45):2765–2778, 2010. M.J. Borden, T.J.R. Hughes, C.M. Landis, and C.V. Verhoosel. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering, 273:100 – 118, 2014. B. Bourdin, C.J. Larsen, and C.L. Richardson. A time-discrete model for dynamic fracture based on crack regularization. International Journal of Fracture, 168(2):133–143, 2011. M. Hofacker and C. Miehe. Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation. International Journal of Fracture, 178(1-2):113–129, 2012. M. Hofacker and C. Miehe. A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns. International Journal for Numerical Methods in Engineering, 93(3):276–301, 2013. A. Schl¨ uter, A. Willenb¨ ucher, C. Kuhn, and R. M¨ uller. Phase field approximation of dynamic brittle fracture. Computational Mechanics, 2014. C. Miehe and L.-M. Sch¨anzel. Phase field modeling of fracture in rubbery polymers. part i: Finite elasticity coupled with brittle failure. Journal of the Mechanics and Physics of Solids, 65:93–113, 2014. C. Miehe, F. Welschinger, and M. Hofacker. A phase field model of electromechanical fracture. Journal of the Mechanics and Physics of Solids, 58(10):1716–1740, 2010. Z.A. Wilson, M.J. Borden, and C.M. Landis. A phase-field model for fracture in piezoelectric ceramics. International Journal of Fracture, 183(2):135–153, 2013. A. Abdollahi and I. Arias. Phase-field modeling of fracture in ferroelectric materials. Archives of Computational Methods in Engineering, doi:10.1007/s11831-014-9118-8, 2014. C.V. Verhoosel and R. de Borst. A phase-field model for cohesive fracture. International Journal for Numerical Methods in Engineering, 96(1):43–62, 2013. J. Vignollet, S. May, R. de Borst, and C.V. Verhoosel. Phase-field models for brittle and cohesive fracture. Meccanica, doi:10.1007/s11012-013-9862-0, 2014. J.A. Cottrell, T.J.R. Hughes, and Y. Bazilevs. Isogeometric analysis: Towards Integration of CAD and FEA. John Wiley & Sons, 2009. T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194:4135–4195, 2005. F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, and G. Sangalli. Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 20(11):2075–1077, 2010. D. Schillinger, J.A. Evans, A. Reali, M.A. Scott, and T.J.R. Hughes. Isogeometric collocation: Cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations. Computer Methods in Applied Mechanics and Engineering, 267:170–232, 2013. F. Auricchio, F. Calabr` o, T.J.R. Hughes, A. Reali, and G. Sangalli. A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 249–252:15–27, 2012. T.J.R. Hughes, A. Reali, and G. Sangalli. Efficient quadrature for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 199:301–313, 2010. D. Schillinger, S.J. Hossain, and T.J.R. Hughes. Reduced B´ezier element quadrature rules for quadratic and cubic splines in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 277:1–45, 2014. F. Auricchio, L. Beir˜ao da Veiga, T.J.R. Hughes, A. Reali, and G. Sangalli. Isogeometric collocation for elastostatics and explicit dynamics. Computer Methods in Applied Mechanics and Engineering, 249–252:2–14, 2012. F. Auricchio, L. Beir˜ao da Veiga, J. Kiendl, C. Lovadina, and A. Reali. Locking-free isogeometric collo-

37

[44]

[45] [46] [47]

[48] [49] [50]

[51]

[52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69]

cation methods for spatial timoshenko rods. Computer Methods in Applied Mechanics and Engineering, 263:113–126, 2013. L. Beir˜ao da Veiga, C. Lovadina, and A. Reali. Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods. Computer Methods in Applied Mechanics and Engineering, 241– 244:38–51, 2012. L. De Lorenzis, J.A. Evans, T.J.R. Hughes, and A. Reali. Isogeometric collocation: Neumann boundary conditions and contact. ICES report 14-06, The University of Texas at Austin, 2014. H. Gomez, A. Reali, and G. Sangalli. Accurate, efficient, and (iso)geometrically flexible collocation methods for phase-field models. Journal of Computational Physics, 262:153–171, 2014. J.A. Evans, Y. Bazilevs, I. Babuˇska, and T.J.R. Hughes. n-widths, sup-infs, and optimality ratios for the k-version of the isogeometric finite element method. Computer Methods in Applied Mechanics and Engineering, 198(21–26):1726–1741, 2009. D. Grossmann, B. J¨ uttler, H. Schlusnus, J. Barner, and A.H. Vuong. Isogeometric simulation of turbine blades for aircraft engines. Computer Aided Geometric Design, 29(7):519–531, 2012. J.A. Cottrell, A. Reali, Y. Bazilevs, and T.J.R. Hughes. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 195:5257–5296, 2006. T.J.R. Hughes, A. Reali, and G. Sangalli. Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k -method NURBS. Computer Methods in Applied Mechanics and Engineering, 197:4104–4124, 2008. T.J.R. Hughes, J.A. Evans, and A. Reali. Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems. Computer Methods in Applied Mechanics and Engineering, 272(15):290–320. M.H. Sadd. Elasticity. Theory, Applications, and Numerics. Academic Press, 2009. A.A. Griffith. The phenomena of rupture and flow in solids. Royal Society of London Philosophical Transactions Series A, 221:163–198, 1921. G.R. Irwin. Fracture. In S. Fl¨ ugge, editor, Elasticity and Plasticity / Elastizit¨ at und Plastizit¨ at, volume 3 / 6 of Encyclopedia of Physics, pages 551–590. 1958. E.A. de Souza Neto, D. Peri´c, and D.R.J. Owen. Computational Methods for Plasticity: Theory and Applications. Wiley, 2008. C. Miehe. Comparison of two algorithms for the computation of fourth-order isotropic tensor functions. Computers & structures, 66(1):37–43, 1998. L. Kachanov. Introduction to continuum damage mechanics. Springer, 1986. D. Krajcinovic. Damage mechanics. Elsevier, 1996. J. Lemaitre and R. Desmorat. Engineering damage mechanics: ductile, creep, fatigue and brittle failures. Springer, 2005. L. Ambrosio and V.M. Tortorelli. Approximation of functional depending on jumps by elliptic functional via γ-convergence. Communications on Pure and Applied Mathematics, 43(8):999–1036, 1990. K. H¨ ollig. Finite Element Methods with B-Splines. Society for Industrial and Applied Mathematics, 2003. L. Piegl and W. Tiller. The NURBS Book. Springer, 1997. P.M. Prenter. Splines and Variational Methods. Wiley, 1989. D.F. Rogers. An Introduction to NURBS with Historical Perspective. Morgan Kaufmann Publishers, 2001. B.A. Finlayson and L.E. Scriven. The method of weighted residuals - a review. Applied Mechanics Reviews, 19:735–748, 1966. B.A. Finlayson. The Method of Weighted Residuals and Variational Principles. Academic Press, 1972. C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Fundamentals in Single Domains. Springer, 2006. R.D. Russell and J.M. Varah. A comparison of global methods for linear two-point boundary value problems. Mathematics of Computation, 29:1007–1019, 1975. D. Schillinger, J.A. Evans, F. Frischmann, R.R. Hiemstra, M.-C. Hsu, and T.J.R. Hughes. A collocated

38

[70] [71] [72] [73] [74] [75] [76] [77] [78] [79]

[80]

[81]

[82]

[83] [84] [85]

[86] [87]

[88]

C0 finite element method: Reduced quadrature perspective, cost comparison with standard finite elements, and explicit structural dynamics. International Journal for Numerical Methods in Engineering, submitted, 2014. C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Springer, 2007. R. De Borst, L.J. Sluys, H.-B. Muhlhaus, and J. Pamin. Fundamental issues in finite element analyses of localization of deformation. Engineering Computations, 10(2):99–121, 1993. Z.P. Bazant and J. Planas. Fracture and size effect in concrete and other quasibrittle materials. CRC press, 1997. M.S. Gelli and G.F. Royer-Carfagni. Separation of scales in fracture mechanics: From molecular to continuum theory via γ convergence. Journal of Engineering Mechanics, 130(2):204–215, 2004. T.J.R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, 2000. G. Strang and G.J. Fix. An Analysis of the Finite Element Method. Prentice-Hall, 1973. D.J. Benson, Y. Bazilevs, M.C. Hsu, and Hughes T.J.R. Isogeometric shell analysis: The ReissnerMindlin shell. Computer Methods in Applied Mechanics and Engineering, 199:276–289, 2010. D.J. Benson, Y. Bazilevs, M.-C. Hsu, and T.J.R. Hughes. A large deformation, rotation-free, isogeometric shell. Computer Methods in Applied Mechanics and Engineering, 200(13):1367–1378, 2011. D.J. Benson, S. Hartmann, Y. Bazilevs, M.C. Hsu, and T.J.R. Hughes. Blended Isogeometric Shells. Computer Methods in Applied Mechanics and Engineering, 255:133–146, 2013. A.V. Vuong, C. Giannelli, B. J¨ uttler, and B. Simeon. A hierarchical approach to adaptive local refinement in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 200(49– 52):3554–3567, 2011. D. Schillinger and E. Rank. An unfitted hp adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry. Computer Methods in Applied Mechanics and Engineering, 200(47–48):3358–3380, 2011. D. Schillinger, L. Dede’, M.A. Scott, J.A. Evans, M.J. Borden, E. Rank, and T.J.R. Hughes. An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces. Computer Methods in Applied Mechanics and Engineering, 249-250:116–150, 2012. D. Schillinger. The p- and B-spline versions of the geometrically nonlinear finite cell method and hierarchical refinement strategies for adaptive isogeometric and embedded domain analysis. Dissertation, Technische Universit¨ at M¨ unchen, http://d-nb.info/103009943X/34, 2012. M.A. Scott, D.C. Thomas, and E.J. Evans. Isogeometric spline forests. Computer Methods in Applied Mechanics and Engineering, 269:222–264, 2014. M.A. Scott, X. Li, T.W. Sederberg, and T.J.R. Hughes. Local refinement of analysis-suitable T-splines. Computer Methods in Applied Mechanics and Engineering, 213–216:206–222, 2012. Y. Zhang, W. Wang, and T.J.R. Hughes. Solid T-spline construction from boundary representations for genus-zero geometry. Computer Methods in Applied Mechanics and Engineering, 249–252:185–197, 2012. Y. Zhang, W. Wang, and T.J.R. Hughes. Conformal solid T-spline construction from boundary T-spline representations. Computational Mechanics, 51:1051–1059, 2013. E.J. Evans, M.A. Scott, X. Li, and D.C. Thomas. Hierarchical t-splines: Analysis-suitability, b´ezier extraction, and application as an adaptive basis for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2014. DC. Thomas, M.A. Scott, J.A. Evans, K. Tew, and E.J. Evans. B´ezier projection: a unified approach for local projection and quadrature-free refinement and coarsening of nurbs and t-splines with particular application to isogeometric design and analysis. arXiv:1404.7155, 2014.

39

Isogeometric collocation for phase-field fracture models

Jul 8, 2014 - [24] leads to higher regularity in the exact phase-field solution, ...... there is no need to exchange additional information across patches during.

1MB Sizes 0 Downloads 240 Views

Recommend Documents

A note on fracture criteria for interface fracture
criteria in (9) and (13) are employed for comparison and exhibited in Figure 10. ... more contact and friction induced in the specimen employed by Liechti and Chai (1992) ... elastic solution into (20). ... length parameter L to center the curve.

A note on fracture criteria for interface fracture
e-mail: [email protected]). Received 4 January .... arc in order to produce compressive residual stresses at the specimen edges. Residual curing stresses ...

Probabilistic Collocation - Jeroen Witteveen
Dec 23, 2005 - is compared with the Galerkin Polynomial Chaos method, the Non-Intrusive Polynomial. Chaos method ..... A second-order central finite volume ...

Simplex Elements Stochastic Collocation for ...
uncertainty propagation step is often computationally the most intensive in ... These input uncertainties are then propagated through the computational model.

Requirements for a Virtual Collocation Environment
Keywords. Virtual collocation, team work computer ... Boeing organizes its development programs as hierarchies of IPTs ..... socially before settling down to business, When annotating ... materials are the office applications used in the phase of ...

Collocation = Word partnership -
Mid = Middle : Midway. 9. Mis = Wrongly : Mistake. 10. Non = Not : Nonsense. 11. Over = Over : Overlook. 12. Pre = Before : Preview. 13. Re* = Again : Return.

Requirements for a Virtual Collocation Environment
Boeing Information and Support Services. P.O. Box 3707 M/S .... virtual collocation environment that will reduce or eliminate the need .... team finds the number for Bob and gives him a call. Bob answers. ... They decide that the new approach is best

a collocation approach for computing solar sail lunar ...
He notes that, as a general rule, the order of the integration scheme is ...... presented in a previous study by the authors.6 The data summary in Table 3 indicates.

2.5 Fracture Me.pdf
Addie, despierta, Addie-“. Me doy la vuelta con un gruñido y un estirón, froto ambos ojos con la palma de mi. mano. Es demasiado temprano para esta mierda.

2.5 Fracture Me.pdf
proses. Sudah sebagian tapi ini kan. berkembang, database e-commerce juga. berkembang terus. Jadi, kita sekarang. kembangkan database-nya, baik untuk. pemain dalam negeri ataupun yang. OTT,” kata Yon. z dit. Akhir September, Skema Pajak Toko Online

Comparison of Stochastic Collocation Methods for ...
Oct 30, 2009 - ... of 58800 volumes with a double cell size in the direction tangential ...... Factorial sampling plans for preliminary computational experiments, ...

Source Domain Determination: WordNet-SUMO and Collocation
Su, Professor Hintat Cheung and Professor Sue J. Ker for commenting on this paper. We would also like to thank Professor Chu-Ren Huang's. Shen-gen Project ...

A methodology for measuring interface fracture ... - Springer Link
Tel Aviv University, Department of Solid Mechanics, Materials and Structures, The ..... where ε and H are given in (3) and (8), respectively, R and t are the radius ...

PDF Online Fracture Mechanics: Fundamentals and Applications ...
PDF online, PDF new Fracture Mechanics: Fundamentals and Applications, Third Edition, Online PDF Fracture Mechanics: Fundamentals and .... worldwide.

hip fracture 16.pdf
mechanism of injury or bone pathology. For 2001, many of. the medical records and radiographs could not be retrieved,. and the community records including details of bone den- sity measurement and medication were unavailable. In. order to include mai

Fracture Me (ESPAÑOL).pdf
Page 1 of 60. FRACTURE ME TAHEREH MAFI. TRADUCCION HECHA POR: DOÑA TRUJI http://www.tumblr.com/blog/trujiyo. Page 1 of 60 ...