Damage Mechanics-Based Failure Prediction of Wirebond in Power Electronic Module

A damage mechanics-based numerical approach for the prediction of the damage evolution in wirebond structures of the power electronic module (PEM) is presented. A simplistic damage evolution model is developed in an in-house finite element code, with a demonstration focused on the analysis of the wirebond damage evolution by thermally induced stresses in PEM subjected to varying thermal loads. The novelty of the proposed methodology is the damage evolution realized at the level of each discretised mesh element of the finite element model of the PEM structure in the numerical approach and the associated impact of damage on the mechanical material properties of that element. A simplified PEM structure is utilised as a case study to demonstrate the proposed damage evolution modelling. The thermal load of each discretised element of the PEM structure was imported from an external thermal code. From the thermally induced stresses, plastic strain rates were approximated and then, using these metrics a damage evolution metric was derived. The damage distribution plot of the wirebond structure for the applied load in the case study indicates that maximum damage accumulation at the heel structure reaches 2.4% of the total damage after 3 seconds. By extrapolating the trendline of damage evolution in wirebond, the time of the structural failure was also predicted. The maximum von Mises stress was observed on the busbar which reaches 64 MPa. The extreme stresses found at the busbar are attributed to the high value of the coefficient of thermal expansion of the busbar material.


I. INTRODUCTION
As illustrated in Figure 1, power electronic modules (PEMs) play a vital role in the conversion, control of alternative energy generation, and distribution.PEM typically uses switching electronic circuits to control the change of voltage/current and frequency level.PEM devices are vital in connecting renewable energy resources such as wind turbines with power grids and the transportation of energy.PEM is an essential technology in all future sustainable energy scenarios.Modern PEM devices are everywhere, for example, in electric trains, motor drives, lighting equipment, etc.The The associate editor coordinating the review of this manuscript and approving it for publication was Ching-Ming Lai .
global power electronics market size, which in 2007 was $9.8 billion, is expected to reach $47 billion by 2027 [1].The percentage of electrical energy controlled by PEMs increased from 40% in 2000 to 80% in 2015 [2].
PEMs are highly inhomogeneous structures assembled with components such as semiconductors, ceramic substrates, copper conductors, aluminium wire/trace, and Printed Circuit Board (PCB) composite materials [3].At the heart of the PEM device are the semiconductor devices (MOSFETs, diodes, IGBTs, etc.).Silicon semiconductors remain widely used and have been proven to be sufficiently reliable.Nevertheless, (a) the limitations of packaging technologies, (b) passive and peripheral components, (c) solder material reliability, and (d) the cost currently limits the allowable junction temperature of the PEM device to 175 • C [4].In automotive drive trains, rail traction systems, aerospace, electric cars, renewable energy generation interfaces, etc., PEMs are subject to significant thermal load and environmental temperature cycling.Modern high-power density PEM devices have an increased temperature and an increased thermal cycling range; both tend to reduce the reliability of PEM devices.Among the new power devices, insulated gate bipolar transistor (IGBT) devices are widely accepted and are increasingly used in traction applications such as locomotives, elevators, trams, and subways.Higher temperature cycling of the modern power device can initiate thermo-mechanical failure mechanisms in the device due to the coefficients of thermal expansion (CTE) mismatch between various materials of the components.Widely known failure modes of modern power electronic modules are (a) wirebond fatigue at the aluminium wire/silicon chip interface, (b) solder fatigue, (c) cracks in silicon chip and substrate, (d) corrosion of interconnections [5], and finally, electromigration induced failure in the interconnects [6].
In order to manage electric currents of up to 75A per chip, it becomes necessary to use multiple wires.Each wire is bonded on the chip using thin aluminium pads or a plate of a few µm thickness.Two major failure modes often occur in wirebond: (a) heel/bond crack and (b) lift-off failure mechanism.These two failure modes occur at the two weakest points (the bond heel and tail (foot)).The lifting of a wirebond can be regarded as a contagion effect because it leads to a non-homogeneous current distribution on the power module chip and subsequently, the higher local and average temperature which accelerates the lifting process of more wires [7].This study addresses, from a numerical analysis perspective only, the wirebond failure due to varying thermal loads acting on a PEM structure.In this work, we focus on a modelling approach where the accumulation of the damage directly influences the material properties of the structure in the analysis.At each numerical iteration time step of the finite element analysis, the accumulation of damage for each mesh element of the modelled material domain decays the relevant mechanical material properties.

