Obtaining initially rigid cohesive finite element models that are temporally convergent

Abstract

Based on a time continuity condition proposed in our earlier work, we develop a fairly general framework, free of regularization issues, for producing initially rigid cohesive finite element models that are temporally convergent in an explicit dynamics setting. These models are adaptive in the sense that an interface is inactive until the traction across it reaches a critical level. We require that the traction across the interface be a continuous function of displacement, and therefore of time, i.e., the traction before and after activation should be the same. We demonstrate with examples that failure to satisfy this condition leads to non-optimal convergence or failure of convergence as the time step tends to zero. We also investigate the predicted crack path produced by a time continuous model versus a discontinuous model in a simulation of a laboratory experiment. The time-continuous model produces crack paths that are stable with variations of the time step and are close to laboratory measurements.

Introduction

Fracture processes can generally be represented either as weak or strong discontinuities. Cohesive zone modeling, which was originally proposed by Dugdale [6], Barenblatt [2] and Rice [26], is in the latter category. In such a model, the separation of bulk material is resisted by cohesive forces, governed by the corresponding cohesive constitutive model.

It is natural to apply phenomenological cohesive modeling to problems where defined interfaces exist and, indeed, this approach has been successfully used in problems concerning static [23], [32], [36] or dynamic [28], [34] bi-material interfacial debonding and fracture; delamination of composite materials [1], [12]; dynamic fracture in fiber reinforced composites [7]; intergranular dynamic fracture under impact loading [8], [9], [40]; void nucleation [19] and many other settings. Many of these numerical simulations were performed using the finite element method by introducing cohesive interface elements. In such analyses, the fracture path is assumed a priori and can usually be justified by the nature of the problem (e.g. physical weak interfaces in delamination problems) or by experience gained from experiments (e.g. observations that fracture occurs frequently along intergranular polycrystal boundaries).

Many cohesive constitutive models, usually expressed in the form of potential functions, are available in the literature [18], [20]. A widely used cohesive model is of a bilinear form, pure normal and pure shear realizations of which are shown in Fig. 1. In the figure, t_c, βt_c represent critical peak normal and tangential tractions, and δ_nc, δ_tc are the corresponding critical opening displacements; δ_nf, δ_tf are failure opening displacements, at which the traction drops to zero, and k_c is a penalty parameter to resist material interpenetration between surfaces. Note that the parameters in the pure normal and pure shear cases are related in such a way that the area under the curve in the first quadrant is equal to the modes I and II fracture energies, G_Ic and G_IIc, respectively. Throughout this paper, we call this type of cohesive model initially elastic (also called intrinsic [16]) because of its finite initial stiffness. Due to its simplicity and ease of incorporation in the finite element method, this model is often used in the literature.

In the past decade, there have been increasing applications of cohesive interface finite elements in the area of dynamic fracture including, among others, branching [38]; and fragmentation [4], [18], [24], [41]. In most of these simulations, the crack pattern is not known in advance. This precludes the possibility of prepositioning cohesive surfaces. Instead, every edge of the bulk elements must be considered as a potential fracture surface and the crack propagation path must be resolved as part of the solution of the governing equations.

Two different approaches are taken to this end. In the first approach, e.g., [18], [38], every edge of each bulk element is the site of a cohesive surface. In such models every element has its own copy of a connecting node and the copies are allowed to separate according to the cohesive constitutive law. The number of degrees of freedom is thus greatly increased, e.g., a four- to six-fold increase in the analysis of [38]. The increase will be even greater in three-dimensional analyses. The more fundamental problem is that, as the finite element discretization is refined, the effective modulus of the material is non-physically decreased, depending on the level of refinement. In fact, a simple one-dimensional analysis [10] shows that the effective modulus of a system of elastic blocks connected by cohesive surfaces is

$$E_{\mathrm{eff}}=E\left[1-\frac{1}{1+k_{\mathrm{i}}h/E}\right]\text{,}$$

where E is the elastic Young’s modulus of the blocks, k_i is the initial stiffness of the cohesive model and h is the linear spacing of cohesive surfaces. Obviously, as h approaches zero, the material effectively loses all its stiffness. The degree of reduction is expected to be more serious in 2D and 3D cases. This means that the wave speeds of the material in dynamic simulations, which are the focus of this paper, will not be captured correctly. One rudimentary solution is to pick k_i as high as possible, but this will introduce very stringent requirements on the time step size in explicit dynamics.

An analysis of the initial slope is carried out by Espinosa and Zavattieri [8], [9] in a two-dimensional setting. They reach the same conclusion, namely, that the initial slope k_i must be much larger than E/h. In the problem they consider, the mesh is much finer than the grain size, and they use cohesive elements only on grain boundaries. In this setting, the spacing h between cohesive elements does not decrease as the mesh is refined, so the initial slope may be selected independent of mesh refinement. If one is interested in the situation in which every element boundary is a potential interface location, then the choice of k_i is not immediately obvious and, in general, requires some heuristics and experimentation. Furthermore, it will be a function of the mesh size.

Alternatively, an adaptive approach [4], [21] is followed, in which cohesive surfaces are only inserted as needed, i.e., the effective initial stiffness of the cohesive model is infinite. We call this type of model initially rigid or simply rigid (also called extrinsic [16]), see Fig. 2. With this model, the effective Young’s modulus of the bulk elements before any crack growth occurs is independent of mesh refinement. Hence, the correct wave speeds can be captured at this stage.

Although the rigid model can eliminate some of the drawbacks of the elastic model, special considerations in its finite element implementation are essential to its numerical behavior. Papoulia et al. [22] proposed that a specific consideration, called the time-continuity condition, is required to obtain correct time convergence results in dynamic simulations using rigid cohesive models, e.g., second-order convergence in time with the central difference method. In this paper, we propose a general framework for producing initially rigid cohesive finite element models that are temporally convergent. The procedure is free of regularization issues present in [22] and in other works [10]. We focus on explicit solution techniques, specifically a central difference time stepping scheme for explicit dynamics.

The time-continuity condition is described in Section 2. Section 3 presents our newly proposed model. Section 4 describes our rigid cohesive modeling methodology. Two numerical examples of a beam bending problem and of a Chevron-notch specimen test are investigated in Section 5, including convergence studies and comparisons with experimental data. Finally, we offer our concluding remarks regarding rigid cohesive modeling in Section 6.