Thermal Delamination Modelling and Evaluation of Aluminium–Glass Fibre-Reinforced Polymer Hybrid

This paper aims to propose a temperature-dependent cohesive model to predict the delamination of dissimilar metal–composite material hybrid under Mode-I and Mode-II delamination. Commercial nonlinear finite element (FE) code LS-DYNA was used to simulate the material and cohesive model of hybrid aluminium–glass fibre-reinforced polymer (GFRP) laminate. For an accurate representation of the Mode-I and Mode-II delamination between aluminium and GFRP laminates, cohesive zone modelling with bilinear traction separation law was implemented. Cohesive zone properties at different temperatures were obtained by applying trends of experimental results from double cantilever beam and end notched flexural tests. Results from experimental tests were compared with simulation results at 30, 70 and 110 °C to verify the validity of the model. Mode-I and Mode-II FE models compared to experimental tests show a good correlation of 5.73% and 7.26% discrepancy, respectively. Crack front stress distribution at 30 °C is characterised by a smooth gradual decrease in Mode-I stress from the centre to the edge of the specimen. At 70 °C, the entire crack front reaches the maximum Mode-I stress with the exception of much lower stress build-up at the specimen’s edge. On the other hand, the Mode-II stress increases progressively from the centre to the edge at 30 °C. At 70 °C, uniform low stress is built up along the crack front with the exception of significantly higher stress concentrated only at the free edge. At 110 °C, the stress distribution for both modes transforms back to the similar profile, as observed in the 30 °C case.


Introduction
Fibre metal laminates (FMLs) are increasingly used in aeronautical and marine applications for their combination of advantages from two different materials. FMLs involve adhesively bonding alternating metal and composite layers to form a laminated sandwich structure. While pure composite materials are prone to brittleness, temperature effects and barely visible internal damage, pure metals are susceptible to fatigue loading and lower yield strength. However, the hybrid nature of the FML allows beneficial properties from both metal and composite counterparts, including lighter weight, ductility, fatigue resistance and damage tolerance [1][2][3]. The multiple layers of the FML structure prevent crack propagation from damaged layers to subsequent layers via crack bridging and prevent catastrophic failure that would otherwise occur unimpeded in monolithic materials. Similarly, multiple alternating laminates are capable of suppressing the spread of thermal and moisture absorption across the layers that allow potential applications of FMLs in hand, Gunther et al. [22] developed a temperature-dependent cohesive zone model for bi-metal laminates from 23 to 150 • C by performing tubular butt joint tests both in tension and torsion setups, a tapered double cantilever beam and tapered end notched flexure tests. Tubular butt joint tests were performed to extract adhesive peak loads with the means of calculating cohesive peak traction stresses, while concurrently obtaining penalty stiffness from 10% of the total strain.
There are only a few known studies on delamination with simultaneous temperature effects on dissimilar laminates. Qin et al. [23] experimentally and numerically researched the effects of temperatures at −40, 23 and 80 • C on a variety of adhesive joints between CFRP and aluminium. They applied fitting quadratic stress criteria to obtain cohesive failure stresses, estimated the fracture toughness values from the literature and assumed that the ratio of Young's and shear moduli to be the same across different temperatures. Fibre metal laminates will continue to be prevalent in broad ranges of vehicles, including exposure to high temperatures. A methodology is required to procedurally obtain and effectively analyse the Mode-I and Mode-II delamination properties between dissimilar materials of FMLs at higher temperatures. The findings aid in ensuring proper suitability of materials and adhesives based on different applications.
In this study, the temperature effects on Mode-I and Mode-II delamination of hybrid glass fibre-reinforced polymer laminated with aluminium alloy are numerically investigated. Explicit nonlinear code LS-DYNA is employed to construct a finite element model of the experimental double cantilever beam (DCB) and end notched flexural (ENF) tests. The adhesive interface is modelled based on cohesive zone modelling defined with the bilinear traction separation law. The main cohesive parameter, fracture energy, is characterised by data reduction of experimental DCB and ENF load-displacement results by using the Irwin-Kies equation. The finite element models are then compared and validated with the experimental tests. Temperature effects on the stress distribution of the delamination crack front are investigated through validated FE models.

Material Modelling Constitutive Equations
Materials of the FML hybrid comprise metal and a composite part which, in the current study, are aluminium 2024-T3 and glass fibre-reinforced polymer (GFRP), respectively. Dissimilar mechanical properties require separate material modelling to accurately capture their behaviour and changes across different temperatures. A cohesive zone model is used to represent delamination between aluminium and GFRP.

