Application of augmented finite element and cohesive zone modelling to predict damage evolution in metal matrix composites and aircraft coatings

This paper reports on the application of a special purpose finite element analysis tool that combines augmented finite element methodologies (AFEM) and cohesive zone model (CZM) methods to simulate initiation and propagation of both cohesive and adhesive cracks. The constitutive behaviour of an aluminium/silicon carbide metal matrix composite was predicted and compared with experimental data as an example of a material system controlled by cohesive cracks. The simulation allowed to determine the strain level, at which particle fracture was initiated and illustrates how the overall material response is dominated by particle fracture beyond that strain level. The effects of silica filler particles on the lifetime of polyurethane matrix aircraft coating systems were investigated in a second example in which adhesive cracks at the filler/matrix interface are a dominant failure mechanism. The influence of particle volume fraction and particle/matrix interface adhesion strength on coating lifetime predictions were investigated and the results show that low filler particle volume fraction and high interface adhesion strength improve coating durability. In general, the paper demonstrates the potential of combined AFEM and CZM micromechanical damage simulation to gain improved understanding of damage mechanisms in heterogeneous materials and to support analysis and design of advanced material systems.


INTRODUCTION
Many modern design solutions include the application of advanced materials such as polymer, metal or ceramic matrix composites, because conventional homogeneous materials cannot meet all performance requirements.But many mechanical and material issues remain unresolved for these complex materials, which are naturally in homogeneous and display most of the time anisotropic material behaviour.In addition, they frequently exhibit initial microscopic flaws after manufacturing.
The application of damage tolerant design principles is an essential component for the safe use of these materials.Currently, this relies primarily on experimental *Corresponding author.Email: m.veidt@uq.edu.au.
Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License testing to determine failure loads and damage limits; an approach that is extremely costly and time consuming.In addition, the possibility to evaluate the in-service growth of damage is essential to guarantee the safe performance of safety-critical components during their lifetime.
These are the reasons why substantial efforts were made over the last two decades to develop virtual testing methodologies and tools (LLorca et al., 2011;Yang et al., 2011;De Rosis et al., 2014).The motivation was not only to reduce the need for physical testing, but also to enlarge the design flexibility by providing capabilities for extended parameter studies.Even today, the challenge remains to provide effective and accurate models to describe the complex failure mechanisms in heterogeneous, anisotropic materials.
Finite element (FE) modelling is the most common numerical tool to develop virtual testing methods.Many approaches exist how to treat strain localisations or discontinuities, for example caused by cracks in a FE setting and a lot of research has been done to improve the existing techniques.One promising technique to model weak as well as strong discontinuities in heterogeneous materials is the augmented finite element (AFEM) or phantom node method (PNM), (Fang et al., 2011;Hansbo and Hansbo, 2004;Wells, 2001).It introduces additional degrees of freedom in elements containing discontinuities, which is then modelled as a superposition of two independent continuous elements.The main advantage of the technique compared to other formulations, in particular the widely used extended finite element method (XFEM), is that only standard shape functions are used and elemental locality is preserved.A cohesive zone model (CZM) is applied along the discontinuity surface to incorporate non-linear damage behaviour (Moës et al., 1999;Xie et al., 2006).This paper reports on the application of a special purpose FE implementation that incorporates AFEM and CZM formulations to investigate the damage initiation and evolution in two inhomogeneous materials with very different failure modes.The first application investigates the case of spherical silicon carbide (SiC) particles in an aluminium alloy matrix.In this case the primary failure mode is ductile failure in the matrix followed by brittle particle cracking in the higher strain regime.The interface between particles and matrix is considered to be perfect for this metal matrix composite (MMC) material.The second application explores damage initiation and evolution in aircraft coatings consisting of a polyurethane (PU) matrix filled with silica particles.Debonding between particles and matrix is a major failure mode for this heterogeneous material.Applying the new damage simulation FE tool that combines AFEM and CZM formulation to these two types of materials with these fundamentally different failure modes demonstrates the potential of these damage simulation methodologies to evaluate the damage characteristics of a broad range of materials.
The paper is organised as follows.The basic formulation of the FE method is introduced in the Methodology section including the description of the process to calculate the stress state at the Gauss points of the augmented element containing the tip of the discontinuity.The application of the program to the MMC material is discussed next including a detailed description of the determination of the CZM parameters for the aluminium alloy matrix material from experimental data.The Aircraft Coating Degradation section reports on the results of lifetime predictions of aircraft coating systems and the influence of system parameters such as filler volume fraction and matrix/filler particle interface strength on the lifetime in accelerated degradation testing.Conclusions of the work are presented in the final section.

