Review on the damage behavior of metal laminated composite

Metal laminated composite (MLC) have attracted widespread attention for their potential in achieving excellent properties. However, damage and failure are more likely to occur during the plastic deformation due to the existence of the interface region, leading to serious safety disasters. Thus, it is essential to research the damage mechanism accurately in order to instruct the practical design of the composite structures. This paper provides a comprehensive review of the experimental methods, various failure criteria and the finite element analysis for studying the damage behavior of the MLC. Some advanced failure criteria such as the cohesive theory and the Gurson-Tvergaard-Needleman (GTN) model which are used to predict the damage behavior of interface region and metal matrix of the MLC are discussed. The establishment of exponential cohesive/shear modified GTN mixed model are expected to accelerate the damage mechanism analysis of the MLC.


Introduction
With the rapid development of aerospace, fuel vehicle and machinery industry, the demands of lightweight, good corrosion resistance and excellent mechanical properties for components are increasing. Metal laminated composite (MLC) have attracted widespread attention for their potential in achieving high performance. These excellent properties are achieved by the combination of the ductile metal matrix and the stiff intermetallic compounds layers [1][2][3][4][5].
For instance, Ti/Al MLC is being considered as an important material for structural and military applications because of its high specific strength, relatively good toughness and lower density [6][7][8]. Cold-rolled Cu/Cu MLC shows an enhanced tensile ductility, twice larger than individual cold-rolled Cu sheets without sacrificing strength [9]. Lu et al found that the strain of the MLC occurred in a larger region rather than confined to the necking region [10,11]. Similar conclusion was also reported by Lhuissie et al, who studied the tensile deformation behavior of the high carbon steel/low carbon steel MLC [12]. When the strain reached 20%, there was no crack in the MLC, only showed a large strain concentration. However, damage evolution and failure may be provoked during the complex plastic deformation [13]. Therefore, it is essential to research the damage mechanisms of the MLC.
Under the action of external load, voids and microcracks will nucleate in the internal structure of the components. These defects gradually propagate and eventually form macroscopic fracture [14,15]. This process is called damage, which was firstly proposed by Kachanov et al to study the creep rupture of metals [16]. Armstrong et al found that the macroscopic structural integrity of the components did not indicate their performance was well [17]. When the macroscopic plastic deformation was small, some regions inside the material had reached the limit of plastic deformation. In these places, damage would occur preferentially, leading to the initiation of the microcracks. MLC is composed of the metal matrix layer and the intermetallic compound layer. The MLC with special layered structure is more prone to produce strain localization, leading to the damage during plastic deformation. According to the regions where the damage initiates, the damage behavior of the MLC can be divided into the interface region damage caused by intermetallic compounds fracture and the metal matrix damage caused by metal matrix fracture [18,19], as shown in figure 1(a). The occurrence of damage leads to serious degradation of the mechanical properties of the components, restricting the wide application of the MLC.
Up to now, various experimental methods and failure criteria have been proposed to study the damage evolution of the MLC. Besides, with the development of the computer technology, the finite element analysis has become a powerful tool to deal with the complex nonlinearity problems. In this paper, the experimental methods, various failure criteria and the finite element implementation for studying the damage behavior of the MLC are reviewed. It provides a reliable fundamental for the research of the damage mechanism of the MLC.