II. LITERATURE REVIEW
The two most commonly used lifetime/damage prediction models reported in the literature are (1) models that use predicted damage metric for one cycle of the repetitive load to predict the number of cycles to failure, and (2) models that predict the continuous damage accumulation for any loading profile over a period of time.The former model can normally be used for regular cyclic loading only.The latter model is based on damage mechanics.It predicts damage as a function of time.In damage mechanics-based models, damage parameter D is used as a metric to describe the extent of damage in the structure over time.D is a continuous scaler variable that varies from 0 when there is no damage to 1 when a complete failure occurs.The advantages of damage-based models over lifetime-based models for cyclic loading are that they can predict the process of fatigue damage for arbitrary loading profiles and predict lifetime in the time domain.For example, Lu et al. [8] applied a damage mechanics method to predict crack propagation in IGBT solder joints.By tracking damage evolution in solder joints, crack propagation paths and rates can be calculated.
Fatigue lifetime models that are used to predict the number of cycles to failure for regular cyclic loadings have a power-law relationship that relates lifetime (cycles to failure) to damage metrics such as the deformation range, inelastic strain, and/or plastic work density.To use these models, the damage metric values for the temperature cycle need to be evaluated [9].For irregular cyclic loading conditions, a widely used technique for predicting lifetime is to use a cycle counting algorithm to sort irregular loading profiles into cycles with different loading amplitudes, and then use a fatigue lifetime model for regular cyclic loading and the Miners' linear damage accumulation rule to predict the damage accumulation in the solder layer.This includes the rainflow cycle counting method [10] for counting the cycles in an irregular loading, and then the Coffin Manson or other nonlinear fatigue models [11] for predicting the lifetime of each counted cycle.An application of this methodology to the SnAg solder joint crack development and propagation process was demonstrated by Kostandyan and Sorensen [12] for PEM devices.This method is useful for irregular loading, but it cannot predict lifetime in the time domain.The model does not capture the effects of loading rate changes over time, nor can it take into account the changes in material properties as damage accumulates.Hence, damage-based models that can capture these changes are required in order to predict the reliability of PEM devices when subjected to real environmental temperature loadings.The damage mechanics-based models in the literature are functions of material physical behaviour (accumulated creep strain, plastic work density, etc.).

