Dynamic debonding in layered structures : a coupled ALE-cohesive approach

A computational formulation able to simulate crack initiation and growth in layered structural systems is proposed. In order to identify the position of the onset interfacial defects and their dynamic debonding mechanisms, a moving mesh strategy, based on Arbitrary Lagrangian-Eulerian (ALE) approach, is combined with a cohesive interface methodology, in which weak based moving connections are implemented by using a finite element formulation. The numerical formulation has been implemented by means of separate steps, concerned, at first, to identify the correct position of the crack onset and, subsequently, the growth by changing the computational geometry of the interfaces. In order to verify the accuracy and to validate the proposed methodology, comparisons with experimental and numerical results are developed. In particular, results, in terms of location and speed of the debonding front, obtained by the proposed model, are compared with the ones arising from the literature. Moreover, a parametric study in terms of geometrical characteristics of the layered structure are developed. The investigation reveals the impact of the stiffening of the reinforced strip and of adhesive thickness on the dynamic debonding mechanisms.


INTRODUCTION
uring the last decades, layered structures in the form of laminates or thin films have employed extensively in many engineering fields, ranging from nano to macro scale applications.Typically, such materials are formed by strong layers and weak interfaces, in which internal material discontinuities may evolve, producing relevant loss of stiffness [1].Moreover, the crack evolution is strongly affected by the time rate of the external loading, which typically produces high amplifications of the fracture parameters.As a matter of fact, the measured crack tip speeds, during crack propagation, are relatively high, ranging also close to the Rayleigh wave speed of the material [2,3].Therefore, in order to predict the interfacial crack growth, models developed also in a dynamic framework are much required.

D
In order to simulate debonding phenomena in layered structures, several approaches have been proposed in the literature.However, among the most important ones, Fracture Mechanics (FM) and Cohesive Zone Model (CZM) are widely utilized in practice [4].In FM, the total energy release rate and its individual mode components need to be evaluated to predict delamination growth.For general configurations, the energy release rates can be computed by using a very accurate mesh of solid finite elements and the Virtual Crack Closure Method (VCCM) [5].Such models calculate the energy release rate as the work performed by the internal traction forces at the crack faces during a virtual crack advance of the tip.Moreover, in dynamic Fracture Mechanics, the VCCM is applied by using the modified form, in which the ERR, during the time evolution, is evaluated by the product between the reaction forces and the relative displacements at the crack tip and at the nodes closer to the crack tip front, respectively, [6,7].The prediction of the energy release rate is strictly dependent from the mesh discretization of the crack tip.Alternatively, cohesive models propose an easy way to simulate debonding phenomena including also crack onset.Distributed or discrete interface elements are introduced between continuum elements based on traction separation damage laws.However, such modeling is strictly dependent from the mesh discretization since the direction of crack propagation is restricted by the element size and orientation adopted by the user [8].Moreover, the presence of an initial finite stiffness may produce in brittle solids an excess of compliance and, in those cases in which a high stiffness is introduced, spurious traction oscillations [9].Such problems may be partially circumvented by introducing a very fine discretization at the crack tip front, to obtain a high resolution of the characteristic fracture length of the interface [10].However, the resulting model is affected by computational complexities, because of the large number of variables and nonlinearities involved along the interfaces.In order to avoid such problems, combined formulations based on fracture and moving mesh methodologies are proposed [11,12].In particular, the former is able to evaluate the variables, which govern the conditions concerning the crack initiation and growth, whereas the latter is utilized to simulate the evolution of the crack growth by means of ALE formulation [13,14].It is worth noting that the use of moving mesh method, combined with regularization and smoothing techniques, appears to be quite efficient to reproduce the evolution of moving discontinuities.However, existing models based on ALE and FM are based on a full coupling of the governing equations, arising in both structural and ALE domain.In this framework, material and mesh points in the structural domain produce convective contributions and thus nonstandard terms in both inertial and internal forces.In the proposed formulation, the use of a weak discontinuity approach avoids the modification of the governing equations arising from the structural model and thus a lower complexity in the governing equations and the numerical computation is expected.Despite exiting numerical methodologies based on pure CZM, the present approach reduces the nonlinearities involved in the governing equations to a small region containing the process zone, leading to a quite stable and efficient procedure to identify the actual solution in terms of both crack initiation and evolution.In order to verify the consistency of the proposed model, comparisons with existing formulations for several cases involving single and multiple delaminations are developed.The outline of the paper is as follows.At first, the formulation of the governing equations for the ALE and interface approach is presented and, subsequently, the numerical implementation of the finite element model is reported.Finally, comparisons and parametric results to investigate static and dynamic behavior of the debonding phenomena are proposed.

