Peeling of metal foil from a compliant substrate

ABSTRACT Large displacement peel was studied for cases where a compliant substrate leads to a large value of the root rotation. An existing simplified beam model to calculate the peel fracture energy was modified to allow for a kinematic hardening beam model of the foil. The steady-state peel force and the root rotation were used as input data to the resulting analytical beam model. Test results from the literature were analysed. A more elaborate finite element model was also studied, using cohesive elements for the interface layer between the foil and the substrate. The cohesive zone parameters used were the fracture energy, the cohesive strength and a shape parameter. An optimization scheme for the cohesive zone parameters was developed and optimized against experimental steady-state peel force and root rotation. The optimization scheme was effective to characterize the cohesive parameters. The method yields similar values of fracture energy for the two peel angles, with the one for being slightly higher than for . The difference in fracture energies for different peel angles suggests that the fracture energy can be mode dependent.


Introduction
During the last decade finite element (FE) simulation has been more widely used in packaging industry to facilitate development in the early design phase. FE models are producing accurate prediction at application level in different length scales. [1] Specifically, thin flexible laminates which are important in liquid food industry are experimentally and numerically studied by many authors. [2][3][4][5][6][7] To begin with, it is necessary to know material properties for constituents of laminate. Additionally, the cohesive properties of the interface between layers are important and it was the focus of the present study. Specifically, adhesion between an aluminium foil (Alfoil) and low-density polyethylene (LDPE) film is studied. Several existing studies addressed the interface of similar material laminates. [2,8,9] The peel test is a simple setup, which makes it one of the most widely studied test methods used to assess the adhesion between two layers of a thin flexible laminate. Despite the simple test geometry, the evaluation of an experiment in terms of the fracture energy is not always as easy. The task does not become easier if a complete set of cohesive properties is requested. In an early contribution, Rivlin [10] presented a method to determine the fracture energy from a steady-state peel test. The method was based on a model where the peel arm was assumed to be inextensible and to have a small resistance to bending. For this case, the fracture energy can be determined by measurement of the peel force and the peel angle. This work was later [11] extended to allow for a linearly elastic peel arm. With this model, also the elastic modulus is needed to determine the fracture energy. Later works have considered more elaborate models for an extensible peel arm having a small resistance to bending, cf. e.g. [12,13] From knowledge of the peel angle, the peel force and the stress-strain curve of the peel arm material, the fracture energy is determined in a straightforward manner. For cases where the bending resistance is not negligible and the material is elastic-plastic, one needs to consider also the plastic energy dissipated due to elastic-plastic loading, unloading and reversed plastic loading. Moreover, one also needs to consider the equilibrium equations in order to know the deformation field from which the plastic energy dissipation is determined. A method to achieve this for an elastic perfectly plastic material was developed in Aravas et al. [14] The method was further developed in Kinloch et al. [2] by modelling the material in the peeling arm according to a linear isotropic hardening rule. In the present study, a similar solution for linear kinematic hardening is presented. This is believed to be a more accurate model for a peeling arm in the form of a metal foil. Many other authors [15][16][17][18][19][20] have studied plastic energy dissipated in the peel arm for different thin laminates.
The modelling of the interaction between the peel arm and the substrate is often performed using cohesive zone models. Analytical solutions for peel are presented in e.g. Lu et al. [21] and Begley et al. [22] Through development of cohesive elements in FE-codes, the cohesive modelling technique has become more versatile. The peeling behaviour for the case of an elastic peel arm was studied in Sauer [23] using a van der Waals adhesion model for the cohesive behaviour. Other work [19,20,24,25] studied plastically deforming peel arms. Use of a bi-linear cohesive law was the most popular choice. A bilinear model has effectively two parameters; the fracture energy and the cohesive strength. In some studies [19,24,25] , the cohesive strength is arbitrarily chosen. In the current study, a trapezoidal cohesive law is assumed, by which an additional shape factor is introduced. FE-simulations show that, for the present case, the peeling behaviour is governed not only by the fracture energy, but also by the cohesive strength and the shape factor.
Although FE-simulations of delamination using cohesive zone modelling is fairly standard, the experimental determination of cohesive laws is certainly more challenging. For comparatively stiff adherends, test geometries allowing for purely elastic deformation can be designed. This makes it possible to measure the energy release rate (ERR) as a function of cohesive deformations. Cohesive relations are then obtained by differentiation of experimentally obtained ERR with respect to measured cohesive deformations, cf. e.g. Andersson and Stigh [26] and Leffler et al. [27] The method also works for plastically deforming adherends, provided that the plastic loading is monotonically increasing, cf. Stigh et al. [28] In Nilsson [29] , a similar method is proposed for 90 peel. However, the method requires the peel arm to be elastic. Thus, the method is not applicable for peeling of thin laminates, as cohesive separation of the layers is often associated with plastic loading, elastic unloading and reversed plastic loading of the layers.
An optimization scheme can be useful to identify cohesive parameters with some confidence. In Lee et al., [30] the authors used design of experiment and the kriging metamodel to optimize cohesive strength by comparing peel simulation and experimental force response at few key points. Cohesive parameter set of fracture energy and cohesive strength was optimized in Xu et al. [31] using a genetic algorithm (GA). In the current study, an artificial neural network (ANN) inside a GA was successfully adopted for the optimization, similar to earlier studies in different contexts. [32,33] The present article is organized as follows: The two different models used are briefly described in Section 2. A beam theory solution, for a linear kinematic hardening peel arm material model, is presented in Section 3. The finite element model and the optimization procedure to determine cohesive properties is outlined in Section 4. Experimental observations from the literature are recapitulated in Section 5. In Section 6, the experiments are evaluated in terms of fracture energy and other cohesive properties using the methods described in Sections 3 and 4. The results are discussed and the paper ends with some conclusions.