Aluminium Material Modelling
The aluminium part is modelled with a Johnson-Cook material model that can sufficiently capture the effects of temperature. Flow stress, σ y , is expressed as [24]: where A, B and n are the plastic strain material constants, ε p is the effective plastic strain, C is the strain rate effect constant, . ε is the effective total strain rate, . ε 0 is the quasi-static threshold rate, m is the temperature effect constant, T is the effective room temperature, T room is the room temperature and T melt is the melting point of the aluminium. The fracture strain with the effect of pressure, strain rate and temperature is: where d 1 , d 2 and d 3 are the pressure effect constants, η is the ratio between pressure and effective stress, d 4 is the strain rate effect constant and d 5 is the temperature effect constant.

GFRP Material Modelling
On the other hand, GFRP is modelled using Chang-Chang failure criteria. Each failure mode is initiated accordingly depending on whether each of its corresponding equation shown below has been met [25]. However, since the Chang-Chang model does not include temperature effects, GFRP material properties at each temperature are included accordingly based on the respective temperature.
For the tensile fibre mode: For the compressive fibre mode: For the tensile matrix mode: For the compressive matrix mode: where σ 1 and σ 2 represent the effective stress tensor components in respective directions, τ 12 represents the effective shear tensor, X T and X C are the longitudinal tensile and compressive strengths, Y T and Y C are the transverse tensile and compressive strengths, S C is the shear strength, E A and E B are the longitudinal and transverse Young's modulus, G AB is the shear modulus, υ BA and υ AB are the Poisson's ratio in respective directions, and β is the weighting factor for shear term.

Cohesive Zone Modelling
The interface is defined as 8-point cohesive elements with the same thickness of 13 µm as the Teflon material in the laminate pre-crack. By implementing the bilinear traction separation law in the cohesive elements, the process of delamination can be captured correctly. As shown in Figure 1, delamination behaviour, including increasing traction forces, followed by reaching the peak of interface strength, then the degradation of the interface and finally the full separation of the laminate can be captured by the traction separation law. Similar to Chang-Chang failure, the cohesive model implemented does not capture temperature effects; hence, the inputs of each temperature are defined separately. To ensure proper stress transfer between cohesive and neighbouring material, the cohesive elements are aligned to their adjacent material elements with the same mesh size to share the same node positions. Subsequently, the coinciding nodes between the cohesive and adjacent material nodes are merged into a single node. The Mode-I energy release rate is defined as [26]: where and δN are the normal direction peak traction and ultimate displacement. EN is defined as the Mode-I penalty stiffness. The Mode-II energy release rate is defined as [26]: where τ and δT are the tangential direction peak traction and ultimate displacement. ET is defined as the Mode-II penalty stiffness.

Finite Element Modelling
Commercial LS-DYNA/Explicit software is used in this research to develop finite element models of Mode-I and Mode-II delamination with a combination of material models for aluminium and GFRP and a cohesive zone model. Experimental tests were conducted beforehand, where they are reported thoroughly in previous literature [27,28], with dimensions, materials and boundary conditions matching the finite element model. The results are then compared with numerical predictions to validate the models.
Constant stress solid elements are used to define the aluminium and GFRP parts. Aluminium 2024-T3 mechanical properties based on the Johnson-Cook material model are shown in Table 1, while the material properties of GFRP previously extracted from experimental tests at each temperature are in Table 2. Base material properties of aluminium, including density, Young's modulus, shear modulus and Poisson's ratio values, are acquired from Wen et al. [29]. The Johnson-Cook yield surface parameters and failure parameters are adopted from Lesuer [30] and Li et al. [31]. Material properties of GFRP at each temperature are adopted from previous work [28]. The Mode-I energy release rate is defined as [26]: where σ and δ N are the normal direction peak traction and ultimate displacement. E N is defined as the Mode-I penalty stiffness. The Mode-II energy release rate is defined as [26]: where τ and δ T are the tangential direction peak traction and ultimate displacement. E T is defined as the Mode-II penalty stiffness.

