Fracture mechanics based assessment of manufacturing defects laying at the edge of CFRP-metal bondlines

This study presents a simplified approach for the evaluation of the bond line load capacity of CFRP-metal hybrid structural parts, bearing manufacturing defects. The proposed approach overcomes the need of modeling debonding with the use of advanced non-linear FEM, and relies on the post-processing of the stress and strain fields developed at the free edges of the hybrid material (critical locations in the presence of manufacturing defects), with the use of interfacial fracture mechanics. First, the mathematics, related to the analytical calculation of the bi-material fracture toughness, is presented and applied to a simple CFRP-metal geometry. Corresponding FE simulations are performed and the J-integral is employed for the numerical calculation of the SERR magnitude of different defect sizes. In the following, a simplified calculation of the SERR magnitude that sources from FE results and implemented in geometries modeled by shell finite elements that typically require advanced modeling techniques, is provided for the simulation of the de-bonding process.


Introduction
The global demand for the minimization of the environmental impact on automotive and aerospace vehicles continuously pose requirements for a lightweight design of the associated structural assemblies. The geometry of components conventionally produced by metallic materials has already achieved a high degree of optimization, in terms of weight savings, even with the introduction of modern stronger alloys. A solution to achieving higher mass reduction is through the replacement of metals with carbon fiber reinforced polymer (CFRP) composites that have superior strength-to-weight and stiffness-to-weight mechanical properties. Typical examples are those found in the aerospace industry with Boeing 787 [1] and Airbus A350 XWB [2] having 43 and 52% of the metallic structure constructed by composite materials, respectively. In the automotive sector, the BMW 7 series has replaced several metallic components with composite materials [3], as shown in Fig. 1. However, there are certainly components whose imposed design requirements are not met if the metal is fully replaced by composites. A good example is that of the B-Pillar component (side part of the Body In White for protecting the passengers in the case of a side collision), which if entirely made of composites, then its crashworthiness is not equivalent to the respective metallic one. A promising solution is by using a hybrid CFRP/metal combination (see Fig. 1), where the CFRP layers are adhesively bonded onto selected parts of the metallic component.
Besides research and development carried out in production technologies that will enable the automated manufacturing of CFRP/metal hybrid parts, it is of high importance that specialized design methods and tools be developed. Recently, the authors have published relevant work [4] on the design of hybrid parts by employing optimization techniques for the calculation of the design variables (number of layers, layer orientation, reduced metal thickness), where in the proposed formulation it is assumed that the bond area remains intact and only failure in the metal and composite has been considered. This study focuses on developing a simplified yet effective modeling technique for the calculation of any failure propagation in the CFRP/metal bond line 1 . Failure analysis of bond lines, formed in adhesively bonded joints with dissimilar materials, has received the interest of several researchers, particularly during the last decade. However, there is still ongoing work on the field, since there do exist complicated failure mechanisms 2 that are not yet fully understood and assessed. The loading of the joint substrates leads to the development of peeling and shearing stresses in the adhesive bondline. Analytical models, sourcing from closed form solutions, have already been developed [5,6]. In these models, the peeling and shearing stresses (or/and strains in some cases) are correlated with the externally applied loads, and hence, when the former are compared with material limits, the design loads may be calculated. Such a strength analysis approach is defined as a stress-based failure criterion and applies well to the prediction of debonding initiation [7]. Alternative load bearing assessment methods are based on the application of fracture mechanics, through the use of energy-based fracture criteria, where debonding propagates once the work of fracture is equal to the strain energy release rate [8,9]. The validity of using stress or energy based criteria depends on the joint's geometry (substrates' thickness and overlap length) and the adhesive material type (brittle or ductile). Since analytical models that are either based on stress or energetic failure criteria are restricted to simple geometries (single and double lap/strap joints), strength analysis of the bondlines, involved in complicated geometries, is performed through finite element (FE) modeling techniques [10,11]. However, detailed FE modeling requires the development of highly dense 3D mesh of the adhesive layer in order to accurately predict the developed stress field for the application of stress criteria or energy criteria. The solution of these models is time consuming and hence not preferred to be employed in the design phase of industrial complicated geometries, where the geometric parameters are continuously modified until the optimal solution has been achieved. On the other hand, cohesive zone modeling techniques that combine stress and energy criteria and can accurately calculate the bondline failure initiation and propagation loads [7,12,13], are highly non-linear in nature and therefore, the computational time required for complicated geometries once again limits their efficient use in an industrial environment. The approach proposed in this study requires a linear static solution of the modeled geometry without having to include material non-linearities that come from the cohesive zone modeling analysis of bonded connections.
The detailed debonding analysis of adhesively bonded CFRP/metal material systems, requires the employment of sophisticated approaches, based on stress and fracture mechanics, due to the different failure modes involved in the problem, as schematically shown in Fig. 2a. This study aims at providing a simple yet effective approach for the calculation of the maximum loads that can be carried by a CFRP/ metal hybrid system, until the manufacturing defects 3 at the bond line, propagate and the two constituent materials separate. The proposed approach considers the problem shown in Fig. 2b, where the diverse failure mechanisms (Fig. 2a) 1 Bondline is defined as the interlayer between the CFRP and metal and is responsible for holding together the two substrates either through an adhesive medium or through the resin matrix material of the composite laminate. 2 The behavior of the adhesive material when used as thin layer to form a bondline in relation to the strength of the metal/adhesive and of the adhesive/composite interface and the way these interact are considered as complicated failure mechanisms. takes place and effective macro-properties may be used as failure parameters.
The equivalent bond line refers to a bond line with effective material properties (max stress and fracture toughness) that from a macro-scopic point of view, involves one or several failure mechanisms that lead to fracture. In a fracture mechanics experimental characterization test, one usually obtains the effective properties of a bond line. In the study performed herein, the thickness has not been taken into consideration since the bond line is represented by an interface (with zero thickness) with effective fracture properties.
The defect in the equivalent bond line approach may be one or a combination of the following: a debonding between CFRP and adhesive, a debonding between steel and adhesive and a failure within the adhesive. In that meso-and microscale it is very difficult for the exact defect mode combination to be distinguished, even with very precise measuring equipment. Therefore, it is convenient from the engineering mechanics point of view that an equivalent macro-debonding be introduced for the assessment of the joint's integrity, based on the length of that defect.