A. DAMAGE MECHANICS BASED MODELS
Over the years, many damage mechanics models were developed for mostly solder materials since solder interconnects are highly vulnerable to stresses exerted by environmental and system temperature, which resulted in a thermal mismatch caused by the difference in the thermal expansion coefficients of the PEM packaging materials.
• The damage process due to the degradation of the microstructure is irreversible.Under mechanical loading, the extreme plastic strain and the resultant damage can concentrate at various locations in the material.
A damage evolution function proposed by Tang and Basaran [13] for Pb37/Sn63 solder interconnect material of the BGA electronic packages under thermal cyclic loading was evaluated with experimental data.The damage model D was defined in equation ( 1) where e, φ, N 0 , k, T, and m 0 are respectively increments of internal energy, free energy, Avogadro's constant, Boltzmann constant, absolute temperature, and average molecule quantity/mol.
• A damage model that incorporates the cohesive zone constitutive model was proposed by Abdul-Baqi et al. [14].The cohesive zone model is a numerical approach to model the crack initiation and propagation [15] of interconnect materials.The damage evolution rate in Schreurs's damage model is defined in equation ( 2) where k is the initial stiffness of the cohesive zone, D(t j -1) is the previous damage, c is the constant that controls the damage accumulation , r and m are constants which control the decay of the reaction force, is the relative opening of the cohesive zone, and σ is the cohesive zone endurance limit.Schreurs's damage model is supplemented with an evolution law to account for the gradual degradation of the solder material and the resultant accumulation of damage during a mechanical cyclic process.
• A damage model which utilised a creep-plasticity constitutive equation to capture the behaviour of SnAg solder was proposed by Stolkarts et al. [16].Stolkart's damage model is defined in equation ( 3) where k is a material constant, f is a function of stress, strain, and their time derivatives.This damage model can manage the loading with and without dwell times and incorporates various thermal cycling dwell times and ramp rates.For lead solders such as SnPb, the value of k is two [16].
• For thermo-mechanical cyclic behaviour, Lemaitre and Desmorat's creep fatigue model [17], is the most widely used.Based on Lemaitre's creep model, a new continuum damage model was proposed by Xiao et al. [18] as in equation ( 4): where k is the damage exponent which corresponds to the damage evolution rate.Nf is the number of cycles at failure, and D 0 is the initial value of the damage.The parameters of Xiao's damage model were determined by strain measured from the thermo-mechanical cycling experiments.For SnAgCu solder, Xiao's damage model parameters N f and k were 3876 and 0.154, respectively [19].For other solder materials such as Sn3.5Ag, Xiao's damage model parameter values k and N f are not available in the literature.
• Wen et al. [20] proposed a continuum-based damage model which includes the McDowell creep plasticity constitutive equation for lead-free solders.Wen's damage model is defined as a function of the physical damage metric, ω and it holds a power law relationship with the damage metric ω as in equation ( 5) where the parameter D e is the critical damage parameter, and η is the mechanical characteristic of the damaged solder layer.
The damage metric ω is a function of the number of cycles at which cracks are initiated and ω c is the function of the limit value of the number of cycles at which the solder layer is structurally damaged.
• Towashiraporn et al. [21] proposed a continuum damage model for the solder layer under cyclic isothermal mechanical and anisothermal loadings.The damage model in their work was defined in equation ( 6) where ξ c , β and γ are material constants and ξ D is the equivalent inelastic strain.Using the critical damage value (0.85) for the SnPb solder layer, crack advancement is predicted by extrapolating the damage variable using Taylor expansion.Most of the damage mechanics.models discussed above were developed for solder interconnect materials.Damage mechanics-based models often rely on computationally expensive numerical modelling to predict material behaviour, which makes it very time-consuming.Yang et al. [22] proposed a time integration damage model for aluminium wirebond in PEM devices that does not require complex numerical modelling.For aluminium wirebond, Yang's model is among the very few that can be found in the literature to the best knowledge of authors.The next section addresses the methodology of the proposed damage mechanics-based failure prediction of the wirebond structure.

III. METHODOLOGY
This paper presents a theoretical methodology for the application of simplistic damage mechanics-based damage evolution model to the wirebond structure.Since Yang et al. [22] model for wirebond was developed based on the assumption that the wirebond structure is a simplified tri-layer shape structure, Yang's model is not suited to wirebond in PEM structures with complex topology.Therefore, we opted to use a simplistic damage mechanics-based damage evolution model for the wirebond structure of a PEM device in the context of a numerical discretisation perspective.The novelty in this study is that the accumulation of the damage directly influences the material properties of the structure.Hence, at each numerical iteration time step, the accumulation of damage of each element decays the relevant mechanical material properties such as stiffness, as illustrated in Figure 2.

A. UNDERLYING EQUATIONS
The numerical discretisation scheme of the partial differential equation is based on the finite element analysis (FEA) approximation.The damage methodology was implemented in the software code PHYSICA [23].The PHYSICA software code was designed to interact with other modelling codes in a co-design modelling concept.It read two input files: (a) a geometry file with element, node, and boundary condition details; and (b) a text file with analysis instructions called ''inform.txt''.The ''inform'' file consists of material properties, boundary constraints, analysis details, etc.Additionally, each element of the model temperature data can be read by the code.Generally, the temperature data of each element can be generated by other numerical modelling codes.The PHYSICA code consists of mechanical constitutive laws for creep, rate-dependent plasticity, linear elasticity and damage evolution model.The linear elastic equation is defined as where The total overall strain consists of elastic, plastic and thermal strains which are defined as The thermal strain caused by the coefficient of thermal expansion mismatch is given as ε th = α T , where α is the coefficient of thermal expansion and T is the temperature difference between the reference temperature and the structure temperature.ε pl is the viscoplastic strain which is described in section (III-C).The strain energy stored within a volume V is defined as The next section addresses the numerical discretization approach employed in this paper.