METHODOLOGY
A special purpose FE analysis tool with AFEM and CZM formulations was implemented in Matlab®.The AFEM formulation follows the approach of Hansbo and Hansbo (2004).Mathematical elements are introduced to account for the cracked sub-domains.Cohesive tractions are applied on mathematical elements to model the plastic process zone.The cohesive tractions follow a material dependant, cohesive constitutive law which describes the relationship between cohesive tractions and the relative displacement of crack surfaces.The CZM formulation follows the method of Discrete CZM (DCZM) reported by Xie et al. (2006).The interface behaviour is modelled by 1-D interface elements.

Implementation of AFEM
Figure 1 shows the flow chart of the FE tool.The FE mesh and boundary conditions are generated with a commercial FE software package (e.g.Ansys) and read into the Matlab ® program automatically through an interface script.After the element stiffness is assembled into the global stiffness matrix, equilibrium equations are solved.Due to the existence of mathematical elements and interface elements, the problem becomes non-linear.A Newton-Raphson iteration method was implemented to solve the non-linear equilibrium equations.Strain and stress are then calculated based on the solved displacement field using standard FE procedures.Next, the failure criterion is examined to evaluate if any discontinuities initiate or propagate in elements and in what direction the discontinuities advance.If the failure criterion is met, the elements are augmented.Then the FE analysis continues to the next loop to calculate and assemble modified element and global stiffness and solve the non-linear equilibrium equations.If there is no new crack initiation or propagation detected, the next load step is initiated until the final load step is applied and no more discontinuity initiation or propagation is detected.Displacement, strain and stress results can be visualized in Matlab ® using special purpose post-processing scripts.The formulation has currently 2-D capability, which enables the evaluation of conservative estimates of damage evolution in many important practical cases of throughthickness micro mechanical damage initiation and propagation such as matrix cracks and delamination in composite laminates, bondline damage in adhesive joints, damage in coatings, etc.

Initiation and propagation criteria
The description of micro-crack dominated failure mechanisms of materials requires the knowledge of three major processes, viz.i) a crack initiation criterion; ii) the calculation of the crack propagation direction; and iii) a cohesive law that describes the relative displacement between crack surfaces.Many criteria were proposed to determine crack initiation and propagation direction, e.g.maximum principal strain criterion, maximum principal stress criterion and maximum strain energy criterion for crack initiation; maximum circumferential (tangential) stress criterion, maximum principal stress criterion and maximum energy release rate criterion for crack propagation direction.Particular criteria are favourable for certain types of materials.In this study the maximum principal strain criterion (Chai and Chiang, 1996;Strawbridge and Evans, 1995) was employed for the MMC as well as the aircraft coating system to determine crack initiation and propagation, that is, a crack will initiate or an existing crack will propagate when the maximum strain at the crack tip exceeds the strain-at-break εb of the material.The maximum principal stress criterion (Maiti and Smith, 1983) was employed for these two types of materials to determine the crack propagation direction.The CZMs used for the two material systems are discussed in the following two sections, respectively.

Nonlocal stress evaluation
The AFEM implementation assumes that newly initiated or propagating cracks create discontinuities in entire elements, that is, the crack tip is always located at the boundary of an element (Figure 2).Therefore the non-local stresses at the discontinuity tip and the corresponding principal stress directions to determine the propagation direction of the discontinuity have to be calculated by nodal averaging of stress tensors as shown in Figure 2.Each Gauss point contributes to the stress values at a node as indicated by the dashed boxes and the arrows in the figure.The non-local average stress is obtained by averaging the stress components from the Gauss points that contribute to the same node in the surrounding elements.For example in the case of Gauss point ②, the average stress is obtained by arithmetic averaging of the stress tensor values from all Gauss points within the corresponding dashed box.The same procedure holds for point ③.The situation for points ① and ④ is slightly more complicated, because they each involve one Gauss point in a mathematical element (ME) indicated by the red "x".Since the ME represents only a physical sub-domain compared to a non-augmented standard element (SE), the contribution from the Gauss point within a ME has to be reduced compared to those located in SE.Hence, a weight factor wf for the Gauss point in a ME was introduced, which equals the ratio of the area of the sub-domain compared to the area of the whole element.Therefore, the averaged stress tensor σ at Gauss point ① is where σ represents the stress tensor at Gauss point ⑤ and i σ are the contributions of the stress tensors at the three Gauss points in the corresponding SEs.This procedure allows the calculation of the principal stress direction at each Gauss point in the crack tip element.
The AFEM methodology has been validated a number of times (Ling et al., 2009;Mergheim et al., 2005;Zhou, 2010).The performance of the special purpose Matlab® tool used in the case studies presented in the following sections was verified through comparison of numerical examples that were used in the validation studies, viz.a) crack propagation in a double-cantilever beam to validate the prediction of crack length; and b) crack propagation in a beam under three point bending load to validate the prediction of crack propagation direction.The results of these verification experiments were presented in Han (2014) and Han and Veidt, (2012).