Research on the damage behavior of interface region
The interface region of MLC is composed of various intermetallic compounds [20][21][22], as shown in figure 2. Their mechanical properties and crystal structures are different from the metal matrix. The metal matrix layers are bonded together by the intermetallic compounds. When the deformation load exceeds the critical bonding strength, damage will inevitably occur in the interface region. Interface region plays an important role in the MLC during plastic deformation. An interface region with high bonding strength can redistribute the stress and reduce the local stress concentration during the deformation, making the initiation and propagation of cracks more difficult and improving the toughness of the MLC [23,24].
The bonding strength of the interface region is affected by many factors such as the crystal match degree, the voids volume fraction and the interface morphology. Zhang et al studied the microstructure and the deformation behavior of the Cu/X (X=Au, Cr) MLC [25,26]. The property of the interface region was determined by the match degree between the metal matrix and the intermetallic compounds, which played an important role in the plasticity and fracture behavior of the MLC. Pei et al researched the microstructure of the interface region of Ti/Al MLC and its effect on the tensile fracture behavior [27]. The bonding strength of the interface region depended on the crystal match degree, and an interface region with good bonding strength could delay the occurrence of cracks. Khalid et al found that the formation of voids in the interface region was the main factor for the failure of Cu/steel MLC, and the reduction of the voids volume fraction after annealing could delay the damage initiation [28]. Ning et al prepared Zr/Ti/steel MLC with a wave interface region using explosive welding [29]. Shear tests on the Zr/Ti/steel MLC showed that the shear strength of the MLC with a wavy interface region was high because of the mechanical locking effect. Besides, the wave interface morphology can make the stress field distribution more uniform at the crack tip, increasing the fracture toughness of the interface region [30].
To research the role of interface region during plastic deformation, Thiyaneshwaran et al prepared multilayer Ti/Al MLC by solid-state diffusion bonding and performed compression tests on the MLC parallel and perpendicular the interface region [31]. Cracks were only formed in the interface region regardless of the loading direction. Pei et al used in situ tensile test to observe the fracture process of Ti/Al MLC [32]. Cracks were initiated in the interface region, then propagated to the metal matrix in which lead to the fracture of the MLC. Yanagimoto et al investigated the bending and tensile formability of stainless steel/stainless steel MLC [33]. It was found that the ductility of the MLC was enhanced. The necking was inhibited by the interface region and the stress field was relocated after the onset of necking. Huang et al designed multilayered Ti/Al MLC with superior ductility than individual Ti or Al sheets [34]. The stress/strain evolution of the MLC was analyzed by digital image correlation (DIC) technology, as shown in figure 3. It was found that the interface region could transfer strain from Ti layer to adjacent Al layers and relieved the strain localization of Ti layers. Besides, the crack propagation is constrained by the interface region, which improves the ductility of the Ti/Al MLC.

