Influence of Dynamic Loading on Fracture Behaviour of DCB Sandwich Specimen

Numerical simulations of dynamic fracture behaviour of a double cantilever sandwich beam subjected to uneven bending moments in plane conditions are carried out using the dynamic finite element analyses with the ABAQUSTM code. The strain energy release rate was evaluated by means of the finite element model developed within the two-dimensional (2-D) linear elastodynamic theory. This demonstrates the capability and the reliability of the finite element modelling as an extremely useful numerical tool for solving dynamic fracture mechanics problems. Also, the dynamic behaviour of fracture parameters and interface crack progression is discussed.


Introduction
Sandwich materials have opened new possibilities for designing in aerospace, automotive, civil, medical, sports and other industries. The sandwich concept bonding highly dissimilar materials layers are competitive in comparison to conventional metallic materials. This assembly structure has many advantages such as weight gain at high bending stiffness, sound insulation, protection against corrosion, etc. [1]. However, the performance of a tri-layered sandwich structure is highly influenced by both initial bonding processes (surface treatment, roughness, rheology, chemicals) and in-service parameters (loading, temperature, time, moisture). Also, a general non-isotropic nature of sandwich structures leads to discontinuous stress fields inside them. All these factors make premises for an inherent susceptibility of sandwich materials to damage. Herewith, debonding the face sheet from the core is the most often encountered issue among others. Failure processes in sandwich materials often initiate from such damage. Moreover, dynamic loading can speed up the debonding onset and cause its catastrophic propagation [2,3]. To ensure the reliability and damage tolerance of sandwich structures, an effective analysis method is needed to evaluate their performance mainly provided by the face sheet-to-core interface strength.
The interface strength can be quantified using the concept of interface fracture toughness [4]. In this respect, fracture specimens are used to supply useful information regarding the fracture resistance of the material [5,6]. The near-tip fracture parameters such stress intensity factors or strain energy release rates controlling the fracture event can be inferred from the specimens through experimental, analytical or numerical methods. However, while testing methods for fracture research of composite laminates are at a high level of maturity, the development of standard fracture tests for sandwich composites has not completed yet [7]. This is mainly due to a bi-material character of the face sheet/core interface and a non-symmetric geometry of the sandwich specimen, the analysis of which is required measurements of the energy release rate as a function of mixed mode loading [8]. An effective and relatively simple from standpoints of both fixing and loading conditions test method among many others, proposed during the last two decades, is a double cantilever sandwich beam (DCB) [9]. Motivated by the need to perform testing over as much as a possible wide range of mode mixity on the same specimen geometry, various modifications of DCBs have been considered as well. One of them is the DCB subjected to uneven bending moments (DCB-UBM) [10]. This test configuration for quasi-static loading and stationary cracks allows simple analytical considerations and, also, it produces a stable crack growth since the crack loading does not change with the crack length. However, when the loading and subsequent crack growth take place in time, the analytical procedures for the determination of fracture parameters involve complexities associated with inertia effects and stress waves as a result, numerical solutions are more suitable to analyse such problem.
Recent developments in the finite element method (FEM) related to fracture predictions such as cohesive finite elements (CEs) and the growth in computer power give a possibility to simulate the dynamic fracture behaviour in a reasonable amount of time. For example, finite element analyses with CEs have been successfully carried out to simulate the dynamic fracture of bi-material and composite plane plates in [11] and [12], respectively. The use of the CEs for modelling progressive impact-induced delamination in laminated and sandwich panels have been performed, e.g. in [13,14]. A propagation of penny-shaped debonding in sandwich plates under general dynamic loading and penny-shaped delamination in laminates under fatigue have been simulated with CEs in [15] and [16], respectively.
A comprehensive literature search shows that studies on the analysis of interface cracks in sandwich composites are mostly limited to quasi-static cases and isotropic constitutive materials. The aim of this paper is to develop the FEM model of a sandwich fracture specimen consisting of orthotropic face sheets and a core with high 'material mismatch' around the interfacial crack. For this aim, a DCB-UBM sandwich specimen is used and the dynamic fracture test of a sandwich material is simulated. Further, the influence of dynamic loading on the overall strength of the sandwich material is gained from this model for discussions.