SIMULATING THE CONSTITUTIVE BEHAVIOUR OF AN AL-2080/SICP MMC
This section reports on the application of the AFEM tool to predict the constitutive behaviour of MMCs.

Material and failure modes
The material used in this study was a 2080 aluminium alloy (3.6% Cu, 1.9% Mg, 0.25% Zr) reinforced with SiC particles.The average particle size was 5 µm.The predicted results were compared with the experimental data taken from Chawla's work (Chawla et al., 1998;Chawla and Shen, 2001).Figure 3 shows the representative tensile stress-strain curves of Al-2080/SiCp composites with varying volume fraction.
In general, different failure modes exist for MMCs depending on the constituents and processing method.Mortensen and Llorca (2010) observed that interface decohesion was the main damage mechanism in the composite material made from a pure aluminium matrix, while the damage mechanism changed to particle failure  when a harder aluminium alloy matrix was used.Ayyar et al. (2008) identified the factors to control the failure modes of Al/SiCp systems.Interface decohesion took place prior to particle fracture in Al/SiCp composites manufactured by liquid infiltration techniques, while damage was mainly controlled by particle fracture in Al/SiCp composites formed through powder metallurgy and extrusion, which generally form a strong mechanical bond at the particle/ matrix interface.The MMC investigated here falls into the second category.Hence, interface decohesion is not considered in this section, and the interface is assumed to be perfectly bonded.The case of matrix/inclusion interface failure is studied in detail in the Aircraft Coating Degradation section for the PU matrix/silica filler particle material system.

Finite element model
The representative volume element (RVE) method introduced by Hill (1963) has been widely used in the materials research community to derive homogenized constitutive laws.A RVE model is the smallest material volume with explicit considerations of microscopic material and geometrical details that represents the macroscopic material characteristics.
The recommended size of a RVE is that the RVE has to be at least five times the largest inclusion (Zhou, 2010).In our case this requires that the edge length of the RVE had to be larger than 25 μm.Therefore, the model shown in Figure 4 was used in this study.
This model had a physical dimension of 28×28 μm 2 , and contained 12 randomly distributed, evenly spaced 5 µm diameter particles (purple areas).This gave a particle volume fraction of 30%.The highest particle volume fraction is selected, since it represents the most challenging configuration for the numerical simulation.Two elastically equivalent homogeneous regions were added (indicated by the red elements) in order to limit the influence of boundary conditions.The size of these regions has been selected based on extended parameter studies to insure that the  damage evolution results are not influenced by them.Linear elastic material properties of the homogenized composite material were allocated to these regions (Table 1).The Young's modulus of the equivalent homogeneous material was determined as the tangent modulus of the linear part of the stress-strain curve for a small stress/ small strain simulation where no damages were initiated in the inhomogeneous model.The model was constrained in the horizontal direction along the left boundary, and displacement controlled load was applied in the same direction along the right boundary.The two lower corner nodes were also constrained in the vertical direction to avoid rigid body motion.The three materials were all treated as linear elastic, isotropic materials.The elastic properties of the material constituents in the RVE and the homogenized boundary layers are listed in Table 1.The plastic behaviour of the Al alloy matrix was accounted for by cohesive tractions applied between crack surfaces.The key to correctly predict the mechanical behaviour of the composite is to determine an accurate description of the cohesive constitutive law for the matrix.This is discussed in the following sub-section.