Finite Element Modelling
Commercial LS-DYNA/Explicit software is used in this research to develop finite element models of Mode-I and Mode-II delamination with a combination of material models for aluminium and GFRP and a cohesive zone model. Experimental tests were conducted beforehand, where they are reported thoroughly in previous literature [27,28], with dimensions, materials and boundary conditions matching the finite element model. The results are then compared with numerical predictions to validate the models.
Constant stress solid elements are used to define the aluminium and GFRP parts. Aluminium 2024-T3 mechanical properties based on the Johnson-Cook material model are shown in Table 1, while the material properties of GFRP previously extracted from experimental tests at each temperature are in Table 2. Base material properties of aluminium, including density, Young's modulus, shear modulus and Poisson's ratio values, are acquired from Wen et al. [29]. The Johnson-Cook yield surface parameters and failure parameters are adopted from Lesuer [30] and Li et al. [31]. Material properties of GFRP at each temperature are adopted from previous work [28].
According to Figure 2, Mode-I delamination is modelled based on DCB while Mode-II is based on ENF experimental setups, as described in depth in previous publications [27,28]. Based on DCB and ENF experiments, the thickness of GFRP is modelled as 2.2 mm while aluminium is 2 mm. The primary fibre direction of GFRP is along the specimen length. The thickness of the main cohesive elements is modelled based on the thickness of nonadhesive Teflon in the experimental test, 13 µm. Hence, the total thickness of the hybrid model, 2 h,  Table 3, where the slope, k, and peak loads, F P , are used in data reduction to obtain fracture toughness.     For the DCB model, shell elements are used to model the plates of hinges. Based on Figure 2a, horizontal plates are fixed to the surface of the laminate while being hinged with vertical plates. Therefore, the edges of two shell planes perpendicular to each other are joined and constrained, permitting only revolutions along the joined tip axis to simulate the rotating hinge. The bottom vertical plate is fixed to disallow any movement or rotation, while the top vertical plate is only allowed to move in the axial direction where it pulls the joined horizontal plate via the hinge. On the other hand, the 3-point bending loading is modelled using 3 rollers, defined as solid elements, as shown in Figure 2b. Two bottom rollers are fixed while the top roller can only move downwards along the vertical direction, as shown.
The normal penalty stiffness, E N is calculated from [10]: where E is Young's modulus of the adhesive with 6.03 GPa and t is the adhesive thickness at 0.013 mm. The in-plane penalty stiffness, E T , is calculated from [10]: where G is the shear modulus of the adhesive at 0.90 GPa. Based on the Irwin-Kies data reduction equations, slope, k, and peak load, F P , from Table 3 are converted into the respective fracture toughness required for cohesive parameters at each temperature. Peak traction stresses, σ and τ, are estimated empirically from compilations of previous literature on Mode-I and Mode-II interface strength done by Zhao et al. [14], where the values are similar to a literature review of aluminium and GFRP laminates. To acquire the cohesive properties at higher temperatures, a ratio approach between experimental data and cohesive properties was implemented in this study.
The values of penalty stiffness, E N and E T , at higher temperatures are ratioed to the ratio of the slope, k, from a load-displacement curve based on Equations (11) and (12) [32].
where k I is the slope for Mode-I and k II is the slope for Mode-II; x is the respective temperature based on the required parameter.
On the other hand, the value of peak traction stresses is ratioed concerning the fracture toughness ratio between different temperatures based on Equations (13) and (14), where a similar method was used by Qin et al. [23] to acquire low and elevated temperature properties.
The resulting cohesive parameters at each temperature are shown in Table 4.
To simulate the quasi-static rate of loading with minimal kinetic effect and vibrations, the rate of loading must be set as low as possible. The top plate of the DCB simulation is set to move upwards at an increasing rate up to a maximum of 1 m/s. On the other hand, the ENF top roller moves downward at an increasing rate up to a maximum of 0.5 m/s. The loading rate at the time of crack initiation below 0.5 m/s ensures minimal effects on the cohesive properties [19]. To minimise computational time, the elements are reduced by using half models since both cases are symmetrical across the length. Boundary constraints are utilised with elements along the cut section of the half models. To control the hourglass effects of elements, Flanagan-Belytschko stiffness for hourglass control is defined for the aluminium. A mesh convergence study is conducted by comparing the load-displacement curve of the FE model to experimental results at room temperature. It was found that an element size of 0.5 mm was able to produce good results within acceptable computational time. Mass scaling is applied for the models within an acceptable range with negligible dynamic effects to increase the simulation time step and reduce the total simulation run time.