Debonding modes
Apart from the quality of the raw materials, the joining processes involved in the production of hybrid parts, are responsible for the quality, reliability and integrity of the formed bondline, between the two substrates. Since the two substrates (CFRP and metal) are attached together through adhesion and mechanical interlocking mechanisms, any sort of contamination on the substrate surfaces or any process related divergence from the optimal settings, has high potential to lead to poor adhesion quality in the micro-and meso-scale. Porosity, micro-cracks, un-cured adhesive, lack of adhesion and other manufacturing defects in the bondline, may accumulate and act as crack starters. The latter will eventually result in debonding propagation in the macroscale, once the component has been exposed to in-service conditions. The ability to carry further loads is hence minimized and the overall integrity of the hybrid part is jeopardized. Although manufacturing defects may be existent anywhere at the bondline, the ones being critical for propagation are adjacent to the edges, as shown in Fig. 3. A defect lying far from the joint edges cannot propagate (no stress concentration in the presence of in-plane stresses at the substrates), unless one of the substrates has buckled and hence the developed stresses (strain energy) can be released.
On the other hand, the substrates at the edges can deform and release strain energy in the case of an edge effect being prone to propagation. In this study, only edge debonding will be studied, since buckling debonding may be avoided if the substrates alone are designed against buckling. It must be noted that a defect far from the edges will not propagate in the case of in-plane stresses, unless buckling occurs. However, in the case of out-of-plane loading, localized peel stresses might cause the propagation of a defect. In this study though, focus is placed on hybrid materials that develop inplane stresses and hence edge defects are of interest.

Debonding failure criterion
Within the framework of linear elastic fracture mechanics (LEFM), defect propagation can be tackled with methods Fig. 2 The failure modes of a CFRP/metal hybrid system (a) and the proposed analytical and numerical approach for calculating the macro-debonding of the equivalent bondline (b) based on the energy approach to failure. When external loading is applied to the hybrid component, containing an initial defect of size α 0 , then the involved substrates deform, develop stresses and therefore store a certain amount of strain energy, U (units in Joules). The stresses are being transferred from one substrate to the other through the bondline, which load bearing capacity is directly related to the work of fracture or fracture toughness G c (units in J/m 2 or N/ mm). This magnitude is experimentally measured, through specialized fracture tests (three-or four-point bending, double cantilever beam, mixed-mode tests, etc) in laboratory conditions. According to the energetic approach to failure, a manufacturing defect in the bondline will not propagate until the strain energy release rate (SERR), G, given by is equal to or becomes higher than the macro-scale fracture toughness of the bondline G c . Hence, the following fracture criterion can be used for the load bearing assessment of bondlines: Experimental techniques have not yet been standardized for gauging the fracture toughness of bondlines, formed in multi-material systems, despite this field undergoing important research. Tests such as the end notch flexure (ENF), the single leg bending (SLB) and double cantilever beam loaded with uneven bending moments (DCB-UBM) have been used in research studies to gauge the Mode II fracture toughness of bondlines, formed between composites and metals [14].
The key advantage of this criterion (Eq. 2) is that it does not consider the complex stress field in the vicinity of the crack tip but only the energy stored in the material by the remote elastic stresses. Therefore, it is practical to be used (2) G ⩾ G c in the design phase of complicated geometries, as it will be shown in later sections.