CZM Parameters for the Al-2080 Matrix
As the first step to predict the behaviour of the MMC, the cohesive properties of the ductile matrix were determined.They were obtained by applying reverse engineering procedures and parametric fitting of the experimental result depicted in the 0% curve in Figure 3.A matrix only FE model with homogeneous material properties was created with the same size and boundary conditions as the composite model.
For ductile fracture, the most important parameters of CZMs are the tensile strength and the fracture energy.However, since we are aiming at predicting the entire stress-strain curve, the shape of the CZM curve also plays an important role for accurately matching the prediction with the experimental stress-strain curve.
The CZM curve determined for the Al-2080 alloy matrix material is shown in Figure 5, which plots the traction force between the crack surfaces as a function of the relative displacement of the surfaces.The curve consists of two parts.The first part has an exponential form and controls the transition from the linear to the plastic part in the stress-strain curve.The second part is a straight line with a small positive slope, which models the strain hardening of the material.The entire plastic behaviour of the material is accounted for by one CZM curve.This CZM was developed from the basic form of CZM proposed by Dugdale (1960) for yielding of thin ideal elastic-plastic steel sheets containing slits.It is noted that the CZM curve used to model the material is a phenomenological representation, and the model parameters are not directly related to any physical processes in the damage zone (Wang, 2010).
The first part of the CZM curve follows the exponential form: The predicted stress-strain curve using the CZM curve stated earlier is as shown in Figure 6.In general, the numerical prediction agrees very well with the experimental data.The experimental data in Figure 4 shows that the elastic limit is reached at a strain level of approximately 0.006, that is, the cohesive zone model is activated in the AFEM simulation tool at that load level.To accurately capture the transition from the non-damaged to the damaged state always requires some compromise between the parameters in the CZM model.This explains why the discrepancy between the experimental and predicted stress-strain curves is larger at that strain level compared to the rest of the curve.But even in the transition regime the relative error is smaller than 4%.

RESULTS AND DISCUSSION
The simulation was divided into two steps.In the first step, the FE model in Figure 4 was used.The particle strength was artificially increased so that only matrix yielding happened without particle cracking.In the second step, a bigger FE model was employed and particle cracking was allowed at a tensile load level determined in step one.This two-step approach had several advantages.First, the use of the small RVE model was computationally efficient to predict the stressstrain curve at lower load levels when particle cracking had not yet occurred.Second, incorporating particle cracking in the small model resulted in predictions that were lower than the experimental data, since cracking of a single, brittle particle released a substantial amount of the totally stored strain energy and made the simulated MMC too soft for the small RVE model.
According to Zhou (2010) a RVE model with 40 µm length and 25 particles was used to include particle cracking as a possible failure mechanism in the simulation (Figure 7).A third reason was that crack merging was not implemented in the current version of the special purpose AFEM formulation.There is no conceptual difficulty to integrate crack merging and branching in AFEM and these methodologies were successfully implemented and studied in (Ling et al., 2009).
The predictions with the small RVE model are as shown in Figure 8.The numerical predictions agree very well with the experimental data when the load is below 1.2%.
A typical strain field at a strain level of 1.2% is as shown in Figure 9.The closed, crack-like features are cohesive zones, which indicate that plastic deformation occurred in these areas.If the crack openings exceeded the critical opening of 1.7 μm, cohesive tractions disappeared and fully-separated cracks formed.For loads higher than 1.2% particle cracking occurs in the real material.Since we artificially increased the strength of the particles this did not happen in the simulation, and the stored strain energy could not be released through particle cracking.This caused the linear increase of the predicted stress values for strains over 1.3%.In addition, since the capability of crack interaction was not implemented in the AFEM formulation as mentioned above, the possibility to release increased strain energy through further matrix cracking was eliminated.The bigger model was employed to investigate how particle cracking affected the predicted constitutive behaviour of the MMC. Figure 10 shows the predicted homogeneous stress at a strain of 2% as a function of the number of cracked particles.The results demonstrate that the predicted stress declined with the number of cracked particles, and after five particles had cracked the stored strain energy was released and the prediction was close to the experimental data.The particles were assumed to be perfectly brittle, that is, no cohesive traction law existed for cracks in particles.Figure 11 shows an example of the strain field when four particles cracked.Note that the cracking particles have been selected carefully to avoid crack interaction.
In summary, the AFEM simulation of the Al-2080/ SiC p MMC demonstrated that the method was able to reproduce the elastic-plastic behaviour of the aluminium matrix with excellent accuracy using a two stage CZM with an initial exponential part modelling the elasticplastic transition and a linear part for larger crack opening displacements modelling the strain hardening of the material.The simulation also showed that particle cracking started at a strain level of approximately 1.2%.It was demonstrated that particle cracking is a major failure mechanism to release stored strain energy.But because crack merging and branching was not implemented in the special purpose AFEM tool the quantitative evaluation of the importance of high-density matrix cracks compared to particle cracking was not possible.The constitutive behaviour of the MMC at high strain levels is a function of the complex interaction of these two major failure modes.