Numerical Results and Analysis
The following section discusses the validation and analysis of Mode-I and Mode-II delamination models at temperatures of 30, 70 and 110 • C. A mesh convergence study is highlighted in Section 4.1 to demonstrate the optimum mesh size that is capable of accurately capturing the delamination mechanism. In Section 4.2, validation is carried out using load-displacement curves obtained from the results. Slopes and peak load between experimental and simulation results are compared to demonstrate the feasibility of models by percentage agreements. Section 4.3 elaborates on the localised stress distribution along the crack tip and how temperature affects this stress. The stress contour shows how stress builds up along the crack front and location of propagation while a detailed plot of mode stress versus distance from the centre of the laminate is used to describe it in detail.

Mesh Convergence
The fundamental process of the mesh convergence study is carried out during FE modelling to ensure the model can produce accurate results. The element size within a model affects the model accuracy where a large mesh size captures exceedingly large stresses. By refining into smaller mesh sizes, the effects of discretisation can be reduced, thus producing a solution closer to real-life specimen counterparts. However, a large mesh number increases the FE simulation time. With a sufficiently fine mesh, the FE model will be able to produce accurate results where further refinements have negligible improvements, therefore preventing the need for finer mesh than optimum where simulation time becomes unnecessarily long.
Mesh sizes of increasing refinement of 10, 5, 2, 1, 0.5 and 0.25 mm were modelled for Mode-I. A similar material model, element type and boundary conditions were defined for all of the mesh sizes. The load-displacement curve of each mesh size is then plotted and compared with the experimental results. Figure 3 shows the load-displacement results of Mode-I models with different mesh sizes ranging from 10 mm to 0.25 mm compared to the experimental results (dotted marker). A noticeable convergence of slope from the coarser size of 10 mm to the finer mesh size of 0.5 mm is visually seen. Mesh refinement from coarse to fine mesh increases the slope to match the experimental gradient and, at the same time, the peak load decreases from exceedingly high towards the average experimental peak load. Such a trend converging towards the experimental curves indicates the feasibility of the current mesh convergence Polymers 2021, 13, 492 9 of 18 method, by comparing load-displacement results. The curves illustrate that mesh sizes above 2 mm are way off bound in both exceedingly low slope and high peak loads. While mesh sizes of 0.7 and 1 mm have peak loads that are comparable with the experimental data, both have slopes which are still comparably low. At the mesh size of 0.5 mm and smaller, the slope coincides with that of the experiment, where 0.25 mm does not exhibit much improvement on 0.5 mm. The convergence from coarser to a finer mesh by comparing the slope value is distinctly illustrated in Figure 4, where 0.5 mm matches the required experimental slope value while 0.25 mm overestimates it. Therefore, a mesh size of 0.5 mm is chosen for both the Mode-I and Mode-II FE models. The mesh convergence at mesh size less than 0.5 mm agrees with the findings by Turon et al. [7]. and compared with the experimental results. Figure 3 shows the load-displacement results of Mode-I models with different mesh sizes ranging from 10 mm to 0.25 mm compared to the experimental results (dotted marker). A noticeable convergence of slope from the coarser size of 10 mm to the finer mesh size of 0.5 mm is visually seen. Mesh refinement from coarse to fine mesh increases the slope to match the experimental gradient and, at the same time, the peak load decreases from exceedingly high towards the average experimental peak load. Such a trend converging towards the experimental curves indicates the feasibility of the current mesh convergence method, by comparing load-displacement results. The curves illustrate that mesh sizes above 2 mm are way off bound in both exceedingly low slope and high peak loads. While mesh sizes of 0.7 and 1 mm have peak loads that are comparable with the experimental data, both have slopes which are still comparably low. At the mesh size of 0.5 mm and smaller, the slope coincides with that of the experiment, where 0.25 mm does not exhibit much improvement on 0.5 mm. The convergence from coarser to a finer mesh by comparing the slope value is distinctly illustrated in Figure 4, where 0.5 mm matches the required experimental slope value while 0.25 mm overestimates it. Therefore, a mesh size of 0.5 mm is chosen for both the Mode-I and Mode-II FE models. The mesh convergence at mesh size less than 0.5 mm agrees with the findings by Turon et al. [7].      curve plot can be divided into three main stages. The first stage is the initial increase in load, indicating a build-up of stress across the specimen. As the cohesive element at the crack tip reaches maximum stress, the peak load is reached, where the crack initiates with the first load drop. After an initial load decrease, the occurrence of crack propagation is associated with a general trend of decreasing load. The fluctuation of load after an initial drop is associated with uneven crack propagation, caused by stress concentration build-up and release at the crack tip [33]. Based on experimental results, the increasing temperature has minimal effects on the slope from 30 to 70 • C, similar to findings by Davidson et al. [20]. However, the similarity ends at the highest temperature of 110 • C, with a large drop of the slope in the current study, which is postulated to be caused by different materials. Differences between Mode-I slopes of each temperature are tabulated in Table 5 with a maximum % difference of not more than 6%.    There is a larger discrepancy between peak loads but, in general, the FE-predicted load curve lies between the fluctuations of the experimental data. For 30 and 70 • C, the experimental results showed sudden large drops of load after the peak, while the drop is more gradual for the finite element models. The abrupt drop during experiments indicates that more energy is accumulated before a sudden release, while, in simulations, the crack grows progressively from element to element. Hence, lower peak loads are predicted for FE. The peak load difference is indicated in Table 5, with 7.41 and 14.76% for 30 and 70 • C, respectively. On the other hand, the 110 • C experimental peak load drops are much smoother compared to the FE counterpart. Therefore, the crack growth region is much better simulated compared to 30 and 70 • C.