B. FINITE ELEMENT DISCRETIZATION
The FEA discretization scheme utilised in the code is to split the entire domain of interest into nodes and elements.The solution of the physics of interest at each node is stored.
As an example, in the case of elasticity, the solution of the displacement is stored.Two types of 3D elements were used in the code, such as the wedge element with six nodes and the brick element with 8 nodes.The values stored at the nodes are interpolated within each element using a shape function where Ni is the weight associated with node i and is a function of the local coordinate position within the element.At nodi i the weight Ni must be equal to 1, and at every point in the element, the sum of all weights must be 1.For 8 nodes, brick element shape functions are defined as The local coordinates of each node are defined in Figure 3.
Given the value of the scalar field φ at each node (φ = φi at node i) of the element, the value of the φ at any location within the element is a linear combination of weight multiplied by the scalar values of each node.
Furthermore, the partial derivative of a scalar field concerning local coordinates can be defined as ∂φ (s, t, u) The derivation of the element stiffness matrix and the loading vector of the system are derived by utilising the shape functions.The next section addresses the plastic/creep models used in the methodology.

C. RATE DEPENDENT PLASTICITY AND CREEP MODELS
The rate-dependent plasticity/creep models, implemented in the multi-physics simulation code PHYSICA, are listed below.In the case study, Section VI, the rate-dependent plasticity, and Hyperbolic sine law creep model were utilized.The solder material of the PEM device has creep deformation characteristics.At elevated temperatures, solder material undergoes plastic deformation even if the applied stress is below the yield stress.Such deformation leads to creep and the creep strain is dependent on not only the applied stress but also on time.A material undergoes a permanent change in shape when the applied stress exceeds the yield stress limit.
• The Perzyna plastic model is defined in the following form.
• Peirce model: which is defined as The parameters are as described in the Perzyna model.
• The hyperbolic sin law creep model is defined as where A -creep rate(1/s), n -stress exponent (dimensionless), σ ref -reference equivalent stress level, Q-activation energy (J/mol) and R is the gas constant.

D. VISCOPLASTIC (VP) LAW
The viscoplastic strain rate tensor is defined as where εVP xy is the creep strain rate or plastic strain rate tensor and the f (σ xy ) is the creep model or plastic model as described in Equations ( 13), ( 14) and (15).In FEA code Equation ( 7) is generalized to provide all 6 components of the strain tensor for each six components of the stress tensor.Hence a viscoplastic law as in equation ( 16) can be converted into a 3D form using equation ( 17) where f (σ eff ) is one of the creep/plasticity flow models, σ eff is the effective stress (von Mises stress) defined as and σ Dev xy is the deviatoric stress tensor defined as The simplest approach to implementing viscoplastic law is the explicit solution approach, where each plastic strain tensor is updated for each time step as where ε VP k+1 xy is the viscoplastic/creep strain tensor at time step k and t k is the time at the step k.The explicit solution approach is susceptible to numerical instability due to the exponential growth of strain with stress.To mitigate the numerical instability and numerical overflow a very small-time step may be required which will increase the computing overhead costs.Another approach is called the implicit solution approach.In the implicit solution approach, the VP rate at the previous time step εVP k xy in equation ( 20) is replaced with VP rate at the current time step εVP k+1 xy

1) IMPLICIT VISCOPLASTIC/CREEP APPROACH
The implicit solution scheme starts by first calculating (a) the displacement and stress over the whole mesh, with an elastic solver, then (b) an approximation of the VP strain rate for each element.This continues until the amount of VP strain rate changes between iterations is smaller than the prescribed tolerance.The element VP rate solver works by assuming that the nodal displacements are fixed and any change in plastic/creep strain causes stress relaxation.The strain is calculated based on the relaxed stress.
The unknown variables are in bold in equations ( 21) and (22).To obtain the new strain rate ε s and the new stress after relaxation σ New a Newton Raphson scheme was utilized.I. Set i = 0 II.The set σ New,i=0 = σ Old III.i = i + 1 IV.First, define F xy to be residual of Equation ( 23) Then Newton Raphson's iterative scheme was employed to better approximate the ε s subject to minimization of the function F (Equation ( 23)