Material and failure modes
External military aircraft coating schemes typically utilize a three-layered coating system, comprised of a chromate pre-treatment, an epoxy primer and a polyurethane topcoat.The topcoat consists of a durable PU polymer filled with silica particles.Figure 12 shows a cross-section of a typical topcoat.Polyurethane provides flexibility to the coating system and holds together all the particles.Silica fillers are added to the coating to increase hardness, change the viscosity, and enhance abrasion resistance on aircraft to counter the impact by airborne debris, such as rain droplets and sand particles (Wicks et al., 2007;Hegedus et al., 1989).Titania pigments (small white dots) provide the primary coloration to the coating.
The durability of the coating system is a primary concern, since environmental factors such as UV radiation, moisture, temperature fluctuations and chemical factors could result in loss of their originally designated mechanical and physical properties, which is known as coating degradation.Fresh coatings have rather good flexibility, whereas after degradation, coatings become more brittle and vulnerable to mechanical loads experienced by the aircraft during flight (Tiong and Clark, 2010).Micro-cracks may form and grow in coatings and accelerate the degradation process.These cracks can allow moisture to come into contact with the substrate and subsequently initiate corrosion.

Finite element model
The AFEM formulation reported in this paper enables the simulation of the evolution of micro-cracks in degraded coating system and hence can be used to predict coating lifetime.Two possible types of cracks may exist in the coatings-cohesive cracks in polymer matrix and adhesive cracks at the matrix/filler interface.These two types of cracks can intersect.In this study, the nano-scale pre-treatment layer was assumed to not affect the crack behaviour on the substrate side of the coating or damage evolution in the micro scale coating system.Thus this layer was not included in the model.The matrix of the primer is a brittle epoxy based polymer.Tensile tests by research partners within the Defence Materials Technology Centre have shown that the average elongation-at-break value is in the order of 4%.Thus under the deformation controlled loading of 15% tensile strain used in these studies, the primer always cracked prior to the topcoat.Therefore, the study focused on the topcoat, which is the main barrier against high strain loading.A typical thickness of 37.5 μm   was taken for the topcoat.Part of the homogenized primer layer was also included in the FE model to provide accurate stiffness in the vertical direction at the topcoat/ primer interface.The effect of the Titania pigments on the matrix properties has been included in the degraded properties of the coating, so the pigments were not modelled.Plane strain conditions were imposed on the PU matrix and primer; plane stress conditions with thickness of 10 μm were imposed on the silica fillers.Vertical constraints were applied on the two lower corners of the model to avoid rigid body motion.Figure 13 shows details of the geometry.Table 2 lists the material properties and the benchmark values of the filler volume fraction, filler diameter, filler/matrix interface fracture strength and fracture energy used in the following analyses.