Debonding of a hybrid beam
This section involves the load bearing assessment of a bondline formed in a CFRP/metal 1D beam geometry that is subjected to in-service conditions and contains an edge manufacturing defect of length α 0 . An infinitesimal element dx, in the vicinity of the edge, is taken for analysis, as shown in Fig. 4a. Assuming that the initial crack length is very small compared to the length dx, then the developed stress field is uniform along the entire length of the element and its magnitude and through-the-thickness distribution do not allow the development of high enough strain energy for the defect to propagate (G < G c ). Magnification of the stress profile will enable the defect to propagate (G ≥ G c ) and hence the debonding will interrupt the load transfer from one substrate to the other, since the part of the composite substrate, behind the crack tip, cannot absorb further load (Fig. 4b).
Considering the composite layers as one isotropic layer with equivalent elastic constants, the problem shown in Fig. 4c is solved in [15] and the obtained formula for the calculation of the SERR magnitude of the bi-layer is provided below: This formula yields G ss , which refers to the constant rate of energy being released, during the steady state debonding propagation. P is the resultant axial force and M is the resultant bending moment acting on the metal substrate, as shown in Fig. 4c. This formula physically expresses the difference in the strain energy density (normalized strain energy to the It 3 c Fig. 3 Possible modes responsible for debonding propagation area of the crack faces) calculated at the beams behind and in front of the defect tip. The composite part behind the crack tip has debonded and therefore, it does not carry any stresses, and consequently its strain energy is zero. M b is the balance moment acting at the composite/metal bi-material at the location of the neutral axis and it is given by the following formula: The location of the neutral axis δ is calculated by the following set of formulas: Magnitudes Ē c and Ē m correspond to the effective plane stress Young modulus of the composite and metal layers, respectively, and are given by: where, E c , E m and v c , v m are the substrate Young moduli and Poisson's ratio, respectively.
The normalized axial stiffness A and bending stiffness I are given by the following formulas: Equation (3), includes only geometric and material parameters alongside with the loading magnitudes, excluding the defect size α 0 . This feature is advantageous for the assessment of the bondline's load capacity within the design process of hybrid parts, since the exact location, density and size of the manufacturing defects are unknown parameters. Additionally, the provided formula considers only the thickness of the associated substrates and not any length related magnitudes, which in turn, result in the steady-state energy release rate G ss , which refers to the constant rate of energy being released during the steady state debonding propagation. As it will be shown in the next section, there does exist a certain defect size beyond which the SERR magnitude becomes constant, G ss .