Research on the damage behavior of metal matrix
The damage behavior of the metal matrix layer of the MLC is similar to that of the homogeneous metal materials. Therefore, the research on the damage behavior of homogeneous metal materials can provide a necessary reference for studying that of the metal matrix layer.
The damage behavior of homogeneous metal materials is determined by the stress-strain field and the voids evolution inside the material. According to the scale describing the damage behavior of homogeneous metal materials, it can be divided into the continuous (macroscopic) damage mechanics and the meso-damage mechanics. Continuous damage mechanics breaks through the concept of fracture mechanics and considers the cumulative effect of damage process by defining damage factors [16,35]. Based on the continuous damage mechanics, several damage models have been proposed. Lemaitre proposed a linear law for damage evolution with plastic strain [36]. Later, Bonora proposed a formulation for the damage dissipation potential, which can obtain a general non-linear law to describe ductile damage [37]. Bonora et al validated a model formulation against ductile damage evolution, experimentally measured in A533B low alloy steel under various stress triaxiality conditions [38]. Majzoobi et al used Bonora damage model to study the nonlinear damage behavior of 2024 Al sheet [39]. Results show that the proposed damage model can predict the damage evolution under different stress triaxialities. However, the continuous damage mechanics ignore the meso-structure evolution during the damage process.
Meso-damage mechanics model has obvious advantages. It considers both the stress-strain fields and the meso-structure of the materials. From the perspective of meso-damage mechanics, the damage mechanism of the metal materials can be divided into three stages, as shown in figure 4: (i) voids nucleate near the second phase particles and the inclusions; (ii) the initiated voids grow under the hydrostatic stress and the plastic strain; (iii) the adjacent voids coalesce with each other.
Generally, voids nucleation is related to the second phase particles and the inclusions in the metal matrix. It can be divided into two ways: (i) the second phase particles or inclusions debond from the metal matrix; (ii) the second phase particles or inclusions break themselves. Stress state has a great influence on the void nucleation. Cox et al used round bar type tensile test to research the effect of stress state on the voids nucleation [42]. The results showed that, compared to the smooth round bars specimen (low stress triaxiality), the voids required a lower strain level when it nucleated in the notched round bars specimen (high stress triaxiality).
Most metals have the initial void volume fraction (porosity) [43]. The primary voids and secondary voids will grow during plastic deformation. The voids nucleated by the debonding of the second phase particles and inclusions are ellipsoidal and will be elongated along the principal stress direction during plastic deformation. The voids nucleated by the break of the second phase particles initially appeared as a tiny crack, which gradually opened up during plastic deformation. If the direction of the void is not along the stress direction or the void is subjected to the shear deformation, the volume of the void will not change significantly. The void only rotates or changes its shape. Beremin et al used round bar type tensile test with different notch radius of the specimens to study the effect of stress triaxiality on the voids growth rate [44]. The results showed that the volume of the voids grew significantly at the high stress triaxiality.
The voids in the metal matrix grow, rotate or change shape during the plastic deformation. When the plastic strain accumulates to a certain degree, the adjacent voids will coalesce with each other and form microcracks. Then the material will fracture quickly. Stress state exerts a great influence on the mechanisms of voids coalescence [45]. At high stress triaxiality, voids coalescence is dominated by the voids internal-necking, accompanied by significant growth of voids volume. At low stress triaxiality, the voids volume does not change significantly, but the voids shape changes (elongation, rotation, twisting). And voids coalescence is dominated by the voids shearing. These two mechanisms were verified by Weck et al, who drilled regularly arranged voids on the Cu sheet and used in situ tensile test to observe the voids coalescence behavior [46,47].

Finite element analysis of the damage behavior of MLC
The above literatures mostly used experimental methods to research the damage behavior of the MLC. The advantage of the experimental method is that it can visually observe the damage process during plastic deformation. However, the dynamic evolution of the stress and strain distribution cannot be presented. In recent years, the finite element technique has become an effective tool for the development of advanced precision plastic forming technology. It can complete a large number of experiments on computer, providing a possibility for research the damage behavior of MLC [48].
Generally, there are overall model and layered model to research the damage behavior of the laminated composite. The overall model can improve computation efficiency, while the layered model can estimate the plastic deformation with excellent accuracy. Kashfi et al used an overall damage model to estimate the damage behavior of the fiber metal laminated composite [49]. The model considered fiber metal laminated composite as an overall material rather than a layered material, which significantly reduced the CPU time. In the finite element software, we can use the partition function to divide the composite into the metal matrix region and the interface region, and use different failure criteria to study the damage behavior of different regions. Huang et al proposed a layered constitutive relationship of TA1/Al1060/SS430 laminated composite by considering the metal matrix layer and the interface region, using the constitutive relationship to study the springback angle of the composite during the bending deformation [50]. The result shows that the layered model correlates better with the experiment compared to the overall model. Therefore, in order to interpret the damage mechanism of MLC during plastic deformation, developing an effective finite element model is crucial.

