Adaptive isogeometric discretizations for diffuse modeling of discontinuities

In this contribution, we discuss the modeling of discontinuities using non‐conforming meshes. For weak discontinuities we propose to regularize the sharp material interface over a characteristic length scale. To define the material properties in the transition region, homogenization assumptions are used that fulfill the kinematic compatibility across the interface and the static equilibrium at the interface. For strong discontinuities we use a phase‐field model that is able to account for failure of adhesive interfaces and extend it to general heterogeneous materials. To improve the efficiency of the approach, adaptive isogeometric analysis is used to discretize the computational domain.


Introduction
In the numerical analysis of heterogeneous materials, weak and strong discontinuities arise in the field variables due to rapidly changing mechanical properties at material interfaces or due to propagation of cracks if a specific failure load is exceeded.
Details about the local topology of heterogeneous micro-structures are typically taken from imaging methods or are given in terms of random distributions. To avoid costly meshing processes, embedded domain methods, that represent the physical domain implicitly in a regular background mesh [1], can be used. While good convergence rates are achieved for homogeneous materials, stress oscillations occur in the heterogeneous case because weak discontinuities are modeled in terms of a continuous basis. Different non-standard discretization techniques exist to overcome this problem as the extended finite element method [2] or a method where for each heterogeneity a separate background mesh is used [3]. In [4], the vicinity of the sharp material interface is adaptively refined which leads to good convergence rates, but did not prevent the stress and strain oscillations.
Here, we follow the modeling approach in [4], but regularize the material interface to smooth the weak discontinuity. The regularization leads to a diffuse interface region of finite width , where the material properties are not defined. Instead of a simple interpolation, a homogenization is applied to compute the effective material parameters in the diffuse interface region, similar to phase field models [5,6].
To account for local damage and failure, the modeling of the heterogeneous material structure is extended by a phase-field model for brittle fracture [7]. In this way, cumbersome topological updates of the analysis mesh to track the crack path, are avoided. Crack initiation and propagation are accounted for by solving an additional scalar field problem that is coupled to the mechanical one.
To simulate interface failure between two materials, we combine the phase-field model with a local reduction of the critical fracture energy in the area of the regularized material interface [8]. The method considers possible interactions between the regularization of the interface and the phase-field and ensures crack propagation with the experimentally observed fracture toughness of the interface.
The variability of phase-field approaches comes at the cost of highly refined meshes that have to resolve strong gradients across the internal characteristic length scales. Hence, adaptive local refinement is advantageous for the efficient application of the mentioned approach. We therefore combine the modeling approach for the weak as wells as for the strong discontinuities with an adaptive spline-based discretizations based on THB-splines [9,10]. In case of weak discontinuities we apply an hadaptive refinement strategy, where the interface width and the element size h are reduced simultaneously. This allows for an improved convergence to the sharp interface solution and a good resolution of the local stress and strain fields [11].

Embedded domains and regularized embedded interfaces in linear elasticity
The problem of linear elasticity is defined by the partial differential equation −∇ · σ(u) = f in Ω phy and σ(u) · ν N = g on Γ N and u| ΓD = u D on Γ D , where f is the volume force , g the prescribed traction at the Neumann boundary Γ N and u D the prescribed displacement vector at the Dirichlet boundary Γ D . The stress tensor σ = C : (u) is given by the fourth order material tensor C and the strain tensor (u) = 1 2 ∇u + ∇u T .  x Ω δ (1) and boundary conditions f = 12x 2 kN mm −1 and uD = 0.01 mm.
In embedded domain methods, the problem above is solved numerically by embedding the physical domain Ω phy into a simple shaped domain Ω = Ω phy ∪Ω fic that can be easily meshed with a Cartesian grid of finite elements or a single rectangular B-spline patch, Fig. 1. To eliminate the influence of the fictitious domain Ω fic , the values of the material tensor in this domain are penalized by a small parameter C fict = C, where 1. If, in addition, Ω phy = n i Ω (i) consists of n different sub domains of different materials C (i) , the basis that is used to approximate the displacement field u should be C 0 -continuous across the material interface Γ ij to represent the weak discontinuity in u that is caused by the jump in C. Because of this contradiction to the paradigm of embedded domain methods, we relax this restriction by a regularization of the sharp interface Γ ij similar to phase field models. For this purpose, in a first step, an indicator functioñ is introduced to describe the heterogeneity of a bi-material body Ω phy = Ω (1) ∪ Ω (2) . In a second step, as illustrated in Fig. 1, the sharp interface Γ 12 is approximated by means of the Modica-Mortola energy [5] that leads to the smooth approximation of (2) defined by c = 1 2 1 + tanh xsd ,where x sd is the signed-distance function that describes Γ 12 with x sd = 0 and is a length scale parameter that controls the size of the regularized interface. Please note that in situations where the geometry description is obtained from imaging techniques, the regularized interface can be also derived directly from gray scale values.
In a third step, the material tensor C is redefined in the diffuse interface region in terms of the indicator function c to provide a smooth transition from the material tensor C (1) to C (2) . For this purpose, homogenization theory is applied at material points where both phases are present, c ∈ (0, 1), Fig. 1. A representative volume element with a straight interface between the phases Ω (1) and Ω (2) is defined in a reference coordinate system (x,ȳ,z). The indicator function c and 1 − c are interpreted as the dimension inȳ direction of the individual phases. Six deformation states for three-dimensional and three deformation states for two-dimensional problems have to be considered to compute the components of the resulting effective material tensor C(c) [11]. To take the orientation of the material interface into account, the material tensorC(c) is transformed from the reference coordinate system to the material coordinate system with basis {µ Γ , ν Γ , ξ Γ }. For the illustrated two dimensional case in Fig. 1, the normalized normal vector ν Γ of the interface can be computed from the indicator function ν Γ = ∇c ||∇c|| . Please note that the resulting homogenized material tensor C(c, ∇c) fulfills the kinematic compatibility across the interface and the static equuilibrium at the interface. Furthermore it contains the Reuss and Voigt type homogenizations as limiting cases and coincides with the homogenized material tensors proposed in [5,6].