Effect of defect's length
In order to evaluate the validity of Eq. 3 in hybrid parts with defects of finite length, numerical simulations based on finite element methods have been performed. The simple hybrid geometry shown in Fig. 5 has been considered in order for the effect of the defect's length α 0 on the SERR magnitude to be evaluated. The length of the metal l m is equal to 200 mm and the length of the composite patch l c is equal to 100 mm. The metal thickness t m is 2 mm. The composite consists of eight layers with a stacking sequence [0/90/45/-45] s and with a ply thickness t l equal to 0.2 mm.
The total composite thickness t c is equal to 8 × 0.2 = 1.6 mm. The width of the hybrid system (normal to the cross sectional view of Fig. 5) is taken as equal to 100 mm. In order for Eq. 3 to be applied to the considered geometry that is formulated according to the 1D beam theory, the effective stiffness of the composite layers has been calculated according to the classical theory of laminated composites. Hence, the effective Young modulus of the composite stack in x direction is equal to 73.9 GPa, given the ply elastic constants, i.e., E 1 = 178 GPa, E 2 = 10 GPa, G 12 = 4.55 GPa, v 12 = 0.3 and v 21 = 0.02. The metal substrate was modeled with a Young modulus of 210 GPa and a Poisson's ratio of 0.3. Although only an axial load P of 50 kN has been applied to the metal free edges, the hybrid area (CFRP/metal overlap) develops secondary bending stresses on top of the axial ones. This effect is due to the fact that the load application point on the metal (t m /2) is different from the location of the neutral axis of the hybrid part (δ), as it is schematically shown in Fig. 4c. The geometry has been modeled in ANSYS commercial FE software in both a 2D and a 3D space with linear quadrilateral finite elements, as shown in Fig. 6. The left end of the metal has been considered as being simply supported (constrained displacement at x and y direction, rotation is permitted) and the load is applied to the right end of the metal. Five models have been generated and solved as linear static by varying the length of the defect (α 0 = 0.2, 0.3, 0.4, 0.5, 1 and 3 mm) and keeping the applied load equal for all cases. A quite dense mesh has been used in all models (2D and 3D) in the vicinity of the defect tip in order for the complex stress field to be accurately captured. The SERR magnitude has been numerically calculated through the use of the J-integral approach, which is a pre-coded feature in ANSYS program. Figure 7 presents the effect of the defect's size, lying at the edge of the bondline, with respect to the SERR magnitude. The numerically calculated G values are normalized by the G ss magnitude as calculated in Eq. 3 and are equal to 0.084 N/mm (84 J/m 2 ). It is clear that SERR is increasing together with the defect's size, in a progressive (non-linear) manner, until it reaches a plateau level, where it remains constant and it is not further influenced by longer defects. This plateau corresponds to the steady state SERR condition. For a given load case, Fig. 5 Simple hybrid geometry considered for load bearing assessment debonding propagation will occur only under the existence of a defect, with a size equal to or greater than a critical size, where in the considered example it is equal to approximately 1 mm. Hence, defects of smaller size will not propagate, since there is not enough strain energy stored into the substrates (< G ss ) to be released, according to the energetic criterion, described by Eq. 2. The existence of a relatively minor defect at the edge will not cause any noticeable stress redistribution that would allow energy release. As the defect's size increases, the stresses redistribute and part of the edge material, being on the upper side behind the crack front gets unstressed, in a way that strain energy may be released. The analytical solution does not model this transition and considers only a fully unstressed material behind the defect front and hence, it provides the steady state SERR magnitude. On the other hand, a detailed numerical model is able to capture this transition in the stress relief process towards the steady state. The analytical method for evaluation G ss is also applicable to other material/thickness combinations, as long as the problem remains one dimensional, where axial loading and bending moment are applied in one direction. In case of in-plane load combinations, then the analytical formula has to be re-worked. On the other hand, numerical simulations, based on FEA, may be used in complex loading conditions by varying the size of a defect, until the plateau level has been found.
Similar conclusions can be drawn from the parametric analysis in the 3D space, where the effect of the strain energy release rate magnitude, along the width of the joint, was calculated, as shown in Fig. 8. The SERR magnitude converges on a constant value that is 8% higher than the theoretical steady state G ss magnitude, for a critical defect size approximately equal to 1 mm. This difference is attributed to the plane stress formulation of Eq. 3, which ignores the arising effects in 3D geometries, as the one studied here. Given that plane strain conditions are assumed in Eq. 3, then, the corresponding theoretical results meet the numerically predicted ones with good accuracy. The distribution of the SERR magnitude along the width is approximately constant when moving from the edges to the center. The predicted peaks at the edges, result from stress concentrations there, well known as edge effects and debonding is prone to propagating from these locations.

Simplified debonding analysis
In this chapter, a simplified yet effective method that can be used in the design of hybrid CFRP/metal parts, for load bearing assessment of the involved bondline(s) is developed. As it has already been mentioned, debonding will initiate and propagate from the part's edges, given the existence of a manufacturing defect, which acts as a crack starter. The location and the size of such defects are stochastic and unknown during the design phase of the hybrid part. Therefore, the engineer should analyze several scenarios with different distributions of defects, based on variable sizes and locations, for the effective evaluation of the bondlines' load capacity under the existence of edge manufacturing defects. The FE modeling of hybrid parts is usually performed on surface based geometries with the use of shell type element formulations, because the thickness of the associated substrates is very or moderately small compared withthe parts' plane dimensions. However, load bearing assessment of bondlines, based on the FE techniques as the ones used in the previous chapter, requires specialized detailed 3D modeling (solid elements) that includes the size and location of defect(s). Considering that the FE mesh used for the accurate calculation of the SERR magnitude in the simple patched problem (see Fig. 5), involves approximately 7000 nodes for the 2D plane case and approximately 100,000 nodes for the 3D solid case (see Fig. 6), then the construction and solution of different defect scenarios, based on such FE models In order to overcome these issues, an efficient method is proposed below for the calculation of the SERR magnitude in hybrid parts of complex geometries.