IV. CALCULATE THE STRESS AT THE NEXT TIME STEP
where f σ New in Equation ( 23) is defined as in Equation (17).For Pierce model f σ New is defined as For the hyperbolic sin law creep model (Equation ( 15)) f σ New is defined as where σ new xy is the deviatoric stress component in the x-y plane as in Equation ( 19) and, σ new eff is the effective stress as in Equation (18).One of the challenging aspects of the implicit solution approach is the evaluation of the in Equation (24), Hence, a chain rule was applied which is defined as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.× By differentiating the Pierce model (equation ( 14)) with respect to effective stress (∂σ new eff ) By differentiating the hyperbolic sine law (equation ( 15)) with respect to effective stress (∂σ new eff ) By differentiating the effective stress (equation ( 18)) with respect to the stress component for shear stress/strain component (j = y) for normal stress/strain component (j = x) (32) Deviatoric stress component (Equation ( 19)) differentiation with respect to stress component where t, E,ν,and δ xj are respectively time, Young's modulus, Poisson ratio and Kronecker delta.By substituting equations (30), (31), and (32) in equation ( 28) the partial derivative of F for the strain rate can be evaluated.

V. DAMAGE LAW
Damage is the value taken between 0 and 1 as the material undergoes transformation due to strain and stresses.A damage methodology is implemented in the PHYSICA code to allow materials with non-linear properties, degradation, and cracking to be captured.The methodology is based on the concept that each element of the material is composed of two parts: an undamaged part and a damaged part.In FEA, the entire domain of the PEM device is discretised into elements.Each element of the materials stores an internal variable D in the mesh system, and the D value of each element starts with 0 (completely intact) and gradually increases to 1 (completely damaged) as damage accumulates due to inelastic strain as in Figure 4.
The intact material and the damaged material have different material properties depending on the D value as illustrated in Figure 4.An increment in the D value causes the stiffness of the element to decrease.For each element, material stiffness property (Young's modulus) is interpolated based on their D value, such as  where ϕ is the Young's modulus.Unfortunately, in numerical computation, Young's modulus cannot be reduced to zero since this can cause a singular stiffness matrix with no unique solution.Hence, the lower bound of the damaged element's Young's modulus (E) is set to 1, then E is defined as Assuming the elastic strain is consistent throughout the element, then stress in the intact part is given by.
where [D] is elasticity matrix as defined in equation ( 8), dependent on E intact instead of E, hence stress of the intact element magnifies with E intact value as where σ is the stress tensor as in equation (7).A completely damaged element is equal to the element that no longer exists in the mesh.A better way of dealing with removing the damaged element from the mesh is to remove the associated nodes from the mesh.The damage evolution D is based on the accumulated effective inelastic/creep strain ϕ acc which is defined as where the effective strain rate is given by εeff The following law is used to calculate the damage index D.
where B is the damaged material constant Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

VI. SIMPLIFIED POWER ELECTRONIC STRUCTURE
The numerical damage evolution model as described in Section III was deployed on the simplified PEM specimen structure as in Figure 5.The dimensions of the PEM structure are illustrated in Figure 6.A thermo-mechanical analysis was undertaken to analyse the effect of varying thermal loads and subsequent residual stress in the PEM structure.The simplified PEM structure consists of the typical layout of PEM devices such as wirebond connection, copper trace, silicon chip, SnAg solder interconnect between chip and copper  trace, substrate, copper base plate, and copper busbar.The thickness dimensions of each layer are listed in 1.Ten wirebonds provide the electrical connection between the chip and the copper trace as illustrated in Figure 6.

A. MATERIAL PROPERTIES OF THE COMPONENTS IN PEM
The mechanical material properties of the components are listed in Table 2. Copper material and solder interconnect materials have a high coefficient of thermal expansion (CTE) compared to other components, hence large expansion, and contraction to varying thermal loads.An alumina substrate has the highest Young's modulus, which implies that the substrate is mechanically stiffer than other components.Solder material has temperature-dependent material properties such as Young's modulus and CTE.Young's modulus of the solder material is inversely proportional to the temperature, which implies that if the temperature increases in the structure, then the solder material stiffness decreases.
In contrast, the CTE value is proportional to the temperature.Hence, at elevated temperatures, the CTE of the solder material has an extremely high value in comparison with other materials.Solder materials hyperbolic sine creep law model (Equation (15) parameters are defined in Table 3.The rate-dependent plasticity models are a function of strain rate (or time).The Pierce model is a rate-dependent plasticity model, as described in Equation (14).The copper plates and the aluminium wire are assumed to obey the Pierce plasticity model.The Pierce model coefficients of the copper plates and the aluminium wire are in Table 4.To summarise, the solder, copper plate, and aluminium wire have rate-dependent plasticity properties in the modelling.For the aluminium wire structure, damage parameter B as in Equation ( 39) is assigned as fifteen.The 3-dimensional structural boundary condition was imposed on the model.A vertical movement restriction was imposed on the base of the copper plate.Furthermore, at one corner of the copper base plate, total movement restriction was applied in all directions.

VII. COMPARISON OF PHYSICA AND ANSYS FOR CREEP ANALYSIS
The creep analysis was undertaken in PHYSICA software and the widely used commercial software ANSYS (Mechanical APDL) to compare and validate the accuracy of our implicit creep solution approach presented in Section (III.D.1)).The creep solder interconnect material with creep model parameters as in Section (V.A) was utilised in the PHYSICA software as well as ANSYS software.The rest of the component's material properties of the PEM structure consist of elastic properties.Additionally, the solder material's mechanical properties are temperature dependent.The boundary conditions of the PEM structure were similar in both software.The temperature loading imposed on both models was the cyclic load, with a dwell time of 760 seconds, a ramp time of 140 seconds, lower and upper temperatures of 0 • C, 100 • C respectively, and a total cycle time of 1800 seconds.The cyclic load was imposed on the model for four cycles (7200 seconds).
Figure (7) illustrates the resultant displacement (m) of the PEM structure in the PHYSICA code and the commercial software ANSYS Mechanical APDL.In ANSYS, the solid185 element was used to build the model.The generalised Garofalo implicit creep model is used in ANSYS APDL.It can be concluded that both models have nearly identical values for displacement.Furthermore, the effective stress (N) and the accumulated creep strain at a node located at the corner of the solder interconnect were analysed from PHYSICA software as well as ANSYS software at the time steps and plotted in Figure 8.The PHYSICA software is slightly underpredicting the effective stress (as in Figure 9) and creep strain for the implicit scheme approach in comparison with ANSYS software prediction, hence the accumulated creep strain is slightly higher in Ansys as in Figure 8.

VIII. DAMAGE ANALYSIS AND DISCUSSION
For damage analysis, as noted earlier, the time-dependent temperature profile of each discretised element was generated by a thermal modeller.A time-dependent temperature profile of each element of the discretised model is set as the initial thermal boundary condition of the structural analysis.
Figure 10 illustrates the time-dependent temperature profile for one element from each component.The element's temperature profile in Figure 10 consists of the temperature profile along with the time steps (0.01 second) of elements for each component (chip, solder, copper plates, substrate, and wire).The temperature value is extracted from a thermal modeller.The temperature distribution of the components is not uniform at any time step.The transient analysis was undertaken in the PHYSICA code for three seconds.The post-process results generated through the numerical simulation, such as damage, strain, and stress, were saved in a suitable software plotting format.
Since the damage of the wirebond is the primary interest in this modelling and the parameter of damage was only provided for the aluminium wirebond, the plot of damage accumulation in the wirebond along the time steps is in Figure 11.From the numerical modelling, the critical locations at which damage accumulation is extreme can be identified.Critical damage accumulation in the wirebond was occurring at the heel locations, as illustrated in Figure 13(a).At the end of 0.5 seconds, the maximum damage accumulation at the wirebond is 1.9%, and at the end of 3 seconds, the maximum damage accumulation at the wirebond is 2.4%.From a theoretical perspective, once the damage accumulation reaches 100%, the wirebond structure is completely damaged.The damage accumulation in the wirebond affects the stiffness of the aluminium wire since 25224 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Young's modulus (E) of the aluminium wire decreases as a result of the damage increment as in Equation (34).As a result, the decrease in Young's modulus value accelerates the accumulation of damage in the wirebond.The effective stress (von Mises) distribution of the PEM structure for every one-second time interval is in Figure 12.The copper busbars were subjected to extremely high effective stresses.At the end of the 3-second time simulation, the maximum effective stress on the busbar reaches 62 MPa.One of the reasons for the higher effective stresses on the busbar is the copper material's high CTE value compared to other materials' CTE values.The high CTE value of the copper busbar causes higher contraption and expansion on the free-standing copper structure.

IX. CONCLUSION
This paper reported the development and implementation of a methodology for damage mechanics-based prediction of damage accumulation in the wirebond of a PEM structure in the context of a numerical modelling approach.The methodology consists of finite element numerical discretization of the linear elastic equation for the mechanical boundary constraint and thermal load at each discretized element along with the time steps.The thermal load for each element of the discretised model was generated from a thermal modelling code.With thermal load and the structural boundary condition, a structural analysis was undertaken, and the resulting stresses were utilised to evaluate the viscoplastic strain rate.An implicit viscoplastic solution approach was implemented in the code to approximate the viscoplastic strain rate.
For the wirebond structure, the Pierce plastic model was utilized.Then a damage model was utilised to predict the damage index of each element of the wirebond structure.The novelty in the proposed methodology is that the accumulation of the damage in the material directly influences the material properties of the structure.The above methodology was developed in an in-house software, PHYSICA, and then demonstrated in a case study.
A PEM structure with a simplified layout, consisting of wirebond connections, copper traces, silicon chip, solder interconnect, copper busbar, substrate, and copper base plate, was used in this study.For modelling simplicity, wirebond structures were modelled as rectangular blocks rather than cylindrical solids.The components, such as the copper base plate, aluminium wirebond, and solder interconnect, have rate-dependent plasticity and creep properties at elevated temperatures.The damage plot of the wirebond structure indicates that after 0.5 seconds, maximum damage accumulation at the heel of the wirebond reaches 1.9% and after 3 seconds, maximum damage accumulation at the wire reaches 2.4%.By projecting the trendline of the plot, the wirebond structure will structurally fail for the continued thermal load as in Figure (13b) after 491 seconds.The maximum von Mises stress in the copper busbar is 62 MPa after 3 seconds of temperature loading.The higher effective stresses in the busbar in comparison to other components are due to the higher CTE value of copper materials.
The numerical framework detailed in this paper enables us to account directly for the effect of damage accumulation on the material properties, thus allowing for more accurate and representative analysis and prediction of the material degradation and failure.

FIGURE 1 .
FIGURE 1.The role of power electronic module.

FIGURE 2 .
FIGURE 2. Modelling methodology presented in this study.

FIGURE 4 .
FIGURE 4. The evolution from intact material to damaged material based on damage D.
component (j = x) (31) By differentiating the equation (21) with respect to the current strain rate component ∂σ new xj

FIGURE 8 .
FIGURE 8. (a) Accumulated creep strain on the corner node of the solder interconnect by ANSYS and PHYSICA software, (b) corner node of the solder layer.

FIGURE 9 .
FIGURE 9. Effective stress (N) prediction on the corner of the solder interconnect by ANSYS and PHYSICA software.

FIGURE 10 .
FIGURE 10.The temperature profile of elements in the components of the PEM specimen structure.

FIGURE 13 .
FIGURE 13.(a) The plot of damage accumulation (a scaler value between 0 and 1) at critical location versus time (s), (b) Critical locations of the wirebond at which Damage accumulation are extreme.

TABLE 2 .
Mechanical material properties of the components in the PEM specimen.