Finite element analysis of the damage behavior of interface region
Based on the finite element technique, the related models for studying the damage behavior of the interface region include virtual crack closure technique (VCCT), extended finite element method (XFEM) and cohesive model. VCCT is the criterion for the strain energy release rate based on the linear elastic fracture mechanics [51]. When the strain energy release rate is bigger than the fracture energy, cracks begins to expand. It is only suitable for predicting the fracture behavior of brittle material. Besides, it requires the pre-set cracks and complex techniques of mesh division near the crack tip, which have convergence problems and only applicable to implicit environments. XFEM can describe the fracture behavior without pre-set cracks [52]. But the energy release rate during the crack propagation cannot be output. It is only suitable for implicit environments and cannot be used well for the simulations of elastoplastic fracture and multi-crack propagation. Based on the cohesive model, the advantages of dealing with the problems of interface damage are: (i) A unified calculation model for crack initiation and propagation is proposed. The traction of cohesive element is a function of crack displacement. It avoids the stress singularity at the crack tip; (ii) Based on the elastic-plastic fracture mechanics, the plastic zone near the crack tip is investigated. It can solve the problem of large-scale yield at the crack tip; (iii) The pre-set crack is not required. (iv) Cohesive model can accurately calculate the crack initiation and propagation process on the crack path. Therefore, the cohesive model is more widely used in predicting the cracking in the interface region of composite. This paper focuses on the research status of cohesive models, as follows.
There are many types of cohesive models such as bilinear, trapezoidal and exponential. It can be selected according to different situations. Table 1 shows the control equations and curves of the three typical cohesive models. Bilinear cohesive model was originally proposed by Mi et al [53]. The traction of the cohesive element increases linearly with the increase of the displacement. After the traction reaches to the maximum value, cracks begin to initial. As the displacement increases, the load bearing capacity decreases, and the cracks gradually propagate. The cracks propagate completely when the traction reduce to zero, and the cohesive element fails at this point. The bilinear cohesive model can better predict the damage behavior of the interface region of brittle materials. But there are errors in describing the large plastic flow of materials under shear deformation, and the sudden transition at the peak point of the traction-displacement curve leads to difficulty in numerical convergence in the finite element analysis [54]. Tvergaard et al proposed trapezoidal cohesive model to study the fracture behavior of the elastoplastic solid material [55]. The model consists of three phases including linear elasticity phase, ideal plasticity phase and linear softening phase. The damage behavior of the interface region is nonlinear [56]. From the mathematical point of view, exponential cohesive model with continuous transition at the peak point of the traction-displacement curve is more suitable for studying the interface damage behavior than the trapezoidal cohesive model [57]. The exponential cohesive model was first proposed by Needleman et al to simulate the crack initiation and propagation until the interface region delamination [58]. Subsequently, Xu et al modified the model and used it extensively for the prediction of interface failure [59]. Chandra et al studied the tangential cracks of interface region of cermet materials by using exponential cohesive models [60]. The results showed that the main parameters and the traction-displacement curve of the model had an important influence on the macroscopic mechanical of the cracking process. Liu et al realized an improved exponential cohesive model using ABAQUS software, and studied the influence of the cohesive shape, cohesive strength and mesh size on the numerical convergence, mesh sensitivity and computational efficiency of carbon fiber composite [61].
The interface damage analysis can be implemented using the famous finite element software ABAQUS. The existing cohesive model in the ABAQUS does not include the exponential cohesive model, so we need to carry out the secondary development of the subroutine, and embed the user-defined mechanical material behavior into the finite element software. User defined mechanical material behavior include UMAT and VUMAT, which are the interfaces of the material property of the ABAQUS. UMAT is applicable to ABAQUS/Standard, while VUMAT can only be called in ABAQUS/Explicit. The UMAT subroutine can complete various nonlinear and complex models, and the subroutine is not difficult. However, when use UMAT to program the cohesive model, there are the following problems. During the calculation process, when the fracture energy reaches the maximum value, the cohesive element completely fails and need to be deleted. But in UMAT, the element does not have the function to delete, it continues to exist when the element fails. Li et al used the UMAT subroutine to program the trapezoidal cohesive model [62]. The crack initiation process can be accurately calculated. But for the crack propagation process at the stage of stress drop, the accurate calculation results cannot be obtained.
By using the exponential cohesive model, the flow chart for the interface damage analysis using finite element method is shown in figure 5. It can be divided into four stages: (i) When the incremental step starts, the program calls VUMAT to transfer the strain increment and the state variables. (ii) The crack displacements of the cohesive elements are calculated from the state variables. The traction values are updated based on the traction-displacement equation. (iii) Calculate the normal and tangential fracture energy and the total fracture energy, and determine whether the fracture energy reach the maximum. (iv) If neither the normal nor the tangential fracture energy reaches the maximum value, the updated state variables is stored for the next incremental step. If the normal fracture energy or the tangential fracture energy reaches a maximum, the cohesive element fails and is deleted.