THEORETICAL FORMULATION OF THE MODEL
he proposed model is presented in the framework layered structures, in which thin layers are connected through adhesive elements.In particular, a shear deformable beam model with a proper number of mathematical elements along the thickness direction is utilized to reproduce a high order zig-zag kinematic (Fig. 1).However, each layer is connected to the adjoining ones by means of imperfect interfaces, in which debonding phenomena may affect the adhesion between layers introducing material discontinuities and traction forces along normal or sliding directions.In order to simulate debonding phenomena, a fundamental task to be achieved is to identify the position, in which the onset of interfacial mechanisms is produced and subsequently to simulate the evolution of the cracked length.The theoretical formulation is articulated into two steps, which are presented separately.At first, when the crack onset condition is not satisfied, the moving mesh elements are inactive and only the structural problem is considered.Subsequently, when the debonding process is triggered, the ALE elements are activated introducing moving traction forces which follows the process zone length.

Governing equations
The governing equations of the FE model are derived by means the principle of virtual works.In particular, for the general dynamic problem presented in Fig. 1 it is required that the total virtual work is stationary: where T  is the virtual work of the inertial forces, U  is the work of the internal forces and the tractions across the interfaces and W  is the work of the external forces.According to the first-order transverse shear deformable laminate theory and multilayered approach, the variational form of the governing equations can be expressed by means of the following expressions: where the subscripts l=1,..,N and i=1,..,N-1 indicate the numbering of the layers (N) and the interfaces (N-1),   , , is the cohesive interface displacement jump vector,  and 0  are the mass and polar mass per unit length of the layer and l f  and l p  ,with , are the per unit volume and area forces acting on the l-th layer, respectively.