Load-Displacement Curves and Model Validation
The load-displacement curves of Mode-II delamination at each temperature are illustrated in Figure 6a-c. Unlike Mode-I, where the load decreases after crack propagation occurs, Mode-II is characterised by load increase after an initial sudden load drop as the crack propagates to the side rollers, preventing any further cracks. Therefore, compared to Mode-I, Mode-II crack initiation does not possess fluctuating loads after a peak load drop. The simulated results appear to agree acceptably with their experimental counterparts, wherein the slope matches with very minor discrepancies. Table 6 shows the corresponding % difference for 30, 70 and 110 • C FE load-displacement curve slopes at 4.12, 7.26 and 3.99%, comparable to the reports by Qin et al. [23] and Gunther et al. [22], despite distinct methods for acquiring penalty stiffness. Thus, this further reaffirms that the penalty stiffness has minor effects on the overall resulting slope validation. For 70 • C, the trend only follows up until the crack initiation, as the load drop is much apparent in the simulation model, whereas the experimental load drop is only up to about 10% before the load increased again.
In similar fashion to Mode-I, the instance of crack initiation is much more difficult to predict, therefore creating a larger divergence between peak loads. The maximum difference in prediction of peak load is for 30 • C, with 15.12%. The larger scatter of peak load is even apparent between the three experimental datasets, only becoming less so at higher temperatures. It should be noted that the FE model slightly overestimates the peak load of 110 • C compared to experimental results for both Mode-I and Mode-II due to a much sharper load drop. The large instantaneous drop of load at lower temperatures could be attributed to brittleness in the composite that leads to a sudden release of the energy during crack propagation. However, for 110 • C, the higher temperatures increased the overall ductility of the adhesive, making crack propagation much smoother [19]. The effects of temperature on the increased ductility and gradual stress release during crack initiation are less well captured in cohesive models. Polymers 2021, 13, x FOR PEER REVIEW 13 of 19      Figure 7 shows the stress contour cut out of cohesive elements near the crack tip, while Figure 8 plots the stress distribution across the crack tip normalised with maximum stress versus the location from the cohesive centre. Due to the FE half model implemented, the stress values are mirrored across the centre. While the stress contours all show the same characteristic of crack initiation from the centre location towards the sides, a substantial change can be noticed by the way stress distribution evolves towards higher temperatures. At 30 • C, elements concentrated at the centre possess much higher stress whilst elements 4 mm and further have much lower stress, about 25% lower than the centre. This stress localisation drastically changes during a temperature of 70 • C where all elements except the outermost (9 mm) from the centre reach maximum Mode-I stress. With the highest temperature of 110 • C, the stress distribution is still akin to that of 70 • C but the stress drops slightly at a distance of 4 to 8 mm (less than 10%). At higher temperatures, stress progression is more rapidly spread from the centre to the outer elements.