Numerical example -Bi-material rod
The modeling approach above is tested and varified in the following numerical example of a bi-material, uni-axial rod under volume load, Fig. 2. The solution is compared against the embedded sharp interface representation under local h-refinement similar to [4]. For a better evaluation, we divide the physical domain Ω phy in a near field Ω δ and a far field Ω phy \ Ω δ . The near field is the domain that surrounds all interfaces Γ ij up to a distance of δ/2 from the interface, Fig. 1. The domain with length 1 mm, boundary conditions and material parameters are given in Fig. 2. The domain is discretized with THB-splines that are C 0 continuous at the boundary of Ω δ to allow for the jump in f and C p−1 -continuous across element boundaries in general. An h -adaptive refinement strategy is applied.
The convergence plot in Fig. 3 shows that the local h -refinement is able to reduce the relative H 1 -error in Ω phy (solid lines) to 10 −3 , an error range that is typical for embedded domain methods and suitable for engineering applications. This error is governed by the error in the near field and the solution converges with a rate that is independent from the chosen polynomial degree of the basis. However, in the far field (dashed lines), the solutions converge optimally due to the choosen homogenization scheme and the main error is bounded to the near field. This behavior is also mathematically justified [11]. Furthermore, the proposed approach avoids stress oscillations that occur for sharp interface representations, Fig. 3 on the right.
Please note, that no optimal convergence rates can be obtained for arbitrary shaped interfaces in multi-dimensional problems in the far field, due to a violation of the homogenization assumptions made in Section 2, that expect constant stress and strain states in the two different phases along the material interface. However, the total error can be clearly reduced under local hrefinement and an improved convergence rate is obtained compared to the embedded sharp interface under local h-refinement. Also oscillations in the stresses and strains are avoided using the proposed approach. [11] 3 Phase-field modeling of interface failure The phase-field approach allows for the proper modeling of crack initiation and propagation by solving an additional scalar field problem for the evolution of the crack phase which is coupled to the mechanical boundary value problem. To account for interface failure between two materials we combine the phase field model with a local reduction of the critical fracture energy in the area of the regularized material interface as proposed in [8]. The Helmholtz free energy for a brittle, linear elastic material reads where s is the phase-field variable indicating intact (s = 1) or fully broken (s = 0) material and s is the characteristic length scale that governs the width of the crack interface. Both, the bulk and the interface fracture toughness are taken into account by multiplying the fracture toughness G s with a scalar function α(x sd , i ), where x sd is a level set function. The function α regularises the formerly discrete adhesive interface over the length 2 i , similar as in Section 2. The influence of the interaction between i and s on the actual fracture toughness of bulk and interface is compensated following the the findings in [8].
To account, beside the heterogeneous fracture toughness, also for heterogeneous elastic properties, we combine the phasefield approach for interface failure with the diffuse modeling of material interfaces from Section 2. For this purpose we redefine the elastic strain energy density in terms of the homogenized material tensor www.gamm-proceedings.com  The specimen has an initial crack, indicated by the black solid line, and is subjected to a shear load. In contrast to the fracture process in a homogeneous body, the crack initiates at the material interface. To increase the efficiency of the computation, the mesh is pre-refined along the material interfaces and adaptively refined where the crack phase field s is smaller 0.5.

Numerical example -Crack propagation in heterogeneous materials
To illustrate the modeling possibilities of the phase-field model, we simulate the fracture process of a notched specimen under shear load. The domain of dimension 1 mm × 1 mm and boundary conditions are illustrated in Fig. 4. Three inclusions (c = 1) with glass fiber like material properties E gf = 79 GPa, ν gf = 0.17 and G gf = 2000 J m −2 are embedded in a matrix (c = 0) with thermoset like material properties E ts = 2.9 GPa, ν ts = 0.34 and G ts = 1000 J m −2 . The material interface has a reduced fracture toughness of G int = 500 J m −2 . The internal length scales are set to = 2 i = 4 s = 0.03 mm.
To efficiently resolve the steep gradients that arise in the phase-field variables, the domain is discretized with THB-splines of six levels. The mesh is pre-refined along the material interfaces and adaptively refined where the crack phase field s is smaller 0.5.
As illustrated in Fig. 4, after loading, two cracks initiate at the weakened material interfaces of two of the inclusions. Subsequently, they propagate into the matrix material, join each other and form a larger crack with an orientation of around 45°with respect to the loading direction. This is in contrast to the homogeneous case, where the phase-field crack initiates at the stress singularity at the initial crack tip [9].