The cohesive interfaces
The cohesive interfaces are introduced between the sublayers in which the crack initiation could be potentially activated.
The crack onset definition is described by means of a mixed crack growth, which is a function of the fracture variables, coinciding with the ratio between ERR mode components and corresponding critical values, as follows: ( ) ( ) ( ) where i represents the generic i-th interface in which debonding phenomena may occur, r is the constant utilized to describe fracture in different material and( ) are the total area under the traction separation law, whereas ( )   The interfaces moving mesh strategy Until the crack onset condition is not satisfied the ALE equations are inactive and the position of the computational mesh points is expressed in the Fixed Material Frame FM  identified by the x1-x 2 coordinates (Fig. 2a), which coincides, at this stage, with the Moving (M) Frame M  (X 1 -X 2) .It worth noting that the difference between FM  and M  is mainly in the mesh element number involved in the numeric model.In particular, in FM  a coarse mesh is required to identify the of the crack onset position, whose accuracy is verified introducing more finite elements by using a remeshing procedure in the M  configuration (Fig. 2b).This procedure leads to the evaluation of the crack initiation position 1 i X which sets to zero Eq. 3. The mathematical description of the moving mesh modeling is defined by a mapping operator Φ , which relies a particle in a Fixed Referential Frame R  and the one in current moving coordinate system (Fig. 2c).The mesh motion in terms of displacement field, at the i-th interface, is described as the difference between material 1 i X and the referential coordinates ( ) , where i i B h W = ´ represents the region in which the debonding mechanisms are produced.In order to reduce mesh distortions, produced by the mesh movements, rezoning or smoothing equations are introduced to simulate the grid motion consistently to Laplace based equation which is, in the case of one dimensional domain for both Static (S) and Dynamic (D) cases, defined on the basis of the following relationships [15]: Eqs.( 5) should be completed by means of boundary equations to reproduce the crack tip motion on the basis of the assumed crack growth criterion, namely i f g .In particular, for a fixed position in which the crack initiation occurs, different boundary conditions should be introduced to enforce internal or external debonding mechanisms.In particular, once the position of the crack onset is determined, i.e.at 1 1 i X X  , a geometrical debonding with length equal to 2 is introduced in the numerical model, producing two potential finite crack tips in which debonding phenomena can be triggered (Fig. 2b).In order to described the evolution of debonding phenomena the following boundary conditions should be introduced, which are completed by the Kuhn-Tucker optimality conditions: However, additional relationships are required for the ALE formulation to reproduce the crack tip motion and boundary conditions.In particular, the displacements of the computational nodes are assumed to be zero at the boundaries of the structure, whereas small portions, close to the crack tip for the left and right debonding mechanisms, namely i   and i   , are assumed to be moved rigidly, enforcing the computational nodes at the extremities to have the same displacements.This choice ensures that the NLs involved in the debonding mechanisms are constrained to a small portion containing the process zone, reducing the total complexities of the model.Consequently, the following boundary conditions should be considered in the analysis (Fig. 2b): It is worth noting that   can be considered as variable quantities for each crack path, that should be determined during the crack evolution.In particular, from the physical point of view, they represent the portion in which the traction separation laws are distributed.Since those regions are assumed to be moved rigidly, by means of the ALE strategy, the nonlinearities involved by the traction forces may be reduced to a small region close to the crack tip, avoiding as a result spurious and oscillatory effects typically documented in pure CZMs.From the numerical point of view, displacements of the debonding regions are determined on the basis of the values of the fracture function at their extremities by enforcing that during the crack growth a null value of the fracture function on the debonding region is achieved.Therefore, by using a linear approximation function along the debonding region, the current crack tip displacement can be expressed by means of the following relationship:

NUMERICAL IMPLEMENTATION
overning equations introduced in previous section are formulated by means of a numerical formulation based on the Finite Element (FE) approach.In particular, the derivation of the FE starts from the principle of virtual works in terms of displacements, integrating the equations on the volume of the elements.A Lagrange cubic approximation is adopted to describe both displacement and rotation fields, whereas linear interpolation functions are utilized for the axial displacements.Moreover, for the variables concerning moving mesh equations, quadratic interpolation functions are assumed to describe the mesh position of the computational nodes.The proposed approach takes the form of a set of nonlinear differential equations, whose solution is obtained by using a customized version of the finite element package Comsol Multiphysics combined with MATLAB script files [16].The model can be solved in both static and dynamic framework, taking into account the time dependent effects produced by the inertial characteristics of the structure and the boundary motion involved by debonding phenomena.In both cases, since the governing equations are essentially nonlinear, an incremental-iterative procedure is needed to evaluate the solution [17].In the case of static analysis, the resulting equations are solved by using a nonlinear methodology based on Newton-Raphson or Arc-length integration procedures.In the framework of a dynamic analysis, the algebraic equations are solved by using an implicit time integration scheme based on a variable step-size backward differentiation formula (BDF).A synoptic representation of the numerical procedure as well as the computational algorithm implemented in the FE environmental program are reported in Fig. 3.

RESULTS
n this section, results are developed with the purpose to verify the consistency and the reliability of the proposed model.At first, a layered structured formed by four mathematical layers and three intact interfaces are investigated in the static framework.The main aim of the present analysis is to validate the proposed procedure to predict the onset conditions and crack growth for a case involving multiple debonding mechanisms.Subsequently, in order to validate the procedure to describe the crack front speed, dynamic debonding mechanisms concerning a FRP strengthened steel beam specimen have been investigated by means comparisons with numerical results arising from the literature.G I Figure 3: Schematic representation of the algorithm for layered structure, crack initiation and evolution.