Crack Initiation and Stress Distribution
Crack front stress distributions of Mode-II FE models at temperatures of 30, 70 and 110 • C are exemplified in Figures 9 and 10. Stress contours from Figure 9 show Mode-II typical stress and a crack initiating from the sides of the ENF specimen comparable to analyses by Koloor and Tamin [11]. The corresponding normalised stress values across the crack-initiating front are plotted in Figure 10. At 30 • C, localised normalised stresses are fluctuating more across the crack tip with values of 0.85 to 0.93, while the edge element has significantly higher shear stress. As is with Mode-I, 70 • C Mode-II stress features the trend where all elements from the centremost to an 8 mm distance have the same stress, only contrasted with all shear stresses exhibiting the lowest values, while the elements at 9 mm have substantially higher shear stress by 25%. Such a phenomenon is intriguing as the same behaviour occurs at 70 • C for both modes. On the other hand, Mode-II shear stresses at 110 • C revert to having a much more gradual increase from the centre to the side, with slight fluctuations at 2 and 4 mm.  Figure 7 shows the stress contour cut out of cohesive elements near the crack tip, while Figure 8 plots the stress distribution across the crack tip normalised with maximum stress versus the location from the cohesive centre. Due to the FE half model implemented, the stress values are mirrored across the centre. While the stress contours all show the same characteristic of crack initiation from the centre location towards the sides, a substantial change can be noticed by the way stress distribution evolves towards higher temperatures. At 30 °C, elements concentrated at the centre possess much higher stress whilst elements 4 mm and further have much lower stress, about 25% lower than the centre. This stress localisation drastically changes during a temperature of 70 °C where all elements except the outermost (9 mm) from the centre reach maximum Mode-I stress. With the highest temperature of 110 °C, the stress distribution is still akin to that of 70 °C but the stress drops slightly at a distance of 4 to 8 mm (less than 10%). At higher temperatures, stress progression is more rapidly spread from the centre to the outer elements.  Crack front stress distributions of Mode-II FE models at temperatures of 30, 110 °C are exemplified in Figures 9 and 10. Stress contours from Figure 9 show M typical stress and a crack initiating from the sides of the ENF specimen compa analyses by Koloor and Tamin [11]. The corresponding normalised stress values the crack-initiating front are plotted in Figure 10. At 30 °C, localised normalised are fluctuating more across the crack tip with values of 0.85 to 0.93, while the edge e has significantly higher shear stress. As is with Mode-I, 70 °C Mode-II stress featu trend where all elements from the centremost to an 8 mm distance have the same

Conclusions
This paper utilises finite element cohesive zone modelling to evaluate tempe effects on Mode-I and Mode-II delamination between aluminium and a GFRP lam hybrid. Separate temperature-adjusted material models are applied for the alumi GFRP and cohesive zone to simulate the effects of temperature to full accuracy. Aft idation, localised stresses across the crack tips are plotted, where the variation of distributions leading to crack initiation are compared between each temperature.

•
The validity of the FE model at each temperature is verified with a slope max difference within 5.73% for Mode-I at 110 °C and 7.26% for Mode-II at 70 °C.

•
Crack front stress is concentrated in the middle for Mode-I, while stress is focu the sides of Mode-II delamination. Results for 30 °C are characterised by a mor tuating gradual stress variation for both modes.

Conclusions
This paper utilises finite element cohesive zone modelling to evaluate temperature effects on Mode-I and Mode-II delamination between aluminium and a GFRP laminated hybrid. Separate temperature-adjusted material models are applied for the aluminium, GFRP and cohesive zone to simulate the effects of temperature to full accuracy. After validation, localised stresses across the crack tips are plotted, where the variation of stress distributions leading to crack initiation are compared between each temperature.

•
The validity of the FE model at each temperature is verified with a slope maximum difference within 5.73% for Mode-I at 110 • C and 7.26% for Mode-II at 70 • C. • Crack front stress is concentrated in the middle for Mode-I, while stress is focused on the sides of Mode-II delamination. Results for 30 • C are characterised by a more fluctuating gradual stress variation for both modes.

•
The stress distribution at 70 • C is very polarised, where all elements except the outermost one of Mode-I have practically similar peak stress, while the opposite is observable in Mode-II. A similar trend of more gradual stress variation can be discovered for 110 • C. • DCB and ENF trends employed from experimental tests successfully obtain temperaturedependent cohesive zone properties. • A Johnson-Cook material model with temperature dependency and Chang-Chang material model properties at each temperature ensured proper modelling of specimen bending and flexure.

•
The validated temperature-dependent cohesive zone model demonstrates the applicability of the current methodology to analyse laminates at high temperatures.