Model assumptions
The silica filler particles play an important role defining the coating properties, and they have a major effect on damage evolution during coating degradation.The AFEM tool was used to study the effects of filler particles on the lifetime of coating systems.For these fundamental investigations idealised round and elliptical filler particles were used and these filler particles were evenly spaced in the models.The two particular aspects studied in this paper are the effect of filler volume fraction and interface adhesion on the predicted lifetime of the coating system.
The following points related to the simulation studies have to be noted.a) Lifetime prediction: Coating systems are considered to have failed, if cracks propagate through the thickness of the topcoat.Under such conditions, the metal substrate is exposed to the environment and subject to corrosion.b) Mechanical loads: Tiong and Clark (2010) have shown that strains in paint coatings around fastener heads and sealant seams may reach up to 15% during service.This value was used in this study to determine the displacement controlled mechanical load.Hence, the predicted lifetimes presented below apply to coatings at high strain areas and therefore represent worst case estimations.c) Degraded material properties for PU matrix: Mechanical property changes during degradation of the PU matrix were taken from (Skaja et al., 2006).In their studies, unfilled coating samples were exposed to cycles of four hours light (0.55 W/m 2 at 340 nm) at 60°C and four hours dark with water spray at 25°C according to ASTM D 4587-91.The tensile properties were measured after each exposure increment at room temperature.Figure 14 shows the changes of tensile modulus and tensile strainat-break during these accelerated exposure experiments and the interpolation functions for the data used in the simulation studies.It is generally observed that coatings become more brittle upon exposure.The tensile modulus increases almost linearly during the initial 4 weeks, and remains relatively unchanged afterwards.The tensile strain-at-break value decreases dramatically within the first two weeks from approx.100% (the exact initial value is not given in (Skaja et al., 2006) to only 8%, followed by a moderate drop to 3.5% over the next 4 weeks.Reduced elongation-at-break is a common observation for polymers subjected to weathering, where brittle fracture initiates from flaws generated by photo oxidation (Skaja et al., 2006).Because of this embrittlement of the PU matrix during degradation CZM for the matrix material was not considered in the simulation.For each simulation case, material properties at different degradation levels were selected to check if cracks can grow through the entire topcoat layer.If cracks go through and the last cracked element had the similar strain level as the PU's degraded strain-at-break value, then the exposure time at the trial degradation level was deemed as the lifetime for the simulation case.Since the lifetime prediction was based on data determined from accelerated exposure experiments, an "acceleration shift factor" would be needed to determine the corresponding lifetime under natural exposure (Simms, 1987).d) Crack initiation locations: There were two types of possible crack initiation locations considered.The first was to initiate from flaws generated by photo oxidation at the coating surface since the surface suffers from earlier degradation than the bulk (Skaja et al., 2006).This type also includes mechanical damage such as scratches, and repair or refinishing that can introduce non-uniform external surface features which can act as crack initiation points (Jaya et al., 2012).The second possible crack initiation sites were at high strain locations inside the coating, which mostly occur at interfaces between matrix and filler particles.Therefore, a crack starting point was

Effect of filler volume fraction
Research partners from DSTO calculated the volume fraction of silica particles in typical topcoats based on 3-D scanning electron microscopy image analysis and found it to be around 40% (Trueman et al., 2013).In the current studies on the effect of filler volume fraction three values were considered, namely 30, 40and 50%.For all other material characteristics the benchmark values given in Table 2 were used in the simulations.Figure 12 shows that the largest filler particles had a diagonallength between 10 to 15 μm.In this study elliptical particles with long axis of 12 μm and short axis of 8 μm were used.Figure 15 shows the mesh of each model.The average element size was 1 μm.
Figure 16 shows that high volume fraction of filler particles significantly decreases coating lifetime.This result is not surprising because the filler particles decrease the flexibility of the topcoat due to higher packing density and reduced matrix material between filler particles, which results in higher strain level in those areas.Figure 17 shows the cracked topcoat for the 40% filler volume case.The surface flaw developed into a major through-thickness crack by interacting with other matrix micro-cracks and matrix/filler interface delamination.

Effect of matrix/filler interface strength
A major challenge related to crack growth modelling is the determination of the parameters and shape of the CZM that controls crack behaviour.Usually extended experimental work is needed.The primary aim of these  fundamental simulation studies was to compare cases with different filler/ matrix adhesion characteristics.Therefore a reasonable estimation of the range of parameters was adequate, especially since accurate values were not available in particular for coatings at different stages of environmental degradation.The three major parameters in a CZM are fracture energy, fracture strength and separation displacement.Fracture energy and fracture strength are considered the two most important parameters (Borst et al., 2004;Tvergaard and Hutchinson, 1992).As discussed in the context of the simulation studies of the MMC, the shape of the CZM also plays a vital role on crack behaviour, especially for brittle materials (Borst et al., 2004).As the exact shape of the curve was unknown, fracture behaviour was modelled based on the frequently used bilinear CZM curve shown in Figure 18.The bilinear form of CZM was also employed by Yeo et al. (2005) to model demoulding of cured polymer from a rigid polyester mould, which has some similarity to our work.Geubelle and Baylor (1998) also used this form to simulate delamination in composites.
Round, 9.8 μm diameter particles were used to investigate the effect of fracture strength on coating lifetime.Sakamoto et al. (2007) measured interface strength values of 2 to 15 MPa between polyurethane and titanium.Based on these values three cases of fracture strength were considered in our studies, namely, 1, 20, and 50 MPa, respectively.The fracture strength of 20 MPa was the benchmark value used in previous simulation studies as indicated in Table 2. Figure 19 shows that coating durability is improved for matrix/ filler interfaces with high fracture strength.The 20 MPa case does not show a distinct improvement compared to the 1 MPa case.The reason could be that better adhesion not only prevents delamination but also helps to distribute more homogeneously the strain/ stress field in the matrix, which alleviates local strain/ stress concentrations so that more cohesive cracks initiate and grow.This in consequence helps to release stored strain energy.However, in coatings with low interface adhesion there are less cohesive cracks to release the strain energy which will drive the existing cracks to propagate for longer distances.The adhesion in the second case was probably not big enough to fundamentally affect crack behaviour in the matrix.

Sensitivity study
To provide confidence in the results of simulation studies for complex material systems such as aircraft coatings it is essential to perform sensitivity studies to determine how sensitive the simulated lifetime predictions are to changes in the simulation model.Hence filler particle distributions were varied in our study to determine confidence intervals of the calculated lifetime predictions.The two cases of round and elliptical filler particles were used to establish confidence intervals.For each case three different particle distributions were considered.Particle positions in each model were selected randomly, but the average separation between particles was maintained.Figure 20 shows the FE models of the six cases.
The results of the sensitivity studies are presented in Figure 21.The confidence i ntervals were 0.08 weeks for coatings with round particles and 0.12 weeks for coatings with elliptical particles.This corresponds to an uncertainty in the order of 5%, considering that coating lifetimes of approximately 2 weeks were predicted, which is a very satisfactory result for simulating such complex material systems.The slightly larger confidence interval for elliptical particles is probably due to the influence of the random orientation of the elliptical particles, which fundamentally changed the strain fields in the three cases.Coatings with round particles were only affected by the particle distribution.
In summary, the AFEM simulation of the PU/silica particle aircraft coating system demonstrated the capability of the method to simulate the complex interaction of matrix cracks and matrix/ particle debonding.Studies of the effects of filler particles on coating lifetime indicated that low volume particle fraction and high interface adhesion improve coating durability.It has to be noted that the accelerated exposure data used as input data was not specifically for aircraft coatings.It is expected that the lifetime of real aircraft coatings is longer than the current predictions considering that aircraft coating systems are some of the highest performing coating systems available.

Conclusion
This paper reported on the development and application of a special purpose FE simulation method to investigate initiation and evolution of micro cracks in aluminium/ silicon carbide MMCs and polyurethane/silica filler aircraft coating systems.The analysis tool combined AFEM and CZM methods to model initiation and propagation of both cohesive and adhesive cracks.
Results indicated that combined AFEM/CZM simulations enable better understanding of fundamental failure mechanisms in complex, heterogeneous material systems and the roles of various material characteristics on the rate of damage evolution.The studies of these two very different material systems demonstrated the potential of AFEM/CZM micromechanical simulations as design tools to predict and optimize the performance of future materials.
The studies clearly showed that the most challenging part of micromechanical damage simulation is to obtain reliable material properties especially properties at material interfaces, which are required to determine CZM parameters.Extended experimental work is required to fill this gap for any new material system.

Figure 2 .
Figure 2. Illustration of nodal averaging of stress tensors.

Figure 4 .
Figure 4. Mesh and particle distribution in a typical RVE used in the MMC study with 12 particles.

Figure 5 .
Figure 5. CZM for the Al-2028 matrix material with an initial exponential part and a second linear part.

Figure 6 .
Figure 6.Comparison of numerical prediction and experimental data for Al-2080 matrix material.

Figure 8 .
Figure 8.Comparison of experimental and predicted data with small model.

Figure 10 .
Figure 10.Change of predicted homogeneous stress with number of cracked particles at 2%.

Figure 11 .
Figure 11.Strain field at strain of 2% with cracked particles.

Figure 12 .
Figure 12.SEM imaging of cross-section of a topcoat: light grey -silica filler particles (image courtesy of Australian Defence Science Technology Organisation

Figure 13 .
Figure 13.Illustration of geometry and boundary conditions of coating model.

Figure 16 .
Figure 16.Effect of filler volume fraction on coating lifetime.

Figure 17 .
Figure 17.Cracked topcoat of the second case.

Figure 19 .
Figure 19.Effect of fracture strength at matrix/filler interface on coating's life.

Table 1 .
Linear material properties of the composite constituent materials.

Table 2 .
Material properties and benchmark values for coating studies.