High temperature oxidation in polymers is a complex phenomenon, driven by the coupled di usion-reaction process, causing changes in the amorphous network structure and result- ing in property degradation. Prolonged oxidation in polymers results in the formation of a coarse, oxide layer on the outer surface and induces spontaneous cracking inside the material. In this paper, we present a chemical reaction-driven evolving network theory coupled with phase- eld fracture to describe the e ect of oxidation in polymers across di erent length scales. The theory considers the coupling between oxygen di usion, chemical reactions, large deformation of polymers, and phase- eld fracture in a thermodynamically consistent way. Guided by the statistical mechanics, the network theory has been introduced to model the reaction induced chain scissions and crosslinking events causing signi cant changes in the three-dimensional network structure. Further, these microscale events have been consid- ered as the reason behind macroscopic mechanical property degradation, namely oxidative embrittlement. Finally the network theory is coupled with a phase- eld fracture model to capture the macroscale damage and fracture in the polymer under stress-coupled oxidation conditions. We derive the speci c constitutive forms for all the physical-chemical processes based on the thermodynamic inequality conditions, and numerically implement the theory in nite elements by writing ABAQUS user-de ned element (UEL) subroutine. To present the model's capability, numerical examples with standard fracture geometries have been studied. Corresponding author, email: tsain@mtu.edu Preprint submitted to Elsevier October 15, 2020 The simulation results have demonstrated the model's capability of predicting the e ect of oxidative aging on the polymer's response. Keywords: thermodynamically consistent, phase- eld fracture, statistical mechanics, transient network theory, large deformation, oxidative aging, chemo-mechanically coupled, chemical reactions, amorphous polymers, oxygen di usion 1. Introduction High temperature oxidative aging in polymers and polymer matrix composites (PMCs) involves a multi-scale damage and fracture process that starts at a molecular level with multiple chain reactions, leading to the microscale chain-scisson and oxidative crosslinking. These microscale events eventually degrade the material properties at the macroscale and adversely e ect the mechanical response of the polymers [45, 65, 68, 69, 91]. Oxidation in polymers at high temperature is characterized as a coupled di usion-reaction phenomenon driven by the di usion of atmospheric oxygen. The damage induced by oxidation signif- icantly alters the constitutive response of the polymer and notably reduces the life-time expectancy. Polymeric materials are greatly desirable in many applications for excellent ex- hibition of competing properties, (such as light weight, toughness, viscous dissipation, etc.). However, the current literature is still limited in terms of the accurate response predictions of these materials at high temperature oxidative environment, such as, in the components of supersonic jets, pipelines and chemical storage etc. [30]. The experimental studies re- garding the high temperature oxidation in polymers can be categorized into two groups:- i) microscale characterization- mainly focuses on quantifying the bond formation/scissions and the oxidative products through Fourier or Raman spectroscopy [21, 72, 73, 85], ii) macroscale characterization- those reported property changes in polymers and PMCs through combined macro and micrographic observations. Studies were done to measure the oxide layer thick- ness and to estimate an increase in the modulus, decrease in the glass transition tempera- ture, shrinkage strain development and decrease in the failure strain, via indentation test, dynamic mechanical analysis (DMA), tensile and bending tests combined with SEM images 2 [30, 33, 37, 39, 40, 47, 48, 64, 68, 69, 91]. Also, identi cation of the other inevitable e ects of thermal oxidation in case of PMCs as the matrix microcracking and debonding at the ber-matrix interfaces, are reported in [16, 45, 47, 64, 76]. These experimental studies in micro and macro scales reveal that, during oxidation, properties of the polymers experience signi cant alteration, caused due to a competing chain scission and crosslinking events at the microscale [19, 21, 26, 27]. These microscale events, as driven by the closed loop chain reac- tion, interact with the macroscale deformation and stress resulting in damage and fracture of the polymer specimens. Modeling of thermo-oxidative aging is an integral part of the study. While experiments are inevitable to characterize the property changes in the course of oxidation, it may some- times render inadequate to gather all the information regarding material behavior at all the possible scenarios. In reality, the e ect of high temperature oxidation in polymer is a slow process requiring days or months to demonstrate measurable changes in polymer response. However, the experiments are commonly being conducted at accelerated aging conditions [21]. Therefore, modeling the oxidation in polymers has always been a strong focus of the polymer literature. The most popular model is the mechanistic model, which was developed by a thorough investigation of the oxidation reaction kinetics [11]. The property degradation can be linked to the reaction products, which can be calculated based on the mass balance of each individual reactants and products [39, 40]. While using the mass balance, it is im- portant to note that each reactant takes part in the multiple closed loop chain reactions, as described in [9{11, 25{31, 73, 89]. An extension of this model can also be found in the literature, known as the three zone oxidation model, that emphasized in the heterogeneity of oxygen di usion leading to an oxide zone formation followed by an active reaction zone and then an unoxidized core [70, 81, 82]. This heterogeneity in the degree of oxidation results in non-uniform mechanical properties and inhomogeneous stress/strain distribution within an oxidized polymer [46, 47, 81]. In the similar trend, a di usion-limited oxidation (DLO) model is also introduced to predict oxidation within a thick polymer sample [20, 21]. The DLO predictions showed that when the rate of reaction is greater than that of oxygen di u- 3 sion, the reaction gets contained within the surface instead of spreading toward the sample core [20, 21, 71]. Recently, a micro-mechanical model based on the competition between chain-scission and crosslinking events occurring at the polymer network during oxidation was reported to predict the changes in the constitutive behavior of polymers by [66]. On the other hand, a thermodynamically consistent continuum model to predict the constitutive response based on high-temperature oxidation behavior was found in the work of [46, 47]. However, these models do not consider the connection between the micro-scale polymer net- work evolution with the macroscopic fracture and property degradation in a stress-coupled situation, during oxidative aging. Thermo-oxidative aging process leads to a gradual change in the material's con guration at the microscale, before any visible macro cracks or other fracture zones appear. The main contributing factor to this microscale alteration is the evolution of the original polymer network with scissoned and newly- formed crosslinked networks. Describing these changes at microscale through the statistical mechanics description of polymer chains, seems to be more tting as this approach enables a more detailed connection with the statistics of the polymer chains, [34, 57, 80, 83, 86, 88, 90]. Assuming that the statistical distribution of the polymer network can be altered through the chain scisson/crosslinking events, one can readily track the e ect of oxidation via the evolution of a so-called chain distribution function. Further considering an ane deformation of the polymer network, one can establish a direct link between the micro and macro scale kinematics, as done in the entropic elasticity theories. In literature, a number of such statistical models are found, for example, the eight chain model ([8]), the fully network model ([90]), the transient network theory ([88]), a bond deformation energy based model ([59]), to name a few. In order to fully understand the behavior of the polymers undergoing oxidation and their lifetime expectancy during application, a suitable model to capture the damage and fracture across di erent length scales is indispensable and still lacking in the literature. The macroscale damage and fracture in polymers undergoing oxidation is the result of a highly nonlinear, stress-coupled, di usion-reaction process. As the heterogeneous di usion-reaction 4 process coupled with mechanical deformation leads to an inhomogeneous damage propagation within the polymer, a phase- eld fracture model seems to be better suited to predict this oxidative aging process at the macroscale. Phase- eld fracture model has been successfully incorporated to solve various multi-physics problems in the literature [6, 12, 58, 60, 62, 67, 79]. The model is also bene cial to implement numerically within a variational setting of the proposed problem compared to a phenomenological discrete crack or cohesive zone model [14, 15, 42, 61, 63]. An excellent review on the phase- eld modeling of brittle fracture can be found in [2]. In order to overcome the diculties in classical Grith-type theory of brittle fracture, phase- eld fracture was originally derived from the variational methods based on energy minimization as suggested by [14, 17, 32, 42]. The key feature of the phase- eld fracture lies in the approximation of a sharp crack topology in the solid by regularized di usive crack zones governed by a scalar auxiliary variable0d0. This variable can be considered as a phase- eld order parameter which interpolates between the undamaged ( d= 0) and a completely broken state of the material ( d= 1). Further, introducing a crack surface density function, the dissipative fracture energy can be quanti ed in terms of the bulk (volumetric) quantity and be combined with the stored energy of the solids for minimization [13, 14, 61, 63]. In this paper, we develop an oxidation reaction-informed evolving network theory cou- pled with phase- eld fracture to connect the microscale network evolution with macroscopic fracture occurring in polymers during oxidative aging. We use the dynamic or transient net- work theory based on the statistical mechanics framework of the polymer chains, (as used in [74, 80, 83, 86, 88]), to model the microscale network evolution yielded by the chemical reactions. According to this theory, we consider the evolution of the amorphous network in the polymer due to oxidation. The degree of crystallinity in a semicrystalline polymer also gets a ected by the in uence of high temperature but will not be considered in the present micromechanics. The con guration of a polymer network changes as chain scission and new crosslinked bonds are formed. To incorporate these microscale events, the proposed model considers the statistical distribution of the polymer chain con guration, assuming the end-to-end vector of a single chain as a primary kinematical quantity. We assume that this 5 distribution can evolve both as a function of mechanical deformation and oxidation reactions. The evolution of the chain statistics further a ects the stored energy of the polymer system. We link this microscale chain distribution to the continuum kinematics description through the de nition of a chain distribution tensor, as introduced in [88]. Further, the stored energy of the network is written in terms of this distribution tensor [74, 80, 88]. To connect this microscale mechanics with the macroscale damage and fracture under a stress-coupled oxida- tive situation, we use the phase- eld fracture to model the macroscopic damage initiation, crack propagation and ultimate fracture. The framework is developed in a thermodynami- cally consistent way to model the coupling between various physical-chemical processes and to derive the speci c constitutive forms for each of them, that prevails the oxidative aging in polymers. The paper is organized as follows- rstly, we explain the chemistry of oxidation reactions in polymers followed by the statistical fundamental of the polymer chain network used in the standard network theory and how it would be modi ed in the present setting to incor- porate the reaction-driven evolution. Considering the reaction modi ed polymer network, subsequently, we develop a continuum-level chemo-mechanically coupled theory incorporat- ing phase eld fracture to model the macroscale oxidative damage. The speci c constitutive forms are then derived for each individual process. We implement the proposed theory in a nite element setting by writing a user element subroutine (UEL) in ABAQUS [1] and present several numerical simulations to demonstrate the capability of the model. At the end, the concluding remarks are presented. 2. A reaction-dependent network theory coupled with phase- eld fracture to model oxidative aging in polymers As mentioned earlier, during oxidation, polymer undergoes a series of chain reactions that leads to random chain scission and oxidative cross-linking at the microscale (as shown in Fig 1b). These microscale events further interact with the polymer's general state of stress and result into macroscopic property degradation and fracture. In the present section 6 we describe a micromechanically informed, chemo-mechanically coupled theory to connect the microscale modi cations in the polymer network due to oxidation with the macroscopic fracture under mechanical loading. 2.1. Chemistry of oxidation:- the reaction kinetics To understand the multiple chain reactions involved in the case of a polymer undergoing oxidation, we brie y summarize the underlying reaction kinetics. The mechanistic reaction schemes proposed to describe the thermal oxidation in polymers had a consensus in the liter- ature except for some minor variation at the initiation stage [9{11, 25{31, 73, 89]. According to this standard mechanistic scheme, during oxidation, the polymer undergoes six sets of closed-loop chain reactions involving three major steps identi ed as initiation, propagation, and termination, respectively, as explained below: (Ia) Initiation POOH !2P (Ib) Initiation PH + O 2!P+ POOH (II) Propagation P+ O 2!PO 2 (III) Propagation PO 2+ PH!POOH + P (IV) Termination P+ P!inactive products (V) Termination P+ PO 2!inactive products (VI) Termination PO 2+ PO 2!inactive products + O 2(1) where PH represents the polymer substrate; Pis the alkyl radicals; POOH is the hy- droperoxide, PO 2is the peroxy radical and O 2represents oxygen. These set of reactions were originally proposed by [31] and then followed in the polymer chemistry literature, [9, 28, 29, 31], to describe the chemistry of the oxidation process. Incorporating the reaction kinetics in a coupled di usion-reaction-deformation framework a theory has been presented in our recent work [52], to model oxidation coupled large deformation constitutive behavior in polymers. 7 2.1.1. Extent of reaction To connect the reaction kinetics occurring at the molecular level with the micro-scale events such as chain scissions and crosslinkng, we introduce a scalar quantity, as extent of reaction , to describe each individual reaction [52]. These extent of reactions are a measure of aggregates for each reaction happening during oxidation. We de ne the extent of reaction for each reaction ( n) at every continuum point Xin a local dimensionless form as, 0n(X;t)1 (2) where0n0is the number of reaction and varies between 1 to 6 depending on the oxidation kinetics for a speci c polymer. Since polymer oxidation is a multi-component system involv- ing multiple reactions (as described earlier), the chemical concentration for any component would change depending on the extent of all the reactions in which it is participating, as, r =X n Rn _n (3) wherer is the rate of production or consumption of chemical species (measured in moles per volume) andRn is the stoichiometric coecient in reaction n. 2.2. A reaction-dependent network theory based on the statistical mechanics of polymer chains At the outset, we consider most of the polymer chains are attached to each other in an amorphous network via cross-links, representing a 3D solid polymer at its virgin state. We also assume that, the actively crosslinked chains in the network essentially provides the structural integrity for the polymer. In addition, the network also contains some dangling chains representing the unstable hydroperoxides, which could act as the initiator for the oxidation reaction. Hence, the total number of chains in a polymer network ( CT) is the sum of the active and inactive chains. Following [50], under the assumption that individual chains do not di use within the polymer, the total concentration follows the standard conservation equation as, _CT+CTtr(L) = 0 (4) 8 (a) (b) Figure 1: (a) Typical virgin polymer chains with Kuhn segment 'b' and end to end distance vector r, (b) polymer undergoing oxidation developing new cross-links and broken chains. withLrepresenting the velocity gradient as L=_FF1, where Fis the deformation gradient. During oxidation, some polymer chains get broken (via reactions (I) and (III)). Since polymer radicals are highly reactive, these broken chains can get re-attached via the termination reactions (IV, V and VI). These new attachments or bonds are formed with a much lower energy compared to the original bonds in the virgin networks as reported in [66]. Thus, one can consider that, the nal oxidative product has more cross-linked network of smaller chains compared to the virgin polymer with longer cross-linked network. It is important 9 to mention that, oxidation reaction gives rise to a di erent type of network in the oxide zone of the polymer, a "PO-OP" type network compared to the original P-P network. We attempt to describe the polymer oxidation phenomena as a collective response of the chain detachments and attachments at the microscale governed by the chain reactions. With the help of a chain distribution function based on the chain end-to-end vector r, (motivated by the transient network theory as introduced in [80, 86{88]), we use a statistical framework to describe the physical state of an oxidized polymer network. Further, we incorporate the reaction kinetics to formulate an evolution equation for this chain distribution. Finally, introducing the second moment of this distribution, termed as ' chain distribution tensor ' the macroscopic free energy of the polymer can be derived. It is important to note that one of the key limitation of this theory while deriving the free energy is the assumption of Gaussian statistics of the polymer chains, which had been circumvented by [87], with the use of so-called mean eld approximation. In the present study we also assume that this approximation is valid and the rate of attachment/detachment events due to oxidation reactions are independent of mechanical deformation. 2.2.1. Chain distribution function To begin with, we consider the end-to-end vector ( r) of a single chain as a primary kine- matical quantity to describe the polymer chain con guration. Let us consider that every single chain is comprised of N number of freely joined Kuhn segments of length b as demon- strated in Fig. 1a. As postulated in the classical random walk model, we assume that all the actively crosslinked chains are of equal length. Hence, the physical state of the poly- mer can be represented by the statistical distribution of the end-to-end vectors ( (r)) for the entire chain population. Let us consider that, the distribution function (r) can be ex- pressed by a simple relationship between the probability density function P(r) and active chain concentration c( t) as, (r;t) =c(t)P(r;t) (5) wherec(t) has the unit of number of chains per unit of current volume. The active chain concentration is the integral of all the possible chain distributions , (r;t) and can be written 10 as, c(t) =h i (6) where the operator hiindicates the integral over all the possible chain con guration as, hi=Z2 0Z 0Z1 0r2dr sindd! (7) whereand!are the direction of the end-to-end vector rin the spherical coordinates. Further assuming that, at a length-scale much larger than the individual chains, the network is random and isotropic, the stress-free chain distribution can be expressed as a normal distribution P0(r) with a zero mean and a standard deviationp N=3bin each of the three spatial directions. We use the expression from [83] for P0(r) as, P0(r) =3 2Nb23=2 exp 3jrj2 2 (8) which leads to the chain distribution at stress-free con guration as, 0(r;t) =c(t0)P0(r) (9) wherec(t0) is the active chain concentration at the stress-free con guration at time t=t0. 2.3. Evolution of chain distribution due to the combined e ect of oxidation and mechanical deformation The hypothesis is that in a network of polymer chains, the statistical distribution of rchanges as a function of mechanical deformation and chemical reactions. As both the deformation and the oxidative reactions cause changes in the chain distribution (r;t) and active chain concentrations c(t) within an amorphous network, we can write, _ (r;t) =_ (r;t) ox+_ (r;t) F(10) The rst term represents the change in the distribution due to mechanical deformation and the second one considers the change due to chemical reactions holding the deformation con- stant. The subscripts ox and Fdenote oxidation and deformation gradient, respectively. 11 Assuming an ane deformation such that the end-to-end vectors of all chains evolve propor- tionally with respect to the macroscopic deformation gradient [88], we can write the 1st term of Eq. 10 as, _ (r;t) ox=L: (r r) (11) where Lis the macroscopic velocity gradient. Following the de nition of spatial time derivative and the use of transport theorem over the chain distribution space, the derivation of the nal form of Eq. 11 is not new and can be found in literature [88]. In this work, we investigate the 2nd term in further details, which incorporates the change in the chain distribution arising from the oxidation reactions. As we can see from Eq. 1, chain dissociation or scissions occur during reaction (Ia), (Ib) and (III) and chain reattachment happens in reaction (IV), (V) and (VI). We assume that, the chain scission is proportional to the current distribution function (r;t) and depends on the scission reactions. We further assume that, chain attachment happens in a stress-free con guration and dependent on the termination reactions. Also, only a fraction of chains that has been scissioned before gets reattached in the termination process. Hence, one can write a kinetic form to accommodate the change in the chain distribution due to the reactions as, _ (r;t)jF=fa(n)cd(t)P0(r)fd(m) (r;t);withm= 1a;1b;3 andn= 4;5;6 (12) wherecd(t) equals the total number of detached or scissioned chains and faandfdare the reaction dependent functions. Further m;n denotes the reaction number, and P0(r) is the probability of chain con guration in the stress-free state. The number of detached chains can be found as, cd(t) =CTc(t). In addition, we assume that, the reaction dependent functions can be expressed as follows: fd(m) =kd(1a+1b+3) fa(n) =ka(4+5+6)(13) 12 wherekdandkaare rate-like reaction dependent coecients, that denote the chain de- tachment and attachment rates, respectively. Substituting Eq. 13 into Eq. 12, we can get the evolution of the chain distribution due to the reaction as, _ (r;t)jF=ka(4+5+6)(CTc(t))P0(r)kd(1a+1b+3) (r;t) (14) Substituting Eq. 11 and 14 into Eq. 10 yields the complete evolution equation for the chain distribution as, _ (r;t) =L: (r r) +ka(4+5+6)(CTc(t))P0(r)kd(1a+1b+3) (r;t) (15) 2.3.1. Evolution of the concentration of active chains Following Eq. 6, one can see that, the active chain concentration is the integral of all the possible chain distribution function, (r;t). As the chain distribution evolves, the concen- tration of the active chains would also change with time. Using the fact, _ c=_h iand by integrating Eq. 15 over the chain space, we can write the evolution of the chain concentration as, _c=ka(4+5+6)(CTc(t))kd(1a+1b+3)c(t)hP(r)ic(t) tr(L) (16) wherehP(r)irepresents the integration of the chain distribution over all the con guration and value ofhP0(r)i= 1 is used in the 1st term. As per this expression, the active chain concentration would increase while the reaction kinetics are dominated by the termination reactions and would decrease when the kinetics are dominated by the scission reactions. Also, the term, tr( L) vanishes for an incompressible polymer. As the active chain concentrations provide rigidity in the polymer network, any change in the active chain concentration would directly in uence the mechanical properties and stored energy of the polymer. 2.3.2. Active chain ratio and the dissociated chain ratio The chain scission and crosslinking events due to the oxidation reactions in uence the mechanical properties of the polymer network. In addition, the evolution of the active chain concentrations a ect the network elasticity and the stress distribution in a polymer sample. 13 We de ne a non-dimensional parameter as a ratio between the current chain concentration c(t) to the total chain concentration CT0at the virgin state expressed as, crc=c(t) CT0withcrc1 (17) This nondimensional quantity termed as "active chain ratio", attains a value of 1 for a fully crosslinked network. A value of crc= 0 represents a network without any crosslink (dilute chains). In an oxidative aging situation crcjt=0=cr0represents the initial crosslinking density of the network, before any oxidation happened. Further, we also de ne another ratio termed as, dissociated chain ratio crd, which measures the concentrations of chains per mole that has been dissociated from the network ( cd(t)), normalized with respect to CT0as, crd=cd(t) CT0withcrd1 (18) wherecrd= 1 represents a network completely dissociated and crd= 0 represents a fully crosslinked network. 2.3.3. Macroscopic property degradation due to oxidative embrittlement The continuous occurrence of chain detachment/attachment events during oxidation alters the mechanical properties of the polymer network such as the shear modulus ( Gshear), the limiting chain stretchability ( L) and the Grith like fracture energy ( Gc) etc. We assume that the mechanical sti ness or the shear modulus of the network is directly proportional to the active crosslinking density. Hence the shear modulus is expressed as a function of crc, asGshear(crc). In addition, due to chain scissions the molecular weight of the chains reduces and the ductility of the smaller chains drops signi cantly as a function of oxidative reactions [39, 40]. To quantify the drop in ductility we use the limiting stretchability of the network and express that as a function of scissioned chain concentration ( crd). Similarly the fracture energy of the network also reduces as the average molecular weight of the network continues to drop [39, 40]. Hence, we assume the functional relations for LandGc, as,  = i(crd) , where [L;Gc] (19) 14 The actual forms for the function  and  i;i= 1;2; will be discussed later in the speci c constitutive forms. It is important to note that during the set of chain reactions prevailing oxidation, the degraded network essentially formed with smaller chains but more number of crosslinks [40]. Hence, there are studies where molecular weight of the new network has been used to characterize the increasing brittleness [39, 40]. In the present work, we have used the concentrations of the associated and dissociated chains which can directly be quanti ed with the reaction rates to correlate with the mechanical property changes. 2.4. The chain distribution tensor and the free energy of the oxidation reacted network In order to connect the chain distribution function (r;t) with the macroscopic description of the state of the polymer chain network at any stage of the deformation process, the chain distribution tensor ( ) is introduced as the second moment of the distribution, following [87, 88]. As mentioned earlier, assuming a Gaussian distribution for P0(r) at the stress-free con guration, 0can be written as, 0=hP0~r ~ri (20) Evaluating the above expression using Eq. 7, results into identity tensor, I. Other than the stress-free con guration, is written as, =hP~r ~ri. Further assuming the distribution P(r) still remains as Gaussian as the oxidation progresses, with di erent distribution param- eters, can be written as a symmetric (3 3) tensor. As it is well known, the symmetric tensor can be better represented by three of its eigen values i; i= 1;2;3; associated with its principal planes. The eigen values are related to the square of the standard deviation of the current end-to-end distance  ri2as,i= 3ri2=Nb2; whereNb2=3 is the mean-square end-to-end distance of the chains at stress-free con guration. More speci cally, can be directly connected to the average chain stretch as, =r r0=r tr 3(21) where, r0is the root-mean-square (average) distance of the chain end-to-end vector at the stress-free con guration and  ris the current root-mean-square value of the chain end-to-end 15 distance. By the very de nition, the average chain stretch can be related to macroscopic right Cauchy-green strain tensor as follows, =1p 3p tr(C) (22) In elastomeric materials, the stored energy arises mainly due to entropic contribution. Hence the free energy of the material can be expressed as a function of based on the statistical theory of entropic elasticity as, R= R Gshear; L (23) The speci c form of the function Rcan be chosen as reported in [7, 43, 53] to incorporate the limiting chain extensibility of the network. To obtain the evolution equation for , one needs to integrate the evolution equation for (Eq. 15) over the entire chain space. Finally it turns out as, _=ka(4+5+6)CTc(t) c(t) 0(hkdP~r ~ri)(1a+1b+3)_c(t) c(t)+D+D(24) where D= 1=2(L+LT) is the symmetric part of the velocity gradient. The derivation of Eqn. 24 has been presented in earlier literature [87, 88], and hence avoided here. The term (hkdP~r ~ri) can be simpli ed by assuming that the dissociation rate kdis independent of rand hence constant. Hence, we get, ( hkdP~r ~ri) =kd. Replacing it in Eq. 24, we can re-write, _=ka(4+5+6)CTc(t) c(t) 0kd(1a+1b+3)_c(t) c(t)+D+D (25) The two speci c cases for evolution have been provided in the appendix. It is important to note that the chain distribution may not remain as Gaussian for the entire deformation history (especially in the regime, when the chain end-to-end distance attends full contour length as r!Nb). For non-Gaussian distribution, simpli cation of via a symmetric ten- sor is not possible. In the present study we assume that the chain distribution remains close 16 to Gaussian for the entire oxidation and deformation history and the kinetics of the chain association/dissociation events are independent of the mechanical force. 2.5. A phase- eld fracture model to predict macroscale damage and fracture in polymers due to oxidation In the previous section, we have proposed a method for predicting the microscale evolu- tion of the network morphology depending on the oxidative reaction kinetics. In this section we propose a phase- eld fracture theory to model macroscopic damage initiation and prop- agation in terms of cracks during a stress-coupled oxidative aging process. The nal goal is to connect the reaction modi ed polymer network response with the macroscopic frac- ture originating due to stress inhomogenity, during an oxidative aging process. The use of a phase- eld fracture model would help in predicting the fracture initiation and modeling crack propagation in complex conditions such as mixed mode, interacting cracks, complex distribution of initiation sites etc., under a coupled physics driven scenario such as oxidative aging. The model builds upon a coupled mechanical and oxygen di usion response, driven by the displacement boundary conditions and the chemical potential gradients. It further con- siders a reaction dependent amorphous polymer network with constantly evolving network morphology and mechanical properties . 2.5.1. Phase- eld order parameter or macroscale damage variable In the phase- eld fracture theory, we de ne a phase- eld order parameter or damage variabledsuch that, 0d(X;t)1 (26) whered= 0 represents the material being intact and d= 1 corresponds to a fully broken material. Thus, the phase- eld order parameter or damage variable interpolates between the broken state of the material and its unbroken state [2, 61, 63]. We assume that d monotonically increases, such that _d0, implying that the damage is irreversible. In 17 order to de ne a di usive crack topology to approximate the non-smooth crack eld, the phase- eld theory introduces an exponential function for d(X), asd(X;t) =eX lc[61{63]. The length scale parameter lcis introduced to govern the regularization or the width of the di usive zone and gives the sharp crack topology for lc!0 (Fig. 2). The theory further considers dependence on dand the gradientrdto incorporate the spatial inhomogeneity due to localized damage formation. Figure 2: Sharp and di usive crack topology: a sharp crack  embedded into the solid B(left) and the regularized crack surface  lcwith a functional of crack phase eld dand the crack regularization length lc. 2.6. A Chemo-mechanically coupled thermodyamically consistent theory for polymer oxidation:- the balance laws 2.6.1. Mass balance for the di using and reacting species LetPdenotes any arbitrary domain of the body BRin the reference con guration sur- rounded by the boundary @P. Letc R(X;t) denotes the concentration of species in moles per unit reference volume which takes part into the oxidation reactions. The changes in the concentration of c Rfor any species are caused by di usion of that species across the bound- ary@P, which is characterized by a ux j R(X;t), de ned by the number of di using species measured per unit area per unit time. Similarly, r R(X;t) denotes the number of moles for 18 species that is either consumed or produced per unit time due to the chemical reaction. The general mass balance law for any participating species therefore takes the form, _Z Pc RdvR=Z @Pj R:nRdaR+Z Pr RdvR (27) for every part P; where nRrepresents a unit vector in the reference con guration. Bringing the time derivative in Eq. 27 inside the integral and using the divergence theorem on the integral over @P, we nd that Z P(_c R+ Div j Rr R)dvR= 0 (28) SincePis arbitrary, this leads to the local form of mass balance as, _c R=Divj R+r R (29) Using Eq. 3 into 29, we can write, _c R=Divj R+X n[Rn _n] (30) We have assumed oxygen is the only di using species. Also from Eq. 1, O 2is consumed in propagation reaction (II) and produced in termination reaction (VI). Hence, the mass balance equation for O 2has the following form: _cO2 R=DivjO2 R+R22+R66 (31) For the other constituents, mass balance would involve only the reaction dependant changes as there is no di usion, as given by [52] and re-written here. The balance of POOH concentration is thus given by, _cPOOH R =R1_1+R3_3 (32) The rate of change of PH-concentration can be tracked as, _cPH R=R3_3 (33) The termination products due to oxidation are the carbonyls (POOP). We can write the mass balance for POOP as, _cPOOP R =R5_5+R6_6 (34) 19 2.6.2. Balance of forces and moments: Principle of virtual power Using the principle of virtual power, we can derive the macroscopic and microscopic force balance conditions. Following similar notation as mentioned in [49], we conisder the rate-like kinematical descriptors as _ ,_F,_n,_dandr_d. For the macroscopic force balance, we consider that, on the arbitrary part Pof the reference bodyBR, there exists a traction force tR(nR) that expends power over the velocity _ . There also exists a stress tensor, TRthat expends power over the rate _Fof the deformation gradient. Similarly, we consider a scalar microscopic stress $, expending power over the rate of the damage variable, _d. Additionally, a vector microscopic stress is introduced that expends power over the gradient r_d. Finally, there also exists a scalar microscopic traction (nR) that expends power over _don the boundary of the part. Then we can write the external and internal power of the system as, Wext=Z @PtR(nR):_ daR+Z PbR:_ dvR+Z @P(nR)_ddaR (35) Wint=Z P TR:_F+$_d+:r_d dvR (36) where bRis the non-inertial body force per unit volume of the reference body. Introducing the virtual rate elds as ( ~ ;~F;~d;r~d), we can further write the external and internal expenditures of virtual power as, Wext;V=Z @PtR(nR):~ daR+Z PbR:~ dvR+Z @P(nR)~ddaR (37) Wint;V=Z P TR:~F+$~d+:r~d dvR (38) Now the principle of virtual power requires that Wext;V=Wint;V, andWint;V= 0 for rigid virtual velocityV, for any part P. The consequence of these requirements provides us the macroforce balance consistent with the Piola stress as, DivTR+bR=0and TRFT=FTT R (39) Eq. 39 represent the local macroforce and moment balance for the reference body. As is 20 standard, Piola stress is related to the symmetric Cauchy stress in the deformed body by, T=J1TRFT(40) And the microstresses and$are consistent with the microforce balance equation together with the microtraction condition as: Div$= 0 and (nR) =:nR (41) 2.6.3. Energy imbalance Let us consider an isothermal polymer oxidation process in any arbitrary part Pof a reference body BR. Neglecting the inertia e ect, the internal energy in Pcan be written as less than or equal to the power expended on the part due to chemical reactions, damages etc. plus the energy increase due to change in concentration in the species due to the di usion. Thus, denoting Ras the internal energy density per unit volume of the reference body, we can write the energy imbalance as, _Z P RdvRZ @PtRnR_ daR+Z PbR_ dvR+Z @P(:nR)_ddaRX Z @P j RnRdaR(42) where is the chemical potential for any species and j Ris the energy carried by species into partPby the ux j R. Applying divergence theorem to the terms on the integrals over@Pand bringing the time derivative inside the integral, Eq. 42 can be re-written as, Z P_ RdvRZ P TR:_F+ Div TR:_  dvR+Z PbR_ dvR+Z P Div_d+:r_d dvR X Z P  Divj R+j Rr  dvR(43) Reorganizing Eq. 43, we have, Z P_ RdvRZ P TR:_F+:r_dX   Divj R+j Rr  dvR+ Z P(DivTR+bR)_ dvR+Z P(Div)_ddvR(44) 21 Using Eq. 30, 41, and 39 in Eq. 44 and considering Pas arbitrary, we could write the local form of energy balance as, _ RTR:_F$_d:r_d+X j Rr +X  X n Rn _n X  _c R0 (45) Since the stress measure TRis not symmetric, we introduce the 2nd PK stress Sand its conjugate right Cauchy-Green strain Cto express the stress power as, TR:_FS:1 2_C (46) and Eq. 45 can be nally written as, _ RS:1 2_C$_d:r_d+X j Rr +X  X n Rn _n X  _c R0 (47) 2.6.4. Basic constitutive equations In this Section, we consider a theory which allows for an energetic and dissipative e ects associated with the temporal changes in d, and also an energetic e ect due to the gradient rd. We consider the gradient rdas a measure of the inhomogeneity due to localized damage formation. We also consider that scalar microstress, $, can be broken down into an energetic and dissipative part respectively, based on the similar approach proposed in [6, 67] as, $=$en+$dis (48) Guided by the free energy imbalance Eq. 45, a functional form is assumed for the free en- ergy R=^ R(), where = (C;c R;n;d;rd); the second Piola stress S, vector microstress and the chemical potential of species  are then determined by the constitutive equations of the form, S=^S(), = ^ (),=^() (49) Further, the species ux j Rcan be assumed to have the following constitutive form: j R=^j R() (50) Then, by invoking the thermodynamic restrictions it is possible to generate the speci c forms for the constitutive equations. 22 Thermodynamic restrictions:-. Following Eq. 49, we get _ R=@^ R @C:_C+@^ R @d_d+@^ R @rd:r_d+X @^ R @cR _c R+X n@^ R @n_n (51) In view of Eq. 51, the free energy imbalance Eq. 45 is equivalent to the requirement that the following inequality must be satis ed for all the constitutive processes: @^ R @C1 2^S! :_C+ @^ R @d$! _d+ @^ R @rd! :r_d+X @^ R @cR ^ ! _c R " X  X nRn X n@^ R @n! _n# +X ^j Rr 0 (52) To hold the above equation true for any arbitrary _C, _c Randr_d, their coecients must vanish. Therefore, we can get the 2nd P-K stress, S, the microstress, and the chemical potential, as, S= 2@ R() @C; =@ R() @c R(53) We also consider the vector microstress as energetic, and thus, =@ R() @rd(54) We de ne the energetic microstress $enas, $en=@ R() @d(55) so that the disssipative microstress $discan be written as, $dis=$$en (56) The species ux satis es the species-transport inequality as, j R()r 0 (57) In this work, we assume that the uid ux obeys Fick's law, i.e. the ux depends linearly on the gradient of the chemical potential and thus can be expressed by, j R=M r (58) 23 withM is the mobility tensor. Substituting Eq. 58 into 57, we get the species-transport inequality as, r M()r 0 (59) Hence, we are nally left with the following dissipation inequality:- $dis_d+" X  X nRn X n@^ R @n! _n# 0 (60) We assume that each individual term of Eq. 60 satis es the dissipation inequalities as: 8 >>< >>:$dis_d0 P  P nRn P n@^ R @n _n0 ,(61) From Eq. 61, we consider that $disis the driving force for the phase- eld damage d. In addition, from the chemistry literature, we know that there exists a quantity called anity ,An, which is the driving force equivalent, as, An=@ R @n(62) Hence, it is possible to de ne a force of dissipative nature Fnwhich is a conjugate to _nfor each reaction n. Thus, Fn=AnX  Rn ; (63) This re ects the fact that every chemical reaction relates to a dissipative mechanism. The extent of reaction for any reaction nis assumed to evolve according to a state relation _n=^_n(Fn;#;n;c R)0 (64) withF_n>0, whenever _n>0, and#is the reference temperature. 2.7. Specialization of the constitutive equations For deriving the speci c constitutive forms, we consider the materials as initially and continually isotropic. Under such assumption, it can be shown that the response of the poly- mer is invariant under any arbitrary rotation in the reference con guration or intermediate con guration and the constitutive responses do not change [3{5, 22, 23, 56]. 24 To simplify the numerical implementation, we consider a Gent hyperelastic free energy function for the mechanical part [43]. It is important to note that the Gent free energy ac- counts for the limiting chain extensibility in a phenomenological way. However, its statistical basis can be proved based on an approximation of molecular based average stretch for the full network, as reported in [51]. For the chemical part of the free energy, a simple quadratic function of the extent of reaction (n) is considered, (based on the recent literature [52, 56, 75, 84]). Following the formulation in [52], we have assumed only oxygen di usion takes place within the polymer during oxida- tion. Hence, the free energy for the di usion is considered as the energy of mixing between the di used oxygen and the polymer. Polymer being the long-chain macromolecules and the size of the di used oxygen molecules being signi cantly smaller, Flory-Huggin's theory of mixing seems to be better suited to model the di usive free energy [55, 77, 78]. Hence, we use a similar expression for di usive free energy as used in earlier works of uid di usion in the polymers ([22{24, 41]). It is important to note that as the phase- eld damage progresses, the stored or mechanical energy of the material continues to degrade. In order to model that, we incorporate a widely used degradation function g(d) as reported in the phase- eld literature [6, 59, 61, 63, 67]. The function g(d) is a monotonically decreasing function in dand it satis es:- g(0) = 1, g(1) = 0,g0(d)0 andg0(1) = 0. We use the most common form for degradation function used in the literature [6, 61, 63, 67], as, g(d) = (1d)2+ (65) whereis a small positive-valued constant which is introduced to prevent ill-conditioning of the model when d= 1. Free energies:- .Motivated by the free energy representation in [6], we consider that the free energy can be written in a separable form constituted of mechanical energy, chemical energy, energy of mixing and the fracture energy. ^ R(C;c R;n;d;rd) =g(d) mech R(C) + chem R(n) + di R(cO2 R) + frac R(d;rd) (66) 25 Mechanical free energy :- .In what follows, we write the mechanical free energy following Gent hyperelastic model [43] in terms of the average stretch as, mech R(;J) =1 2GshearLln 13(21) L +Kbulk 2(J1)2(67) where, we have used the 1st strain invariant as, I1= 32=tr(C),J=det(F), the shear moduludGshear varies as a function of crc,Kbulkis the bulk modulus of the polymer and L is the limiting stretchability, dependent on crd. Chemical energy :- . chem Ris the reactive part of the free energy in uencing the oxidation reaction and is given by, chem R=X nHn 2(1n)2(68) where the parameter Hnis the chemistry modulus of reaction n. This term in the free energy favors the local state ( n= 1) [84]. Energy driving the di usion of oxygen:- . diff Ris the part of free energy involved in the di usion of oxygen (O 2) into the polymer. As mentioned we have Flory-Huggins free energy to de ne the energy for di usion [22{24, 41] as given by, di R=O2 0cO2 R+R#cO2 R ln cO2 R 1 + cO2 R + 1 1 + cO2 R (69) whereO2 0is the reference chemical potential for oxygen, Ris the universal gas constant, # is the temperature, is the dimensionless Flory-Huggins interaction parameter, and is the volume of a mole of oxygen. Fracture energy :- .The free energy for the fracture is also considered as a quadratic form, frac R=Gc(crd)lc 2jrdj2(70) whereGc(crd) is the Grith-like fracture energy that evolves as a function of chain dissocia- tion ratio,crd. 26 Thus the complete free energy expression can now be written as, ^ R=g(d) 1 2GshearLln 13(21) L +Kbulk 2(J1)2 + X nHn 2(1n)2+O2 0cO2 R+R#cO2 R ln cO2 R 1 + cO2 R + 1 1 + cO2 R +Gclc 2jrdj2(71) Evolution equation for the phase- eld variable d:-.Based on Eq. 55, we can write the energetic microstress as, $en=@ R @d=2(1d) mech R (72) where, mech R is the mechanical free energy of the reaction modi ed polymer network, but without any macroscopic damage. To determine the expression for $dis, let us consider the recent phase- eld literature [61, 63], where the energy required to create a new fracture surface was given by, s=Z PGc1 2lcd2+1 2lcjrdj2 dV (73) In order to be consistent with the energy dissipated to create new surface (as the rst term in Eq. 73), we de ne the dissipative microstress as, $dis=Gc(crd) lc(74) It is important to recall that the modi cation of the Grith energy would occur due to the chain scissions occurring in the polymer network during oxidation, as mentioned in Sec.2.3.2 . Therefore, Eq. 48 can be re-written as, $=$en+$dis=2(1d) mech R+Gc lc(75) Further following Eq. 53, we can derive the vector microstress as, =Gclcrd (76) Substituting Eq. 75 and 76 into the microstress balance Eq. 41 yields the evolution equation fordas, 2(1d) mech RGc lc+Gclcd= 0 (77) 27 where drepresents the Laplacian of d. It is important to note that the evolution of din Eq. 77 does not ensure the irreversibility constraint of the damage ( _d0). In order to ensure the irreversibility condition, we can rewrite Eq. 77 while neglecting the gradient e ect ( rd) on the damage variable as, 2(1d) mech RGc 2lc Gc lcd= 0 (78) It can be shown, that a monotonically increasing evolution of drequires the the term mech R Gc 2lc to be strictly positive for the entire deformation history, [6, 61, 63, 67, 79]. Based on that a monotonically increasing history eld His introduced as, H(t) = max s2[0;t] mech R(s)Gc(s) 2lc (79) So that the evolution equation has the nal form as, 2(1d)HGc lcd+Gclcd= 0 (80) Stress, chemical potential and anity:- .From the free energy expression Eq. 71 and using the guidelines for thermodynamic restriction, it is possible to get the speci c consti- tutive forms for the Piola stress TR, chemical potential of oxygen, O2and anity of each reaction,An. The Piola stress can be expressed in the following form: TR=J2=3g(d)" Gshear 1tr(C)3 L1 [F1 3(tr(C)FT] +KbulkJ(J1)FT# (81) wheretr(C) can be related to the average macroscopic stretch as tr(C) = 3 2, which can further be related to as explained before; the shear modulus Gshear and the limiting stretch- abilityLare functions of crcandcrd, respectively. Chemical potential driving the di usion of O 2,O2, can be obtained by taking the partial derivative of the di usive energy with respect to oxygen concentration as, O2=O2 0+R# ln cO2 R 1 + cO2 R +1 1 + cO2 R+ 1 1 + cO2 R2! (82) Finally, the anity of any reaction ncan be calculated from Eq. 62 and 71 as, An=Hn(1n) (83) . 28 Evolution of extent of reaction n:.Based on the kinetics of chemical reaction for polymer oxidation and our discussion in section 2.1, we choose a thermally activated relation for the evolution of extent of reaction _nas, _n=knexpQn act R# Fn (84) where,knis the pre-exponential rate constant for reaction nand has an unit of1 MPas, Qn actis the activation energy for reaction n. As explained earlier, reaction would occur only when Fn>0. Further, the reaction would be considered as completed when the extent of reaction reaches unity. Following this, we can rewrite Eq. 84 as, _n=8 >>< >>:knexp Qn act R# Fn;whenFn>0 andn<1 0;otherwise(85) Oxygen ux and di usivity:- .From Eq. 58, oxygen ux is given by, jO2 R=mO2rO2(86) withmO2being the mobility of oxygen. Ignoring the e ect of temperature gradient we can write from Eq. 82, rO2=R# cO2 R1 1 + cO2 R cO2 R (1 + cO2 R)22 cO2 R (1 + cO2 R)3 rcO2 R (87) Further, di usivity of oxygen can be de ned as, D=mO2R# cO2 R(88) Using Eq. 88 in 86 the constitutive relation for oxygen ux can be obtained as, jO2 R=D1 1 + cO2 R cO2 R (1 + cO2 R)22 cO2 R (1 + cO2 R)3 rcO2 R (89) 29 Figure 3: Linear trend of fracture energy and limiting chain stretchability paramater as a function of disso- ciated chain ratio Reaction dependent shear modulus, fracture energy and limiting chain extensi- bility parameters:- .According to the description in section 2.3.2, we assume that, shear modulus is directly proportional to the active chain ratio and can be written as, Gshear(crc) =crcG0=crcc(t0)kB# (90) whereG0is the virgin shear modulus of the material which is a function of active chain concentration at virgin state c(t0), andkBis the Boltzmann constant. Fracture energy ( Gc) and the limiting chain stretchability parameter ( L) have been consid- ered as degrading linearly with the dissociated chain ratio (as shown schematically in Fig 3), characterizing the increase in brittleness due to oxidative scissions in the network. Thus, we de ne the reaction dependent fracture energy via a linear relation as, Gc(crd) = 1crd+C1 (91) with 1=Gc;oxGc;un crd;oxcrd;unandC1= [Gc;ox+Gc;un 1(crd;oxcrd;un)]=2; whereGc;oxandGc;un are the fracture energy at completely oxidized and virgin state respectively, and crd;ox and crd;un are the corresponding dissociated chain ratios. Similarly, we de ne the limiting chain extensibility parameter as, L(crd) = 2crd+C2 (92) 30 where, 2=L;oxL;un crd;oxcrd;unandC2= [L;ox+L;un 2(crd;oxcrd;un)]=2; withL;oxandL;un are the stretchability at completely oxidized and virgin state respectively,. 2.8. Governing di erential equations and the boundary conditions There are three governing partial di erential equations required to be solved in this case: 8 >>>>>< >>>>>:DivTR+bR=0 _cO2 R=DivjO2 R+R2_2+R6_6 2(1d)HGc lcd+Gclcd= 0(93) We require the initial and boundary conditions to complete the solutions of these di er- ential equations. Let, S1andS2are complementary subsurfaces of the boundary @BRof a reference body BRsuch thatS1[S2=@BRandS1\S2=?. Similarly, let ScRandSjRare complementary subsurfaces of the boundary @BR=ScR[SjRwithScR\SjR=?. Lastly, we consider that, there exist two complementary surfaces SdandSof the boundary such that @BR=Sd[SwithSd\S=?. For a time interval t2[0;T], three boundary conditions can be described such that displacement is known on S1and traction on S2. Thus, we can write, =  onS18t2[0;T] TRnR=tRonS28t2[0;T] (94) Similarly, another pair of boundary condition can be considered for a time interval t2 [0;T] such that oxygen concentration is known on ScRand oxygen ux on Sjand thus, cO2 R= cRonScR8t2[0;T] (jR):nR=jRonSjR8t2[0;T] (95) Finally, the boundary conditions for the evolution of d are, d= 0 onSd8t2[0;T] (rd):nR= 0 onS8t2[0;T] (96) 31 The initial conditions are, (X;0) = 0;andcR(X;0) =c0 Randd(X;0) = 0 inBR: (97) Thus the coupled set of equations, Eq. 93, together with 94, 95 and 97 pose an initial boundary value problem to be solved for the displacement (X;t), oxygen concentration cO2 R(X;t) and phase eld parameter dsimultaneously. In this work, the system of equations are solved numerically by writing a user element subroutine (UEL) in ABAQUS/Standard (2017) [1]. 3. Material Parameter Estimation To numerically implement the present theory, a number of material parameters are re- quired for modeling the coupled behavior of multiple chemical reactions-di usion, large de- formation and phase eld fracture in polymer. Estimation of the kinetic parameters for the di usion and chemical reactions has been explained in detail in our previous work [52]. Hence, we ignore that part here for the sake of brevity. The kinetic parameters used in the present work are listed in B.2. As described in Section 2.3.2, chemical reactions in the oxidation process cause chain- scissions (associated with the rate coecient kd) and oxidative crosslinking (associated with the rate coecient ka); which in turn, leads to the evolution of the active chain concentration and chain distribution. In order to determine these rate-like reaction dependent coecients we have used experimental data as reported in [39]. A simple curve tting procedure is adopted assuming an exponential relation between the concentration data for chain-scissions and crosslinks as a function of oxidation time, as shown in Fig. 4. As evident from the data in [39], the chain-scission phenomenon essentially dominates over the oxidative crosslinking in this typical thermoplastic polymer. On the other hand, crosslinking event could be more prominent in the case of some typical thermosets during oxidation as reported in [35{37]. Thermo-oxidative aging signi cantly a ects the critical strain energy release rate or the fracture toughness ( Gc) of the material, reducing it by two or three orders of magnitude as reported in [38, 39]. As identi ed, the chain-scission mechanism destroys the polymer 32 Figure 4: Estimation for the kaandkd, the rate-like coecients for polyethylene; the "markers" correspond to the experimental data from [39]. network structure to a point where plastic deformation process do not initiate and hence the embrittlement occurs [39]. In this work we choose a Grith-like fracture energy release rate as 1:44 kN=mm for virgin polymer and 10 N =mm for a completely oxidized polymer, consistent with the values as reported in the literature [38], for a typical polymer. The other parameter to describe the mechanical free energy of the polymer network is the limiting chain extensibility ( L). Oxidative embrittlement process notably reduces the ultimate failure strain of the material, as the stretchability drops from a value of 800% to 10% as observed in the case of oxidative aging of polypropylene, [38]. For the present simulation, we useL= 3:0 for the unoxidized polymer and L= 1:9 for the completely oxidized polymer. TheGcandLvalues used in the simulations are tabulated in the appendix (Table B.1). 33 4. Representative numerical simulations 4.1. Case study-1: Oxidative aging of a dogbone specimen followed by tensile loading In this example, we numerically simulate the oxidative aging of a representative polymer dogbone specimen followed by a tensile testing to imitate a typical oxidation experiment performed in a laboratory setting. In the simulation, we consider the initial shear modulus asG0= 0:58GPa and a nal oxidized modulus as Gshear;ox = 0:7 GPa. The other important mechanical and kinetic parameters are as listed in Table B.1 and B.2 respectively. The geometry of the dogbone sample and the corresponding FE mesh have been shown in Fig. 5(a). Due to the symmetry, we analyze only one-quarter of the sample using 2D plane strain elements and apply symmetric boundary conditions along the mid-planes in the x and y- directions. The length scale parameter to describe the di use damage zone has been chosen as,lc= 500m. In a traditional lab environment, a sample is rst oxidized within an oven at a temperature close or higher than the glass transition temperature of the material for several hours. Subsequently, it is taken out to perform mechanical testings to determine the e ect of oxidation on the material's response. In what follows, we perform a two-step simulation where at rst we let the sample to be oxidized at 80oCfor several hours in an accelerated aging environment. In the second step, we pull the sample at a displacement rate of 2:5105mm=s until the nal failure happens. The corresponding chemical and mechanical boundary conditions are shown in Fig. 5(b). In the rst step, we consider that, the sample surface is surrounded by a constant ow of O 2with partial pressure of 2 atm, which translates to a constant O 2concentration of 2 :8102mol=L. In order to avoid convergence issue, we have applied the chemical boundary condition in a ramp-like manner over a small time period of 2 hours. In the following step, we have applied the displacement boundary condition on the top surface of the specimen in a ramp-like manner as well (Fig. 5(b)). Our aim is to study the e ect of oxidative aging in terms of the mechanical property degradation upon oxidizing the sample for 30, 50 and 100 hours, respectively. We consider two speci c cases to understand the kinetics of the chains- i) when chain scission rate is dominant over the reattachment rate ( ka< kd) and ii) when the chain reattachment is higher ( ka> kd) 34 compared to the scission rate. Since oxidation is a time-dependent process, longer exposure to oxidative environment inevitably leads to the higher values of the extent of reactions, which, in turn, yields harsher property degradation, as demonstrated in the simulations results, as discussed in the following paragraph. (a) (b) Figure 5: (a) Geometry and corresponding 2D plain strain nite element mesh of a dogbone specimen; (b) boundary condition for the oxygen concentration applied at the edges of the dogbone sample (top) and the displacement boundary condition applied at the top surface of specimen in Step-2 (bottom). Fig. 6(a) and (b) show the contour plots of extent of reaction (III) ( 3) and extent of reaction (VI) ( 6) at the end of 30 hours (left), 50 hours (middle) and 100 hours (right), respectively. The corresponding changes in Gshear andLare shown in Fig. 7 and Fig. 8, respectively. The extent of reactions monotonically increase in a continuous oxidation process. As time progresses, the evolution of 3leads to resin depletion by breaking the PH-chains and formation of some new hydro-peroxides (POOH). On the other hand, the termination reaction 35 (a) (b) Figure 6: Contour plots for (a) extent of reaction 3after 30 hours (left), 50 hours (middle) and 100 hours (right); (b) extent of reaction 6after 30 hours (left), 50 hours (middle) and 100 hours (right). (VI),6forms new PO-OP bonds and these new bonds accumulate within the material as the reactions progress. As explained in the theory, the chain distribution changes as a function of oxidation reactions and alters the active chain concentrations c(t). Following Eq. 12, let us consider the chain detachment happens as a function of 3,fd=kd3. Similarly, the chain reattachment happens as a function of 6asfa=ka6. Forkd> kacondition, the number of active chains decreases leading to a decrease in the shear modulus ( Gshear) as seen in Fig. 7(a), (b) and (d). On the other hand, when ka> kd, we observe an increase in Gshear, as demonstrated in Fig. 7(c). The contour plots of the limiting stretchability parameter, L 36 (a) (b) (c) (d) Figure 7: Contour plots of shear modulus after (a) 30 hours, (b) after 50 hours (with ka<kd), (c) 50 hours (ka>kd) and (d) 100 hours (in MPa). (a) (b) (c) Figure 8: Contour plots of Lafter (a) 30 hours, (b) 50 hours and (c) 100 hours (with ka<kd). at various stages of oxidation are also shown in Fig. 8. It is important to note that Lis essentially a function of the dissociated chain ratio, crd. As the oxidation progresses, the number of dissociated chains goes up, leading to an increase in crdfor theka<kdcondition; which in turn causes a reduction in the stretchability, as shown in Fig. 8(a), (b) and (c). The changes in the material properties further lead to a change in the constitutive response of the polymer. Fig. 9 shows the force versus displacement plots for the dogbone samples 37 Figure 9: Force vs displacement plots for the dogbone samples subjected to tensile loading upon oxidation for various hours. Figure 10: Contour plots of the phase eld variable d at signi cant points in the force-displacement plot as marked a, b, c and d in Fig. 9. oxidized for several hours. In general, a signi cant drop in the limiting stretchability and tensile strength has been observed due to oxidation. Approximately 18 :33% drop in the failure 38 strain has been noticed after 100 hours of oxidation. The maximum tensile strength of the material also decreases signi cantly. We also observe a slight change in the sti ness of the material during the initial period of oxidation (between 30 to 50 hours), which signi cantly decreases upon further oxidation (at 100 hours). Fig. 10 shows the damage contours of the dogbone sample corresponding to the signi cant points in the load-deformation curve as denoted in Fig. 9 for 50 hours of oxidation (corresponding to ka<kdcondition). At point a, damage initiates with a very small value at 0:3 mm of displacement. The force continues to increase beyond this point until it reaches the peak at point b corresponding to a displacement of0:4 mm. The respective contour of the phase eld variable is shown in Fig. 10 b. The continued loading of the sample ultimately causes rupture of the material, as shown by c and d in Fig. 9 with corresponding damage contours in Fig. 10c and d, respectively. Oxidation in the presence of mechanical loading (stress-coupled):- .In order to identify the e ect of external loading on the polymer sample simultaneously undergoing ox- idation, we consider the same dogbone geometry as shown in Fig. 5(a). The displacement (mechanical) boundary condition as applied to simulate the stress-coupled situation is demon- strated in Fig. 11(a). We have applied a total displacement of 1 mm on the top surface of the sample in a ramp-like manner over the period of 2 hours and then hold it for a total of 50 hours. The chemical boundary condition remains the same as Fig. 5(b) for the 50hrs period. Hence the sample is allowed to oxidize for 50hrs under a constant displacement of 1mm. Fig. 11(b) shows the load vs. time plot for the simulation. A quantitative comparison between the stress-coupled situation and the two-step analysis (case 1) is not realistic as the boundary conditions are completely di erent. From Fig. 11(b), we see that the material behavior is dominated by the mechanical deformation for the rst 2 hours. Subsequently, the response is predominantly driven by the oxidation process, as the displacement is held constant after 2 hrs. The material eventually fails as an e ect of the degraded mechanical properties caused due to the network embrittlement during oxidation. Since the displacement is held constant at a value of 1 mm, we see the load starts to drop at the corresponding value until complete failure happens in Fig. 11(c). 39 (a) (b) (c) Figure 11: (a) Displacement boundary condition for the dogbone sample for a stress-coupled oxidation situa- tion; (b) the force vs. time plot for a constant deformation+oxidation (stress coupled); (c) the corresponding force vs. displacement plot. 4.2. Case study-2: Oxidation followed by Mode-I loading on a single-edge-notched plate In this example, we attempt to predict the response of a polymer specimen containing a pre-existing crack undergoing oxidation. In reality, the presence of cracks promotes oxygen di usion providing favorable paths. Hence pre-existing aws make the material more prone to oxidation near and around the crack-zones. To simulate one such scenario, we choose to study the e ect of oxidation on a single-edge-notched polymer plate, as shown in Fig. 12a. As earlier, the kinetic and material parameters are chosen from Table B.1 and B.2, respectively. 40 The chemical boundary condition is considered as the same as in case study 1. In this simulation, we assume that, oxygen di uses from the left side of the sample, and all other edges are impermeable. We let the sample to oxidize at 80oCfor 100 hours. Following which, the oxidized sample is pulled at a rate of 1 103mm/s along the top-edge, till the nal fracture occurs. (a) (b) (c) Figure 12: (a) Geometry of the single-edge-notched tension sample showing the displacement and di usion boundary conditions; Contour plots of the (b) extent of reaction (III), 3showing resin depletion; (c) chain dissociation ratio crdafter 100 hours of oxidation. Fig. 12(b) and (c) show the contour plots of extent of reaction, 3and the chain dissoci- ation ratio, crd, respectively after 100 hours of oxidation. As seen in Fig. 12(b), 3, reaches to a value of 0 :946 at the outermost surface along the left boundary as well as around the crack tip, indicating chain-scission would take place predominantly within this zone. As a result, the chain dissociation ratio continues to increase and reaches a value of 50% on the left surface and around the crack tip after 100 hrs as shown in Fig. 12(c). This phe- 41 (a) (b) Figure 13: Contour plots of (a) Gc; (b)Lafter 100 hours of oxidation. nomenon further leads to a reduction in the limiting chain extensibility and Grith-type fracture energy in the oxidized zones, as shown in Fig. 13. The Grith-type fracture energy drops signi cantly in the oxidized zone, from 1 :44 kN/mm to 10 N/mm after 100 hrs of oxi- dation, as seen from Fig. 13(a). The limiting chain extensibility also drops down to 2 :447 on the left surface compared to the unoxidized value of 3 :0 for the virgin polymer as shown in Fig. 13(b). Consequently, these events a ect the load-displacement response of the material as demonstrated in Fig. 14(a). After 100 hrs of oxidation, we see a signi cant drop in the maximum load, compared to the virgin material (as 62%). Similarly, the stretchability of the oxidized material has been dropped compared to the virgin material, as seen from the ultimate displacement values corresponding to nal fracture. Fig. 14(b) shows the phase eld contour plots of the oxidized polymer sample corresponding to the locations marked by a, b and c in the force-displacement curve. The sample reaches to the maximum load 42 (a) (b) Figure 14: (a) Force vs. displacement plot of the single-edge-notch sample shown in Fig. 12(a) ; (b) Contour plots of the phase eld variable d at signi cant points in the force-displacement plot marked as a, b, and c in Fig. 14(a) for the oxidized polymer. just before the point a. Beyond this point, the load continues to drop with a continuous increase in the phase eld parameter value. Finally the sample completely breaks at point c, with the maximum value of corresponding phase eld parameter reaching 0 :95, as shown in Fig. 14(b)-c. 5. Concluding remark In this work, we have developed a reaction-driven evolving polymer network theory cou- pled with phase- eld fracture to model the high temperature oxidative degradation and fracture in amorphous polymers. The proposed theory considers the closed-loop chain re- action kinetics occurring at the polymer chain level and successfully incorporates them into the statistical mechanics based network theory to account for the network alternation and mechanical property degradation. The proposed theory is capable of quantifying and re ect- ing the network level changes (degradation) in the polymer via the active and dissociated chain ratios; both of which evolve as a function of oxidative crosslinking and chain scission. Upon establishing a seamless connection between the micro and macro-scale kinematics the 43 continuum level constitutive response of the material is further derived. In addition, we have used a continuum level phase- eld fracture theory to model the macroscopic damage and fracture originating due to the prolonged oxidation in the material. We further tested the capability of the proposed model by numerically implementing it in a commercial FE-setting and performing several important case studies. It is observed that the proposed model is capable of predicting the embrittlement of the materials caused by the high-temperature oxidative aging, and macroscopic damage propagation and fracture. In the future, we plan to extend the proposed model to capture the e ect of oxidation-induced embrittlement in the anisotrpic deformation setting in polymer matrix composites (PMC) and perform exper- imental validation. 6. Acknowledgment The authors would like to acknowledge the super-computing resources from the SUPE- RIOR computing facility at MTU. Appendix A. Evolution of the chain distribution tensor while no oxidation (or chemical reaction is taking place), pure deformation Eq. 25 can be concisely written as, _=D+D (A.1) with all other reaction-dependent terms being zero and assuming the un-reacted net- work as incompressible. In such case, the solution for turns out to be =FFT, which is equivalent to the left Cauchy-Green strain tensor as, B=FFT. [88]. when the mechanical deformation is held constant ( F=I) and oxidation (only chemical reactions) is allowed to happen Introducing, = (4+5+6)CTc(t) c(t)(A.2) 44 and =1a+1b+3in Eq. 25, we get, _=ka I kd +_c(t) c(t) + (+T) (A.3) Following which, the solution of turns out as, =Fe(kd +cm)tFT+ka kd +cmFFT(A.4) wherecm=_c(t) c(t), is used. Appendix B. Parameters used in the numerical simulation Kbulk 5 GPa ([54])  0.45 ([54]) G0 0.58GPa ([54]) Gc;un 1:44kN=mm Gc;ox 10N=mm L(un) 3 L(ox) 1:9 crd;ox1:0 c(t0)0:58e9mol=m3 CT(t0)0:7e9mol=m3 Table B.1: Material parameters for a representative elastomer undergoing oxidation 45 Table B.2: Kinectic parameters related to di usion and reaction kinetics for a typical polymer undergoing oxidation ([52]) Parameter Unit Value Reference Solubility coecient, S molL1Pa11:4108[73] Di usivity , D m2s18:8105[72] Di usive activation energy, EdkJ=mol 36:5 [73] Reaction rate constants k2 1=MPas 2:84109calibrated k3 1=MPas 2:9109calibrated k6 1=MPas 2:76109calibrated Q2 kJ=mol 30 [73] Q3 kJ=mol 45 [73] Q6 kJ=mol 10 [73] H2 J=mol 128106[18] H3 J=mol 40:52106[18] H6 J=mol 267:26106[18] Gas constant, R J/(mol-K) 8.31446 0.1 [24] 0at STP J=mol 3:88105[44] c(PH) 0 mol=L 20 [72] c(POOH ) 0 mol=L 1104[72] 46