Finite element analysis of the damage behavior of metal matrix damage
Based on meso-mechanics, the related models for studying the metal matrix damage behavior of MLC are shown below. McClintock and Rice et al first applied voids theory to analyze ductile fracture [63,64]. They studied the growth process of isolated spherical voids in an infinitely ideal rigid plastic matrix and the effect of stress triaxiality on the voids growth, establishing a well-known isolated void growth model: the McClintock model and the Rice-Tracy model. When using the two models to predict the ductile fracture, it is generally considered that the material will fracture when the voids growth rate reaches the critical value, which is regarded as a material constant. However, the interaction between the adjacent voids are not taken into account.
Gurson et al proposed a porous plasticity damage model, which analyzed the growth process of cylindrical and spherical voids and used the upper bound analysis to derive the plastic yield potential considering the influence of the void volume fraction and hydrostatic stress [65]. In the original Gurson model, the metal matrix was idealized as a rigid plastic material obeying the Mises plasticity criterion. The appearance of the voids made the Gurson yield surface closely related to the hydrostatic stress. It changed the characteristics of 'the hydrostatic stress does not affect the plastic yield' and the 'plastic volume is incompressible' in the classical plastic theory. However, the original Gurson model did not contain a failure criterion. The material was weakened by the continuous void growth. The bearing capacity of the materials completely lost when the void volume fraction f reaches 1.
The Gurson model was received great attention after it was proposed. Subsequently, many scholars revised and applied it. Tvergaard et al introduced two parameters q 1 and q 2 to account for the interaction between the voids [66,67]. The void volume fraction f still needs to reach 67% to make the stress of the materials reduce to zero. It is difficult to match the calculation results with the experiment results. Chu et al supplemented the criterion for the voids nucleation [68]. Needleman et al considered the rapidly loss of bearing capacity caused by the voids coalescence [69]. The critical void volume fraction f c was used as a sign of voids coalescence, and it was considered to be a constant only related to the material. Leblond et al improved the predictive accuracy of the model when the strain hardening exponent n greater than 0.2 [70]. Benzerga et al modified Gurson model to make it suitable for the anisotropic metal matrix [71][72][73]. Perrin et al proposed a Gurson model considering the influence of 'secondary voids' [74]. Wen et al considered the effect of void size on the Gurson model [75,76].
Among the many modified forms of the Gurson model, the Gurson-Tvergaard-Needleman (GTN) model takes into account the process of the voids nucleation, growth and coalescence [65][66][67][68][69]. It takes the void volume fraction as the damage variable and can better predict the ductile fracture process [77,78]. However, the GTN model is only suitable for the stress state with high stress triaxiality. The damage evolution of the metal materials is not only determined by the increase of the voids volume, but also by the change of the voids shape [79]. Under the stress state with low stress triaxiality, the void volume fraction will not increase, but the change of the voids shape still causes material damage. Nahshon et al considered the change of the voids shape and established a shear modified GTN model to extend the application of the GTN model to the stress state with low stress triaxiality [80]. The modified Nahshon model includes the shear deformation of the voids. At this time, the damage is no longer indicated as the void volume fraction, but is only expressed a damage variable, which can also be considered as the equivalent void volume fraction. The Gurson models for the various corrections are summarized in table 2.
A reasonable method for obtaining the metal damage parameters is an important guarantee for accurately analyzing the metal damage mechanism. The GTN model needs to determine 3 yield function constants and 4 parameters that related to the voids volume fraction. Generally, all parameters can be obtained by electron microscopy experiment, the representative volume element method or the finite element reverse calibration method [81]. The voids volume fraction of metals can be directly obtained by the electron microscope, but there is often a certain deviation due to the limited observation area. The representative volume element method can obtain GTN parameters by finite element analysis of a single cell. However, the analysis process is complicated and the ideal cell is different from the actual material. The finite element reverse calibration method combines the experimental data and the numerical simulation to obtain the reliable damage parameters. The damage parameters are constantly modified by numerical simulation until it consistent with the experimental data.
The existing GTN model in the ABAQUS software does not consider the shear stress, so we need to carry out the secondary development of the subroutine. The most accurate and reliable numerical integration algorithms for elastoplastic constitutive equations is the backward Euler algorithm [82]. The flow chart for the metal matrix damage analysis using finite element method is shown in figure 6. The calculation process is as follows: (i) Read the initial stress, strain and state variables at the t time: σ t , ε t , a H , t Δε t+Λt , H α (α=1,2 represents two state variables f and ε − p , respectively). If f tr <0, the current state is pure elastic state, then go to step (v) to update variables; If f tr 0, the current state is elastic-plastic state. Go to step (iv) and start the plastic correction.