Finite element methodology
A dynamic finite element framework incorporating cohesive elements is used in this study. A brief review of the numerical schemes applied to perform simulations and some specific notes concerning the cohesive element model utilized have been done below.

Linear momentum balance
Let us consider a DCB specimen as a two-dimensional continuum occupying a space V. The motion of the continuum is described by a displacement field u(x, t). At an arbitrary time t, the boundary of the medium ∂V is divided in a boundary ∂V u with prescribed displacement u and the boundary ∂V t with given tractiont . In addition to these, the continuum contains a crack presented by a cohesive surface referred to ∂V c = ∂V + c + ∂V − c . The the principle of virtual work in Lagrangian description without body forces can be stated as follows: where σ σ σ is the stress tensor; ρ is the density of material; Δ Δ Δ is the displacement jump across the flanks of the cohesive surface, and T = σ σ σ · n c along the cohesive surface ∂V c oriented with the normal n c in the reference configuration of the continuum.

Constitutive equations
The continuum is featured by two sets of constitutive equations. The first one defines the behaviour of the orthotropic bulk material in V and the second one is so-called tractionseparation law (TSL), i.e. a relation between the cohesive traction and the displacement jump which is defined on ∂V c . Since we are concerned with the interfacial crack problem, it is assumed that the crack lies between two parts of the continuum, i.e. V = V #1 ∪ V #2 , where #1 and #2 denote those parts which are different in geometry and material properties. Let the material of each parts obey the general 2-D Hooke's law in plane stress [17]: where In plane strain the coefficients of the reduced compliance matrix should be replaced by � S 33 for i, j = 1, 2, 6. The strain tensor ε ε ε associated with the stress tensor σ σ σ is defined in accordance with the assumptions of infinitesimal deformations as ε ε ε = 1 The cohesive element is idealized by a pair of separate top and bottom surfaces, which are discretized in such way that in the finite element mesh the nodes at the interface have the same coordinates, but different node numbers (Fig. 1a). A TSL of the element allows predicting both the onset of the softening process at the crack tip as a result of a strength-based analysis and the crack propagation, conditions of which rely on fracture mechanics criteria. A typical bilinear cohesive law for a single fracture mode is presented in Fig. 1b. The law contains an initial linear region defined by a penalty stiffness k and the softening part starting from the value Δ 0 , where the traction reaches a maximum normal/shear cohesive value T 0 and, then, evolving linearly till Δ f , where complete failure occurs. The irreversibility conditions are assumed to be rendered by unloading to the origin. The area under the lines being the work done per unit area for complete fracture defines the strain energy release rate. Analytically, for each fracture mode (i = I, II, III) the bilinear TSL can be presented as [18]: Here, a damage variable D can be calculated as a function of current separation between the cohesive element faces, i.e.
In the case of mixed modes, an effective displacement jump II, III is the number of modes involved, is introduced and the damage initiation and evolution criteria are to be formulated in terms of interaction between the fracture parameters of each mode, Fig. 2c. In this regard, the equivalent mixed mode separations at damage onset Δ 0 m and failure Δ f m are to be defined. Following [18] the damage initiation based on the quadratic stress initiation criterion takes the form: is a mixed mode ratio, whereas the damage propagation relying on the Benzeggagh-Kenane (B-K) fracture toughness criterion reads Here Δ = s � Δ 2 II + Δ 2 III , G T is the total ERR and η is a parameter obtained by curve-fitting the fracture toughness of mixed mode tests. Once mixed mode separations are known the mixed mode damage parameter D m can be calculated identical to the expression for D i using Δ m , Δ 0 m and Δ f m instead the pure mode components there.

Contact and friction conditions
In general, the crack flanks are assumed to be traction free. On the other hand, upon cracking the cohesive top and bottom surfaces may be subjected to unilateral contact and friction. Moreover, dynamic loading may cause to interact the crack flanks [19,20]. Thus, contact and friction conditions should be implemented into the model to provide accurate simulations. Decomposing the contact traction vector t c over presumed contact surfaces within ∂V c into normal t N and tangential t T components, each of them can be defined in terms of gap g N and slip g T functions of the displacements u, respectively. Then, the impenetrability and friction constraints referred to the Karush-Kuhn-Tucker conditions are stated as [21]: and respectively. Here, t N is the scalar quantity of the normal traction component, i.e. t N = t N n c ; τ crit is a threshold of tangential contact traction due to the tangential slip. In the case of the Coulomb friction model, this value is expressed as τ crit = μt N , where μ is the coefficient of friction.
In the presence of contact and friction, an appropriate term referred to as the work of contact forces should be added into the variational equality (1). Hence, it reads