Layered Structure -Multiple debonding mechanisms
The loading scheme, reported in Fig. 4, is based on clamped end conditions and concentrated untisymmetric opening forces.The mechanical properties assumed for the laminate are reported in Tab. 1, whereas those concerning the cohesive interface are reported in Tab. 2. The numerical model is discretized along the thickness by using one mathematical layer for each sublaminate, whereas, for the interfaces, three ALE elements are introduced between the sublayers, in which the crack initiation could be potentially activated.The analysis is developed under a displacement control mode, to ensure a stable crack propagation.In order to verify the stability and accuracy of the solution, several mesh discretizations, ranging from a coarse uniform distribution a refined one, are considered.In particular, for the proposed model, the following numerical cases are analyzed:  uniform discretization with a characteristic element mesh equal ΔD/L=2/200 (M1) with 1841 DOFs;  uniform discretization with a characteristic element mesh equal ΔD/L=1/200 (M2) with 3633 DOFs; In addition, in order to verify the consistency of the proposed approach, a model based on Pure Cohesive approach, namely PC1, is developed, in which a uniform discretization of the mesh with a length equal ΔD/L=0.2/200involving in 12012 DOFs is adopted.In Fig. 5a, results in terms of resistance curve are reported.The loading curve, obtained by the proposed model, is in agreement with the results obtained by using refined CZM approach.Moreover, in the case of a very low mesh element number (M1), the prediction in terms of resistance curve is not affected by a divergent behavior, but it is always very close to enriched one, namely PC1.In Fig. 5b, the evolution between crack tip and applied displacements for two different mesh discretizations are considered.The results show that the proposed model is quite stable, since the predictions in terms of crack tip displacements coincide with that of the PC1 solution.However, it should be noted that in the case of a pure cohesive approach, the crack tip position is taken as the point where the fracture function of the cohesive interface tends to zero, whereas, in the proposed model, an explicit movement of ALE region is identified, since it corresponds to a variable which enters in the computation.

FRP strengthened steel beam specimen
The analyses are developed with reference to loading schemes based on the 4-point bending, in which the dynamic effects are considered from both onset and evolution mechanisms.The loading, the boundary conditions and the geometry are illustrated in in Fig. 6, whereas the mechanical properties assumed for the steel, the adhesive, the FRP strip and those concerning the potential cohesive zone model are reported in Tabs.3-6, respectively.In the present study, comparisons with results arising from the literature [18,19] are developed.The main model refers to a steel beam, strengthened with FRP strip elements.The model is based on two cohesive interface elements, which are introduced between adhesive-steel and adhesive-FRP strip elements.As a consequence, debonding phenomena may affect the layered structures at two different interface levels.The interface law utilized to reproduce the debonding process is consistent with the model proposed by [20].In order to obtain a stable crack propagation, the structure is loaded under a displacement control mode.In particular, to avoid the dynamic effects due to the external load, a very small loading rate equal to 1 mm/s is assumed.However, time steps are modified during the computation from 1E-3 to 1E-7 sec, before and after the activation of the debonding phenomena, to capture accurately the effects produced by crack growth.In Fig. 7, results in terms of resistance curve and crack speed time histories for different thickness of the FRP strips are reported.At first, the structure presents a linear, stable and quasi-static behavior.Subsequently, when the crack growth criterion is satisfied in the adhesive-steel interface, the ALE interface is activated to reproduce the debonding phenomena.
During the activation of debonding mechanisms, the resistance curve presents an oscillatory and variable behavior which varies very fast.In the same figure, a detail of the resistance curve at the point in which the crack onset is activated is also reported.This trend is quite in agreement with similar experimental results available from the literature [21], which show the importance of the dynamic effects during the crack growth.It is worth nothing that the resistance curves are quite dependent from the thickness properties of FRP strip.In particular, an increase of the FRP strip thickness reveals a similar impact on the critical displacement and load at the onset of the dynamic process (Fig. 7).Increasing the thickness of the FRP strip, the edge debonding strength of the beam is reduced (Fig. 7a).This effect is attributed to the increased amount of energy that is accumulated in the stiffened FRP layer and the corresponding increment of the edge stresses.Once the dynamic process is activated, the influence of the FRP strip thickness produces an increase of the crack speeds, which leads to more severe failure mechanisms.Contrarily to the properties of the FRP layer, which are well documented in the literature, the influence of the adhesive on the debonding phenomena is not completely investigated.To this end, in Fig. 8, results in terms of resistance curves and crack speed time histories for different values of the thicknesses of the adhesive layer are presented.In particular, an increment of the adhesive thickness reveals a different impact with respect the previous analyses in terms of FRP strip characteristics, i.e.Fig. 7.As a matter of fact, the results show how by using thin adhesive layers, an increase of the dynamic debonding strength is observed (Fig. 8a) leading the structure to be affected to a more severe dynamic state (Fig. 8b), since the observed crack tip speeds tend to be increased.From the results reported in Fig. 7 and 8, a good agreement with the data available from the literature is also observed [18].
In Fig. 9 the interfacial tractions across the two cohesive interfaces, i.e.Adhesive-Steel (as) and Adhesive-Frp (af), for different time steps of the delamination process, are reported.At first, in Fig. 9(a), the distribution of the interfacial traction forces is presented for the status A of the zoom reported in Fig. 7a and 8a, which basically corresponds to the peak load of the quasi-static branch.It represents the stage just before the initiation of the debonding process, in which all layers are still bonded together.However, at the point A, the non-linear response of the cohesive adhesive-steel interface shows how the interfacial normal and tangential tractions tend to zero.This reflects the initiation of the dynamic debonding failure.In Figs.9b-d, the representation of the evolution of the dynamic debonding, in terms of interfacial traction, has been reported for different lengths of the debondend region.In particular, the results are referred to the points B, C, D of the zoom reported in Figs.7a-8a, in which the debonding lengths of adhesive-steel region are equal to 25, 50 and 75mm, respectively.It is worth noting that the af does not debond but it is able to provide interfacial tractions between the mathematical layers (Fig. 9).Finally, the consistency of the proposed model has been investigated also in terms of computational efforts.In represent the generalized strains,   N,T ,M are the generalized stresses defined as a function of the classical extensional   A , bending   D , bending-extensional coupling   D and the shear stiffness   H variables,