Literature Description
Gurson [65] The Gurson yield potential for cylindrical voids: The Gurson yield potential for spherical voids: Where q and p are the equivalent stress and the hydrostatic stress, respectively; s 0 is initial yield strength of the metal matrix; f is the void volume fraction. Tvergaard [66,67] Introduced two parameters q 1 , q 1 to account for the interaction between the voids Where s m is the flow stress of the metal matrix; q 1 and q 2 are the introduced constants. Chu and Needleman [68] Added criteria for the void nucleation.
The criteria for void nucleation by strain control: Where e m p  is the equivalent plastic strain of the metal matrix; f N is the limit volume fraction of the void nucleation; ε N is the average equivalent plastic strain at the time of the void nucleation; s N is the standard deviation of the normal distribution; σ m is the average stress; σ N is the average equivalent stress at the time of void nucleation. Tvergaaerd and Needleman [69] Used f * function instead of f in the Gurson yield potential: Where f c is the critical void volume fraction; f F is the void volume fraction at the time of the material failure ; = f q 1 u 1 * is the value of f * when the stress reduce to zero. Leblond [70] Introduced two matrix yield strengths s m1 and s m2 to improve the predictive accuracy of the original model when the strain hardening exponent n is greater than 0.2.  ( ) h i is the anisotropy coefficient.
Nurson [80] Considered the change of the voids shape and established a shear modified GTN model:  Where I is the second order unit tensor; ε p is the macroscopic plastic strain tensor; k s is the shear coefficient; S is the deviatoric stress tensor; σ eq is the macroscopic equivalent stress; J 3 is the third invariant of stress deviator.
(iv) In terms of the associated flow law and yield condition, the plastic correction can be obtained: Where e D s m and e D s eq are plastic strain increment induced by hydrostatic stress and equivalent stress, respectively. The Newton-Raphson iteration is used to solve the e D s m and e D s eq until the tolerance of | m1 | and | m2 |10 −6 .
(v) Update the stress, strain, and state variables at t+Δt time: (vi) Enter the next incremental step.

Conclusions
This paper gives a comprehensive review on the recent developments of the damage analysis of the MLC. The experimental method, various damage criteria and the finite element implementation for studying the interface damage behavior and the metal matrix damage behavior are discussed. Some advanced experimental method such as the DIC technology provide a way to directly visualize the strain distribution and the damage behavior of MLC during plastic deformation. It should be noted that the finite element analysis for damage evolution is an important topic, and proper selection for finite element model are crucial for accurately studying the damage mechanism of the MLC. The cohesive model and the GTN model can better describe damage behavior of the interface region and the metal matrix of MLC, respectively. Therefore, the establishment of exponential cohesive model/shear modified GTN mixed model are expected to accelerate the damage mechanism analysis of the MLC.