Finite element discretization
Following the FEM framework, the actual continuous model is idealized as an assemblage of finite elements interconnected at nodal points. Consequently, (1) is transformed to the discrete system of equations of motion with respect to nodal degrees of freedom, e.g. displacements, at time t as follows [22]: where {R int }, {R ext } , {R coh } and {R cont } are the vectors attributed to the nodal internal, external, cohesive and contact forces, respectively, calculated using the corresponding integrals in (8) where velocities and displacements are calculated as{U} t+ 1 2 Δt ={U} t− 1 2 Δt + Δt{U} t and {U} t+Δt = {U} t + Δt{U} t+ 1 2 Δt , respectively, [22].
[ M] is a lumped mass matrix obtained by the transformation of the consistent mass matrix [M] for the purpose of efficiency.
The dynamic implicit analyses are carried out by using the implicit Hilber-Hughes-Taylor (HHT) temporal integrator in ABAQUS/Standard. In accordance with the HHT scheme, the equations of motion (8) at a particular time point t n+1 = t + Δt can be rewritten in the form: where displacements and velocities at the time point t + Δt are approximated by [22]. The solution of the implicit analysis from t to t + Δt is updated incrementally within the well-known Newton-Raphson iterative scheme for finding the roots of (10).

Dynamic energy release rate
In LEFM the Rice's J-integral is identical to the ERR. The generalization of this fundamental concept on an elastic solid with a crack advancing straightway along the x 1 -axis direction under dynamic conditions can be expressed as [23]: where the path Γ is an arbitrary contour surrounding a crack tip; n is an outward unit normal of Γ; W is the strain energy density and T is the kinetic energy density at a material point, x. Under a steady state crack growth condition, the dynamic integral (11)  and corresponds to the instantaneous energy release rate for any crack configuration including the interface crack between two dissimilar orthotropic materials [23].
The 'domain integral formulation' allows a simple FEM computation of the dynamic J �integral. With using a weight function q 1 (x), the line integral (11) is transformed to a domain integral, i.e. the dynamic ERR can be evaluated by computing the expression [23]: where A is the domain enclosed by the contour Γ, an arbitrary contour C with unit normal m, which embraces Γ and the surfaces of crack flanks, C + and C − between the two contours. The weighting parameter q 1 is chosen as a smooth function of x which takes the values from zero on the C-contour to unity on Γ. One of advantages of this test method is that the specimen allows loading the crack tip by a variety of mode mixities by changing the moment ratio M R = M 1 /M 2 [10]. For this specimen, the steady state ERR can be determined analytically from the specimen geometry, elastic properties and applied external moments inducing a pure bend. Following [24], the J-integral calculated along the outer boundaries of the specimen (Fig. 2b) leads to the expression:

Numerical results
A 2-D finite element model of the DCB-UBM specimen is developed eight-node reduced integration plane strain finite elements (CPE8R) available in ABAQUS [26]. The mesh contains a refinement near the crack-tip region as shown in Fig. 3. The debonded region of the specimen is modelled by a real gap of h 1 100 along the damaged interface. The contact and friction conditions similar to (5) and (6) are introduced between the faces of finite elements along the pre-cracked interface. The contact behaviour under the assumptions of small displacement kinematics was assumed to be governed by the 'hard contact' model with frictionless conditions [26]. In the case of the explicit dynamic analysis (9) the contact constraints were resolved using the kinematic predictor corrector method [19,26], while the penalty contact algorithm was used for tracking contact in the case of the implicit dynamic analysis (10) as described in [20,26]. In the calculations, bending moments or equivalent rotations have been applied to the arms of DCB specimen at the points of the arms' neutral axes, Fig. 3. Coupling kinematic constraints between the nodes of the arm edge and the points of neutral axis are used to enforce equal rotation of the entire edge.  The J-integral approach is a built-in option in ABAQUS fracture analysis [26]. In the analysis, five to six domains were typically chosen to compute G. For dynamic analyses, the dynamic ERR is calculated using (12), where the mechanical fields at an instant of time t account for inertia forces. The debonding propagation was analysed using four-node 2-D cohesive elements (COH2D4) available in ABAQUS [26]. The elements with the bilinear TSL (3) were introduced at the uncracked ligament between the face sheet and the core. Then, at any loading step, the fracture criterion is checked, if it is satisfied, the separation occurs and the deboning advances on an arbitrary length and the contact conditions are activated between newly formed crack surfaces. Herewith, the damage initiation associated with the debonding onset was tracked by satisfying the quadratic stress criterion, and the damage evolution referred to the debonding propagation was modelled using the mixed-mode energybased B-K fracture criterion. It also was assumed that a presumed crack path is confined only the face/core interface, i.e. no kinking is considered. The calculations have been performed under displacement-controlled loading by applying prescribed rotations to the specimen arms.
To verify the developed finite element model, the static ERRs, G s , for a variety of moment ratios have been computed using the J-integral method. The bending moments were either rotated in opposite directions or co-rotated, but they induced nearly the same ERR in all cases. A good agreement between the numerical results and those found by the analytical formula (13) has been achieved as shown in Table 2.

A stationary crack under dynamic loading
First, the effects of impulse loading on the transient dynamic ERRs are examined. The bending moments are applied to the sandwich specimen arms as impulses. Several types of the impulses of different forms and durations have been examined in the calculations. The influence of the impulse loading on the transient dynamic ERRs, G d , is presented in Fig. 4. One can see that the transient dynamic ERRs exceed their static values for all cases of loading, and both the impulse duration and the impulse form remarkably affect the dynamic ERR.
Thereafter, the specimen is subjected to sinusoidal moment loading with driving frequencies which are either fraction or multiplier, ζ = Ω ω 1 , of the fundamental frequency of the same intact sandwich beam ω 1 , which was computed as 150.37 Hz. Fig. 5 shows that the behaviour of long-term dynamic ERRs highly depends on the driving frequency. The amplitude of the dynamic ERRs tends to increase with increasing the driving frequency. Also, the oscillation form of the ERRs changes from a regular one at the lower frequency to irregular one at the higher frequency as illustrated in Fig. 5.

A crack growth under dynamic loading
First, the case of quasi-static bending induced by prescribing rotations is simulated and forcedisplacement curve is extracted from the finite element solution, Fig. 6. From this illustration, one can see that the interface crack propagates in a stable manner as expected for this test method. Also, the stress state near the crack tip during the deboding growth involves both the normal and shear stress tensor components as shown in Fig. 6b.
Next, the impulsive rotations of different durations, but the same rectangular form are applied to the DCB UBM specimen. For all cases of loading, the same rotation magnitude is hold. In Fig. 7a the simulations demonstrate that the total debonding extension increases with decreasing the impulse duration and for the shortest impulse load a complete disintegration of the specimen occurs. In the cases of loading without fracture, the debonding extends with   Fig. 7b. Finally, the sandwich specimen is subjected to harmonic rotations with a given driving frequency accepted as high as 3/2 of ω 1 . The analysis was limited by 3000 increments and lasted at least 100 increments after reaching the steady state oscillation regime. The results of simulation of the debonding behaviour under harmonic loading are shown in Fig. 8. In this case the debonding propagates in a stick-slip manner, i.e. by jumping from one debonded state to another one as demostrated in Fig. 8a. Moreover, the simulations revealed that the interface crack was intensively growing when the detached vibrated face sheet had the form of a concave downward curve, Fig. 8b. This is associated with a mode II dominated regime.

Conclusions
The problem of the dynamic debonding is considered and main features of the problem are highlighted. The finite element model for simulating dynamic debonding of sandwich DCB-UBM specimens with numerical tools available in the ABAQUS code is presented. The influence of dynamic loading on the strain energy release is found out. The CZM with bilinear TSL assumptions for the mixed mode debonding initiation and evolution, which was embedded into the cohesive elements, is used to simulate debonding propagation in the DCB-UBM test under both quasi-static and dynamic, involving impulse and harmonic loads.