.
For each mode components, the Traction Separation Law (TSL) is assumed to be described by the critical cohesive stresses, ( ) , c c t n T T , the critical and initial opening or transverse relative displacements, namely ( ) 0 , c n n D D and ( ) 0 , c t t D D .It is worth nothing that the proposed model is quite general to include other existing cohesive formulation on a different TSL or stress based initiation criteria, just by modifying the analytical expression defined in Eq. 3.

Figure 2 :
Figure 2: Representation of the coordinate systems employed: Before crack initiation (a), After crack initiation, material and moving coordinates systems are coincident (b), ALE formulation: referential and moving configuration introduced to described debonding phenomena (c).


indicates the displacement of the crack tip front of the k-th debonded interface, i f g is the fracture function defined in Eq,(3) and   /    represents the value of the    variable evaluated at left (-) or right (+) crack tip position.

Figure 5 :
Figure 5: Comparisons in terms of loading curve with pure cohesive approach (a); Comparisons in terms of cracks tip position with pure cohesive approach (b).

Figure 6 :
Figure 6: Steel beam configuration and loading scheme.

Figure 7 :
Figure 7: Comparisons in terms of loading curve for different thickness of the FRP strip (a); Comparisons in terms of time histories of the debonding front speed for different thickness of the FRP strip (b).

Figure 8 :
Figure 8: Comparisons in terms of loading curve for different thicknesses of the adhesive layer (a); Comparisons in terms of time histories of the debonding front speed for different thickness of the adhesive layer (b).

Table 1 :
Geometrical and mechanical properties of the laminate.

Table 2 :
Interface properties of the laminate.

Table 3 :
Geometrical and mechanical properties of the steel beam.

Table 4 :
Geometrical and mechanical properties of the adhesive layer.

Table 5 :
Geometrical and mechanical properties of the Frp strip.

Table 6 :
Interfaces parameters of the Adhesive-Steel interface (as) and Adhesive-Frp interface (af).