Models of laminate
The materials tested are the constituent of liquid food packaging with multiple layer of films and foils ( Figure 1). The most widely used packaging materials in liquid food packaging industries are LDPE, Al-foil and paper board. Each of the layers of a package has its own role. For example, the paperboard, which is considerably thicker than the other layers, bears the load when the package is filled, folded and gripped, while the Al-foil isolates the liquid inside from light and diffusion. [8] The outer decor layer of LDPE protects the paper and print on it from moisture-related damage. It is important that the inside layer which is in contact with the product inside does not react or dissolve with it and contaminate the product during its expected lifetime. LDPE layers are used to serve this purpose. The layers are combined together in several steps. The Al-foil is produced by thinning Al-sheets by rolling them several times. Also, during the manufacturing process, LDPE is extruded and rolled with other layers. The laminate of LDPE film and Al-foil are most important among the layers. The adhesion between them dictates the load carrying capacity and allowable strain in the laminate. This makes the study of adhesion between them very important in packaging applications and many previous works have addressed this phenomenon. [3,4,6,7,[34][35][36][37] Firstly, a simplified analytical beam theory model of this multi-layered packaging material was considered for the peel study in this work. This model is depicted in the left lower part of Figure 1. The decoration layer, paperboard layer and laminate layer are modelled as a rigid base. The aluminium layer, inside adhesion layer and inside layer are considered as a single layer peel arm ( Figure 1). The peel arm was assumed to be represented by only stand-alone aluminium layer. It was motivated by the assumption that the LDPE part in the peel arm does not play an important role in the peeling process due to its low Young's modulus and higher ductility. The deformation of the substrate is assumed to be concentrated to the interface between the Al-foil and the laminate LDPE-layer.
Secondly, a corresponding FE-model was considered, cf. the lower right part of Figure 1. The Al-foil was modelled with beam elements and cohesive elements modelled the interface between the Al-foil and the laminate LDPE-layer. The top face of the cohesive elements was connected to the beam elements of Al peel arm by beam connectors. This is needed since the nodes of the beam elements are located at the centroid of the beam.
The beam theory solution presented in Section 3 and the FE-model described in Section 4.1 were both used to evaluate the experimental results given in Section 5. The beam theory solution was used first, to get an estimate of the fracture energy. These results were then used as a starting point in the calibration of the cohesive parameters of the FE-model. The results of this procedure are presented in Section 6.

Beam theory solution
In Aravas et al., [14] a method to evaluate the fracture energy from a peel test was proposed. The method was based on large displacements and the assumption that the peel arm deformed in bending as a perfectly plastic material. The method was further developed in Kinloch et al. [2] by modelling the material in the peeling arm according to a linear isotropic hardening rule. Here, a similar solution for linear kinematic hardening is presented, cf. Figure 2b. This is believed to be a more accurate model for a peeling arm in the form of a metal foil.
As in Aravas et al. [14] and Kinloch et al., [2] it is here assumed that plasticity due to bending of the peeling arm is dominant, i.e. plasticity due to stretching of the peeling arm is not included. This assumption can be inappropriate for small values of the peel angle, θ. However, for the experimental settings studied in the present work, the assumption is appropriate.
Steady-state crack propagation is considered. The fracture energy, G c , is formed as where G e is the energy release rate 1 (ERR) for the present loading configuration, if all points of the system are subjected to monotonically increasing loading. As shown in Kinloch et al., [2] elastic unloading and reverse plastic loading can take place in the peeling arm. Thus, in such cases, the plastic energy dissipation rate, G p , must be subtracted from G e to achieve the fracture energy.
In Appendix A, the J-integral is used to derive G e . With the assumption of small elastic strain at the loading point of the peeling arm where b is the out-of-plane width. Figure 3a depicts the deformation of the peeling arm under steady-state crack propagation. Consider a point, O, located far to the right of the crack tip. At this point the foil is undeformed, i.e. both the bending moment and the curvature are zero. As we move to the left of point O, elastic bending of the foil is experienced and a cohesive zone is formed. At some point, A, within the cohesive zone, the bending stress in the foil has reached the elastic limit, σ y , cf. Figure 2b. As we move further to the left, the foil is plastically deformed and the bending moment increases. Following the methodology in Kinloch et al., [2] the maximum bending moment is assumed to be reached at the current position of the crack tip, denoted B in Figure 3a.
To the left of the crack tip, the foil is unloaded from a plastic state. Between points B and C, the unloading is elastic. To the left of point C, the unloading involves reverse plasticity. At point D, located far from the crack tip, the peeling arm is straightened by the tensile peel force. This means that the curvature is zero at point D. The points O, A, B, C and D are indicated in Figure 3b, where the bending moment, M, is shown as a function of the curvature, 1=R, where R is the radius of curvature due to bending. During steady-state crack propagation, a material point in the foil will experience the entire history of events displayed in Figure 3b. This means that the energy dissipated due to plastic deformation is proportional to the area A OABCD in Figure 3b. As described in Kinloch et al., [2] this means that the plastic energy dissipation rate is given by, The specific form of the moment-curvature plot ( Figure 3b) for a linear kinematic hardening foil is derived in Appendix B. In order to determine the plastic energy dissipation, G p , the conditions at the crack tip, B, must be known.
Here, a modified version of the solution in Kinloch et al. [2] is proposed. The solution is based on large displacement beam theory for the deformation of the peeling arm. As shown in Kinloch et al., [2] the bending deformation of the peeling arm (from point B to point D in Figure 3a) is governed by where θ B is the rotation at the crack tip (root rotation), cf. Figure 3a. The area A DCBF depicted in Figure 4 is a function of the curvature at the crack tip, 1=R B . In the following, a normalized version of the curvature is used, where R p is the curvature at first yield. Thus, in Equation (4), A DCBF on the left-hand side is a function of k B and the right-hand side is a function of θ B . In Kinloch et al., [2] the deformation of the foil to the right of the crack tip ( Figure 3a) is estimated from a model consisting of a linear elastic beam on a linear elastic foundation, where the stiffness of the foundation is given by the transverse stiffness of (half of) the foil. This leads to a linear relation between θ B and k B . Attempts have been made to accomplish a similar solution for the present case, by considering the compliant substrate as an elastic foundation. However, due to substantial plastic strain in the foil and large root rotations, such a linear model has proven to be unsuccessful. Instead, the value of the root rotation, θ B , is obtained from experimental observations. Thus, k B is the only unknown and can be solved from Equation (4).
As mentioned previously, the present solution considers plastic deformation due to bending of the foil only. For small values of the peel angle, θ, plasticity due to stretching of the peel arm can be important. The smaller the value of θ, the smaller the values of θ B and k B . This means that the solution may be inappropriate for small values of k B . However, for the sake of completeness, the solution for all values of k B is presented here. As shown in Appendix C, the solution consists of three separate cases: Case 1: For 0 < k B < 1 peeling involves only elastic bending of the foil, i.e. no plastic dissipation is involved, G p ¼ 0. For this case, points B and A coincide and points C and D coincide with point O. Thus, the area A DCBF in Figure 4 is of triangular shape. As shown in Appendix C, Equation (4) takes the form, after a minor rearrangement. Case 2: For 1 < k B < 2 peeling involves elastic-plastic loading and elastic unloading, but no reverse plasticity. This means that, also for this case, A DCBF is triangular. As shown in Appendix C, Equation (4) takes the same form as for Case 1 above. That is, Equation (7) applies also for Case 2. The plastic energy dissipation rate 2 of Equation (3) takes the form, is used as a normalizing factor. Case 3: For k B > 2 reverse plastic deformation is involved. As shown in Appendix C, Equation (4) takes the form, and the plastic energy dissipation rate of Equation (3) takes the form, Summing up, the procedure to calculate the fracture energy is as follows: Equations (7) and (10) are separately solved numerically for k B . The two roots are denoted k B1 and k B3 , respectively. If 0 < k B1 < 1, Case 1 applies and G p ¼ 0. If 1 < k B1 < 2, Case 2 applies and G p is given by Equation (8) with The fracture energy is finally calculated using Equations (1 and 2).

Optimization of cohesive parameters
A FE-model of the peel test was used in order to obtain more detailed information about the cohesive parameters of the interface in the Al-LDPE laminate. The FE-simulations were performed on the θ ¼ 90 and θ ¼ 180 geometries used in the experimental setup, cf. Section 5. The input parameters of a trapezoidal cohesive law were optimized to the experimental results in terms of the steady-state peel force and root rotation. An ANN was trained based on simulation results and integrated inside a GA optimization code in search of an optimal cohesive parameter set (optimum G c , σ c0 and r). 2 Note that G p in Equation (8) differs slightly from Equation (11) in Kinloch et al. [2] Since Case 2 does not involve any reverse plasticity, they should be equal. It is noted that, in Equation (8) for all values of α, G p ! 0 when k B ! 1, as expected. This is not the case for Equation (11) in Kinloch et al. [2].

FE-modelling of peel
The peel setup was modelled using beam elements for the Al-foil and zero thickness cohesive elements for the interface between the Al-foil and the LDPE-film. Connector elements were used to connect the beam elements and the cohesive elements, cf. Figure 1. The results of this simplified model have been compared to a more detailed FE-model, where all layers present as described in Section 2 and modelled with 2D elements. The two models show comparable results, which motivates the use of the computationally less expensive beam model. For the beam theory solution described in Section 3, the maximum curvature and bending moment were assumed to be reached at the current position of the crack tip. For the FE-model, no assumption regarding the position of the maximum curvature is made. However, θ B denotes the rotation at the crack front also for the FE-model. The dimensions of the FE-model were chosen as small as possible, still capable of achieving steady-state conditions. The peel arm length in the initial configuration was 4 mm and the length of the portion of the Al-foil connected to ground through cohesive elements was 1 mm. The element length in these two portions of the model were 4 µm and 1 µm, respectively. The out-of-plane thickness was 15 mm. Beam elements of Timoshenko type was used (Abaqus element type: B21). An elastic-plastic constitutive model based on Hooke's law, von Mises linear kinematic hardening, was used to model the Al peel arm. The details of the model are described in Section 5.
The cohesive elements (Abaqus element type: COH2D4) of the interface were modelled with a trapezoidal cohesive law. A trapezoidal cohesive law was advantageous over a commonly used bi-linear law in this case as the size of the plateau ( Figure 5) can control the softening and hence root rotation. This was important in the optimization scheme as it could decouple peel force and root rotation to some extent. A short description of the cohesive law is given here for the case of pure normal separation, i.e. no shear separation. The normal cohesive separation is denoted δ and the normal cohesive stress σ c . The key cohesive input parameters in Abaqus [38] were the fracture energy (G c ), the interfacial strength (σ c0 ), the initial stiffness (c) and a ratio (r) governing the shape of traction-separation relation, cf. Figure 5. The ratio is defined as, For a given set of G c , σ c0 and c, the effect of the ratio r is as demonstrated in Figure 5.
The separations defining the cohesive law are expressed in the input parameters according to Equations (13-15) To respect the condition δ 1 < δ 2 < δ 3 , r should be within the following limit, In Abaqus, the cohesive law is defined through a damage variable, D. This means that the cohesive response is written as For the present trapezoidal cohesive law, the damage evolution, for a monotonically increasing δ, takes the form, In Abaqus, the input data are given as tabular values of D as a function of δ À δ 1 . Multiple points of values between the δ 1 to δ 2 and δ 2 to δ 3 ranges were necessary for a smooth simulation response. The cohesive separation in the peel test is predominantly normal separation, i.e. only small shear separation is encountered. In the cohesive model, it is assumed that the traction-separation relation for pure shear separation is the same as for pure normal separation. For mixed mode loading, a damage initiation criterion of the cohesive element was chosen as quadratic nominal stress (Abaqus type: QUADS). The damage evolution and hence the fracture energy is defined to be independent of the mode mix.
At the end of the simulation, the steady-state peel force response was recorded. Measurement of the root rotation involved an evaluation of the steady-state displacement field. First the damage variable, D (In Abaqus: SDEG), was requested along the line of cohesive elements, cf. the solid line in Figure 6. The cohesive zone is identified as the region where 0 < D < 1. The crack tip was identified as the point in the cohesive zone where D reaches unity, cf. the circle in Figure 6. The dashed line in Figure 6 shows the variation of the (normalized) rotation of the beam element nodes along the horizontal direction. The root rotation was obtained as the rotation of the beam element node connected to the cohesive element node located at the crack tip.
When using the FE-model to evaluate the experiments, the steady-state peel force, P, and the steady-state root rotation, θ B , are used as input data. The three parameters, G c , σ c0 and r, of the cohesive law are optimized to achieve the input values of P and θ B . The fact that P is known means that G e is known according to Equation (2). In order to determine the fracture energy, G c , the plastic energy dissipation rate, G p , must be determined, cf. Equation (1). In the evaluation of experiments (optimization), knowledge of the value of θ B , is used to determine the parameters of the cohesive laws and eventually G p . In order for this method to be successful, there must exist a correlation between θ B and G p . To this end, a numerical study was performed to illustrate the correlation. In that study, the cohesive strength, σ c0 , and the shape parameter, r, are varied, while the fracture energy, G c , is held constant. The steady-state peel force is recorded and G e is calculated according to Equation (2). Equation (1) is used to calculate G p ¼ G e À G c . The results are shown in Figure 7, where G p =G c is plotted as a function of σ c0 for two different values of r, cf. the solid and dashed curves. As seen, for a given value of r, G p varies considerably with σ c0 . For a given value of σ c0 , the value of r will also affect the value of G p , but to a lesser extent. This shows that G p is highly dependent on the shape of the cohesive law. It can also be seen from Figure 7 that θ B is highly dependent on σ c0 and r, cf. the dotted and dash-dotted curves. The general trend is that, the higher the values of σ c0 and r, the smaller the value of θ B and the larger the value of G p . Thus, G p and θ B are correlated, which motivates the use of θ B as a means to estimate G p .
The simulation/optimization process was set such that MATLAB could automatically take the cohesive parameters and make the calculations for input data, write the input file, read the necessary outputs, make the calculations for root rotation, take decision (optimize) based on simulation results and finally update the input file again until a desired accuracy was reached.

Optimization framework
The optimization framework of this study was divided into three sequential steps and is summarized in Figure 8. The method of using an ANN inside a GA was successfully adopted. A narrow range of the input parameters were   Table 3 in Section 6. These values were based on the fracture energy estimate obtained with the beam theory solution described in Section 3, followed by a coarse exploration by FE-simulation (step 1 in Figure 8).

Training the ANN
A well trained and validated ANN can replace the necessity of additional costly FE-simulations. The chosen ranges of cohesive parameter inputs (G c , σ c0 , r) represents the design space (DS) ( Table 3). FE-simulation outputs (peel force and root rotation) at points in the DS (that corresponds to sets of three input cohesive parameters) were generated. In order to select the distribution of simulation input variables in the DS, the Latin hypercube sampling method [39] was used. Such 30 different peel simulation input sets representing 30 points in the DS and corresponding outputs were found to be sufficient for initial training of the ANN. Latin hypercube sampling method is good for generating a small yet representative sample of cases. [40] The trained network was hoped to behave as a close enough function to output the steady-state peel force and root rotation response when any set of cohesive parameters is called by. A three-layered feed forward ANN with back propagation can be used to train any non-linear relationship with arbitrary accuracy. [41] Hence, an input layer with three neurons, one hidden layer with six neurons and an output layer with two neurons were adopted. For 30 training simulations with chosen inputs and output numbers, six hidden layer neurons were checked to be ideal in the present case. This is more than the maximum recommended hidden layer neuron number, according to Hagan. [42] Here, N i = number of input neurons, N o = number of output neurons, N s = number of samples in training data set and β = an arbitrary scaling factor usually between 2 and 10. Training, validation and testing ratio of 85/100, 10/100 and 5/100 were used respectively. Levenberg-Marquardt backpropagation algorithm in MATLAB optimization toolbox was found to increase the accuracy of ANN predictions. This ANN training procedure is depicted in step 2 of Figure 8. The simulation cases were run in Abaqus 6.14 with an Intel Xeon processor workstation (6 cores, 2.50 GHz). For beam element model, 30 peel simulations (90 or 180 ) took about 10 hr to complete. When simulation results were checked with a more detailed 2D FE-model, several times higher simulation time was required.

Optimization by the genetic algorithm (GA)
GA is a gradient-free optimization method for global optimization which is based on stochastic approaches. This is an efficient method for non-linear, non-differentiable objective functions. [32,43] GA requires a large number of function evaluations which is a major limitation if the function evaluation involves time costly simulations. [32] Like in this study, performing that much simulation can be an obstacle. One solution to this problem is to replace the simulation by a well trained ANN. The optimization objective was set to minimize the sum of the mean square difference between ANN output values and experimental values (peel force and root rotation), where, obj is the objective values for peel force (obj(1)) and root rotation (obj(2)) ( Table 1). I j is a set of cohesive input parameters where I 1 is fracture energy (G c ), I 2 is interfacial strength (σ c0 ) and I 3 the ratio (r). The function net I j À Á Â Ã is expected to produce simulation alike peel force and root rotation for given input set I j from the trained ANN, net. Finally, the squared sum of the differences between net responses and objectives are normalized by the sum of the square of both objectives to obtain a relative error. This classical aggregative method, which combines all objectives into a weighted sum was used to get a singular optimum solution. GA code from MATLAB Optimization Toolbox [44] was implemented. Scattered type crossover and Gaussian mutation function was adopted. The termination criterion was set to an average change in fitness value less than 1e À7 . The ranges of variables used for optimization were as presented in Table 3.
The new set of optimized parameters after every GA evaluation were used for additional FE-simulation to test the prediction of the current ANN. If this new simulation results were found to be acceptable the optimization procedure was stopped and if not, the new simulation result was added to the ANN training data set and the network was re-trained. Always the default MATLAB weights and biases were used for consistent training of the ANN with additional simulations in this step. New prediction of optimized cohesive input set from the ANN was again tested by FE-simulation and the procedure was repeated until a desired accuracy was obtained, cf. step 3 of Figure 8. After 10 iterations, optimum cohesive parameters for 90 peel were found and results are presented in Section 6. Similarly, optimum cohesive parameters for 180 peel simulation are obtained. Optimized cohesive parameter sets for 90 and 180 peel were different by some margin. Finally, a common cohesive parameter set that satisfies 90 and 180 peel objectives was sought. It took minimization of the summation of the error functions for 90 and 180 peel.

Test results
Peeling of Al-foil away of LDPE-film was performed at the laboratory of Division of Structural Mechanics, Lund University, Sweden, and were reported in Bruce and Holmqvist. [8] An Instron tensile testing machine was used to perform the peel tests. The laboratory environment was controlled at 23°C and 50% humidity for 24 hr before specimen preparation and peel testing. Peel specimens were cut into 15 mm wide strips. For convenience of mounting the specimens, the paperboard that is attached to the laminate during actual package production is kept intact. [8] Three peel angles were studied; 0 , 90 and 180 . The test setup for peeling of 180 and 90 tests are illustrated in Figure 9. The 0 peel test was reportedly not possible to complete because the Al-foil failed before the adhesion could be separated due to the foil stresses becoming too large. Results on peel force response and root rotation for 90 and 180 peel were reported in Bruce and Holmqvist. [8] The results are also given here. Figure 10 shows the force-displacement response. A summary of the steady-state peel force and root rotation is given in Table 1.
The measurement of the rotation at the crack front, θ B , is indicated in Figure 11. The root is defined as the point where the foil separates from the substrate. Note that the substrate exhibits a residual deformation for the points where peeling has taken place. An on-screen protractor was used to measure the root rotation angles from the digital photographs.
It should be noted that the peel arm visible in Figure 11 consists of the aluminium layer, the inside LDPE layer and the inside adhesion layer, cf.   due to deformation of the inside LDPE layer. The inside LDPE layer is located on the compressed side of the peel arm, which leads to an increase of the thickness due to the Poisson's effect. This increase in thickness is larger near the root where the curvature is large.
In order to evaluate the peel experiments, the stress-strain relation for the Alfoil is needed, cf. Sections 3 and 4. Tensile testing of the 6.3 µm thick Al-foil was performed in the work presented in Sharif and Majeed. [45] Due to anisotropy inherited during the production process of the substrates and lamination, the final laminate is anisotropic. The tensile test data were obtained from the machine direction (rolling direction during production). In subsequent analysis of peel test data, all the materials were treated as homogeneous and isotropic with the properties obtained for the machine direction. Figure 12 shows the stress-strain curve. Analysis of the initial linear elastic part gives an elastic modulus 3 of E ¼ 55 GPa. The Poisson's ratio of the specific Figure 12. Stress-strain relation from Al-foil tensile test [45]. Figure 11. Root rotations from experiments. [8] (a) θ ¼ 90 and (b) θ ¼ 180 . 3 The value of E differs considerably from bulk values. The reason for the discrepancy is that the thickness of the foil is less than the normal size of aluminium crystals, cf. Andreasson et al. [3].
Al-foil has been reported to be 0.3. [45] The yield stress of Al-foil was σ y ¼ 63 MPa. For strains higher than 0.008, a behaviour close to linear hardening is displayed, cf. Figure 12. The linear hardening parameter is estimated to α ¼ 0:012, cf. Figure 2. Unfortunately, the Al-foil fails at a relatively small value of the strain, i.e. about 0.014. This is due to local failure emanating from small cracks at the edges of the foil. The strain due to bending of the foil during the peel tests, is higher than 0.014. In lack of information about the stress-strain relation at higher values of the strain, linear hardening is assumed also for strains higher than 0.014.

Results and discussion
The methods described in Sections 3 and 4 are here used to evaluate the peel experiments reported in Bruce and Holmqvist [8] and in Section 5.
In the beam theory solution of Section 3, the experimental values of the steady-state peel force, P, the steady-state root rotation, θ B , and the peel angle, θ, are used as the input data. The results for the two different peel angles are given in Table 2. For both peel angles, the normalized root curvature k B > 2, i.e. reverse plastic deformation is encountered. This means that Case 3 in Section 3 is used to calculate k B and the plastic dissipation rate, G p . The method yields similar values of fracture energy, G c , for the two peel angles, with the one for θ ¼ 180 being slightly higher than for θ ¼ 90 .
Based on the estimated fracture energy in Table 2 and a coarse exploration by FE-simulation (step 1 in Figure 8), a narrow range of the input parameters were selected for further exploration by ANN and GA, cf. Table 3.
Through the FE-optimization scheme presented in Section 4, the fracture energy (G c ), interfacial strength (σ c0 ) and the cohesive displacement ratio (r) of a trapezoidal cohesive law were estimated. The peel root rotation in addition to the steady-state peel force was considered as optimization objective. The use of two objectives, increased the chance of the optimum cohesive parameters to be more unique. A common set of values for (G c ), (σ c0 ) and (r) could not be achieved that satisfied both the 90 and the 180 test results very Table 2. Evaluation results for beam model.  closely, cf. Table 4 and Figure 13. However, when the cohesive parameters were optimized for a single peel angle individually, a close agreement of peel steady-state force and root rotation was achieved between FE-simulation and experimental results cf. Table 4 and Figure 13. The optimized G c for 180 peel in this case was higher than that for 90 peel; similar to the analytical estimation, cf. Table 2.
The cohesive model used here is defined to have a mode-independent fracture energy. The FE-simulations show that the case of 180 is associated with a larger part of the fracture energy being related to shear separation than for the case of 90 . Several earlier studies in different laminates reported shear fracture energy to be more than that of for normal fracture. [46][47][48] Thus, the difference in G c for different peel angles, suggests that the fracture energy can be different in normal and shear mode of fracture.
It is noted from Tables 2 and 4 that higher fracture energies are predicted with the beam theory solution than for the FE-solution. This discrepancy is due to differences in the estimated plastic energy dissipation rate, G p . The source of this discrepancy is the assumption of the beam theory model, that the curvature and bending moment are maximum at the crack front. This assumption is justified for a rigid substrate. However, the FE-solution shows that, for the present compliant substrate, the point of maximum curvature is located well inside the  cohesive zone, cf. Figure 14. Due to the direction of the peel force, this behaviour is more pronounced for 90 peel than for 180 peel. The effect of the fact that the maximum curvature is located ahead of the crack tip is depicted in Figure 15, where the moment-curvature history is shown for the beam theory solution (solid curve) and the FE-solution (dashed curve), respectively. The points corresponding to the crack tip are denoted B and B 0 , respectively. The relation between the two curves can be explained by considering Equation (4). The root rotations, θ B and θ B 0 , are equal for the two solutions. Also, the peel force, P, is the same for both solutions. Hence, the right-hand side of Equation (4) is the same for the two solutions. According to Equation (4), this means that the two curves in Figure 15 are related such that the areas A DBF and A DB 0 F 0 are equal. As mentioned above and seen in Figure 15, the fact that the maximum curvature is located inside the cohesive zone, means that the plastic energy dissipation rate, G p , is underestimated by the beam theory solution, cf. Figure 3. This issue has previously been addressed in Kawashita et al. [49] where the curvature is measured using an experimental method based on high resolution digital photography. The point of maximum curvature was found to be located inside the cohesive zone and that point was defined as an effective crack tip. The measured maximum curvature was used to find the plastic energy dissipation rate and ultimately the fracture energy.
In the present study, a different path is taken in that the root rotation is measured and defined as the rotation of the actual crack tip. This method does not require high resolution digital photography. From optimization of the FE-model, the cohesive parameters are determined using the steady-state peel force and root rotation. Based on the above, we may also conclude that the beam theory solution may well be used to obtain an estimate of the fracture energy, but that the more elaborate FE-solution is needed to obtain more reliable values of the fracture energy. Moreover, the FE-solution provides additional information about the cohesive properties.

Conclusions
In this work, peeling of a metal foil from a compliant substrate is studied. To this end, two supplementing models, analysing test results from steady-state peel, are introduced. The input data from steady-state peel test results are the peel force and the root rotation. The first model is based on large displacement beam theory and provides analytical expressions for the fracture energy associated with peel. It is a modification of the method developed in Kinloch et al. [2] The method presented here is based on a linear kinematic hardening model of the peel arm, unlike the linear isotropic hardening model used in Kinloch et al. [2] As in Kinloch et al. [2] it is assumed that the curvature and bending moment in the peel arm are maximum at the crack tip. The second model is a finite element (FE) model, where beam elements are used to model the peel arm and the substrate is modelled using cohesive elements. With this model, the properties of a trapezoidal cohesive law are extracted from test data. Besides the fracture energy, the properties consist of the cohesive strength and a shape parameter. A method to determine the cohesive properties is developed. In the method, an ANN is trained based on simulation results and integrated inside a GA optimization code.
The two methods are used to analyse data from tests where an aluminium (Al) foil is peeled from a LDPE film. The tests are conducted using the peel angles 90 and 180 and have previously been reported in Bruce and Holmqvist. [8] The more detailed FE-model revealed that the maximum curvature and bending moment is located ahead of the crack tip, inside the cohesive zone. Thus, the assumption of the analytical model introduces errors, which leads to an underestimation of the plastic energy dissipation rate. Thus, the fracture energies determined with the analytical model are slightly higher than for the FEmodel. Both methods indicate that the fracture energy is slightly higher for the 180 peel angle than for 90 . This indicates that the fracture energy is modedependent. Reported orders of magnitude of the fracture energy and cohesive strength are 60-70 J=m 2 and 5 MPa, respectively.
Thus, only the parts of C having either N 1 Þ0 or a non-zero T i contributes, for this case only the loaded end of the foil where The deformation gradient is formed using polar decomposition, where R ik is the rotation tensor and U kj is the right stretch tensor. For the present case, with a coordinate system aligned with the peeling arm, cf. Figure 16a, where λ i are stretches; λ 1 in the direction of the peeling arm, λ 2 perpendicular to the peeling arm and λ 3 out of the plane. Inserting Equations (A9) and (A10) into Equation (A8) and comparing the result with Equations (A2) and (A3) yields Now, Equations (A5-A7), (A11) and (A12) are inserted into Equation (A4), resulting in, where λ 1 and W are evaluated at the loading point. The stretch is defined as where a is the engineering strain. The present study is confined to cases where no plastic deformation takes place at the uniaxially loaded end of the peel arm. That is, a < y and hW ¼ P a =2b. Equation (A13) then takes the form, For small peel angles, the term including a is important. Here, the peel angle, θ, is assumed to be large enough to ensure that 1 À cos θ > > a =2. Under these conditions the energy release rate takes the simple form, The assumptions made regarding the size of a is also motivated by the fact that the method to determine the plastic energy dissipation rate, G p , is based on the neglection of plastic strain due to stretching of the peeling arm.
and the normalized stress, means that the normalized bending moment in Equation (B2) takes the form, For a foil subjected to pure bending deformation, the strain distribution is linear through the thickness of the foil, where R is the radius of curvature due to bending. The curvature of the bended foil is defined as 1=R. A normalized curvature, k, is defined in Equation (5). From Equations (B7) and (5), it follows that For monotonic loading (increasing ), the stress-strain relation can be written as Here, a bilinear stress-strain relation is used, Thus, for monotonically increasing k, the normalized stress is given by Thus, the stress-strain relation σðÞ shown in Figure 2 is translated to σðkζÞ The loading-unloading sequence is divided into four stages related to Figure  (ii) Elastic-plastic loading (AB): When k is monotonically increased beyond k ¼ 1, a plastic zone is formed starting at the outer fibres of the foil. For a given value of k, the current interface between the elastic and plastic zones is situated at ζ ¼ 1=k, cf. This solution is valid until reverse plastic deformation is initiated at the outer surface of the foil. For the kinematic hardening model employed here, the condition for this is σðζ ¼ 1Þ ¼ f ðk B Þ À 2. With Equation (B17) the value of k at the onset of reverse plasticity is (iv) Reverse elastic-plastic unloading (CD): As the curvature is decreased below k C , there are three different zones through the thickness of the foil: (a) the elastic core, 0 < ζ < 1=k B , (b) the zone of elastic-plastic unloading, 1=k B < ζ < ζ, and (c) the zone of reverse plasticity, ζ < ζ < 1. For a specified value of k, the position of the interface between zones (b) and (c) is given by For a given position, ζ, reversed plasticity starts at the following value of the curvature, With these new entities, the stress distribution is written as σðζÞ ¼ kζ for 0 < ζ < 1=k B f ðk B ζÞ þ ðk À k B Þζ for 1=k B < ζ < ζ f ðk B ζÞ À 2 þ ðk À kÞαζ for ζ < ζ < 1  With m OA , m AB , m BC and m CD from Appendix B, Equation (11) follows. The area A DCBF in Figure 4 is mathematically expressed as It may be expressed in terms of the normalized bending moment, m, defined in Equation (B2), where G n is defined in Equation (9). For Case 1 and Case 2 (see Section 3), the area is triangular and the slope of DCB is m 0 ðkÞ ¼ 2=3, Thus,