Method description
The proposed method takes advantage of the fact that the driving force of a bondline fracture, i.e., the steady state strain energy release rate, does not depend on the defect's size for values equal to or higher than a critical size level. It has already been shown that for a given load case, defects of a smaller size than the critical ones will not propagate, unless the loads are magnified or unexpected in-service loads occur. The goal is that magnitude G ss be calculated at the material points, located around the bondline's edges of a hybrid part, modeled in FE environment with the use of shell elements. In such modeling techniques, it is assumed that the material develops in-plane stresses, as illustrated in Fig. 4a, whereas the out-ofplane stresses remain at low levels, hence are ignored. In this regard, the strain energy may be normalized by the element's area, which leads to the strain energy density Ū with units of energy/area (J/m 2 ) that compares with the units of fracture toughness magnitude. The difference between the strain energy density ahead of the defect tip and behind the crack tip yields the steady state SERR magnitude as given by the following formula: The strain energy density is given by the formula: where, σ and ε are the in-plane stress and strain tensors in a Cartesian frame on the surface of an element, Due to tensor's symmetry, σ 12 is equal to σ 21 and ε 12 is equal to ε 21 . Since the post-processing results of the finite element solution provide the stress and strain tensors for each defined layer of each element, then by approximating the integral of Eq. 12 and substituting with Eq. 13, the following equation is generated: where the summation runs through the total number n of the associated shell element layers.
Stress and strain components are calculated at the integration points of the finite element mesh. Quadratic shell elements, with one in-plane integration point, located at the center of the element (following Gauss numerical integration), are well known for their accurate prediction of the developed stress field. The through thickness stress variation is considered as linear and calculated through trapezoidal integration rules.

Method validation and benchmarking
The above mentioned approach, which is modeled with quadratic (8-node) shell elements, has been applied to the geometry presented in Fig. 5. A parametric analysis has been conducted by varying the in-plane dimensions of the shell elements attached to the bondline edge, i.e., 0.5 × 0.5, 1 × 1 and 5 × 5 m 2 . Following the solution of the FE models, Eq. 11 was applied to the stress-strain results at the row of (13) = elements ahead and behind the bondline edge via Eq. 14. Figure 9 presents the distribution of the SERR magnitude to the modeled geometry with element dimensions 1 × 1 mm 2 . It can be seen that the maximum calculated SERR value is 0.087 N/mm that differs only by 3.4% from the analytical steady state magnitude G ss obtained by Eq. 3, i.e., 0.084 N/ mm. Figure 10 presents the results of the parametric analysis, in terms of the distribution of the SERR magnitude along the width of the hybrid system. It is evident that there is no important influence of the element dimensions on the SERR magnitude and the results show that the proposed method calculates with great accuracy the analytical steady state energy release rate magnitude. The extent of defect size comprising the shell element model is directly related to the critical defect size and the shell element size. Hence, the designer has to be aware of the critical defect size in order to mesh the surface geometry with shell elements of accurate dimensions. Two dimensional FE models with the given material and thickness combinations may be used prior to shell modeling.

Conclusions
A simplified method has been proposed for the evaluation of the criticality of manufacturing defects, lying on the bondline, on the structural integrity of the hybrid composite/ metal systems. A simple patched geometry that enables the application of an analytical equation for the calculation of the energy release rate, has been selected for benchmarking purposes. Finite element simulations have shown that there does exist a critical defect size that will propagate once the SERR magnitude is equal to the analytically calculated steady-state release rate. Smaller defects than this critical size will not propagate under the design load case. The proposed method is a post calculation tool that when applied to FE simulation results, the driving force for defect propagation, i.e., the SERR magnitude can be calculated with high accuracy. The validity of the proposed method is yet to be evaluated in complex geometries under the existence of complicated stress fields and diverse load cases.