A thermodynamically consistent chemo-mechanically coupled large deformation model for polymer oxidation

In this paper, we present a continuum-level thermodynamically consistent model for high temperature oxidation in polymers, that incorporates the coupling between diffusion, chemical reaction and large deformation behavior of polymers. The specific constitutive forms are derived based on the thermodynamic inequality conditions and the kinetics of the oxidative reactions are considered. Oxidative shrinkage has also been considered in the kinematics as an irreversible effect. Subsequently, the model is implemented in ABAQUS/Standard to analyze numerically the coupled diffusion-reaction behavior of polymers undergoing oxidation. Several numerical simulations are performed to understand the effect of various material parameters on the oxidative response. The model is capable of predicting the heterogeneous oxidation profile within a thick polymer sample. It can also track the growth of oxide layer in the case of a long-term thermo-oxidative aging process. The model can be used to simulate the oxidation process involving complex geometries (as fiber reinforced composites) fairly easily under various ambient conditions.


Introduction
Polymers and polymer matrix composites (PMC's) are attractive materials to aerospace and automotive industries because of their certain properties, such as light-weight, excellent strength to weight ratio, damping, thermal stability, corrosion resistance and so on. With the continued discovery and widespread use of new PMC's, the challenge lies into finding a reliable and efficient way to predict the lifetime of these materials. The performance of polymers can degrade significantly due to their long-term exposure in extreme environments, such as corrosion, thermal oxidation, chemical environment, etc. High temperature oxidation or thermo-oxidative aging in polymer is one such example, where polymer undergoes a set of irreversible chemical reactions driven by oxygen diffusion. This phenomenon results in formation of the oxide layer on the surface of the material which may lead to spontaneous cracking and eventually cause failure of the materials [15,16,21,25,55,59,60]. This could be a major concern in aerospace and automotive industries because it impacts both the safety and cost-effectiveness of these materials for their long-term applications.
Thermo-oxidative aging of polymers is a highly nonlinear, irreversible phenomena governed by a coupled diffusion-reaction process, resulting changes in the material both at the molecular and macro-scale. These changes in macroscopic properties and molecular structure of polymers due to oxidation have been extensively reported in literature [21, 24, 25, 27, 52-55, 59, 60]. According to these works, at the macromolecular level, chain scission and oxidative cross-linking effect the polymer network and eventually alters the material property during oxidation. At the macroscale, as oxygen infuses into the polymer, it reacts with the material to create oxides and volatile compounds; latter leaving the material to form voids.
Due to the void formation, the material progressively loses weight and simultaneously due to oxidative crosslinking, it becomes denser. Ultimately these two combined factors result in irreversible volume shrinkage within the polymer. More detailed explanation of these can be found in [25]. Quantifiable consequences of oxidation in the polymer are irreversible property changes such as increase in modulus and glass transition temperature, and decrease in failure strain, etc. [13, 30-32, 48, 49, 52, 53]. In addition, other effects of thermal oxidation include d r a f t debonding on the fiber-matrix interface in PMC's and matrix micro-cracking, which might lead to premature failure of the composites [37,54,55,60].
Modeling of polymer oxidation presents a great deal of difficulty as it involves various highly coupled physical and chemical phenomena. In the early studies of thermal aging of polymers and PMC's , the property deterioration and lifetime predictions were performed by conducting accelerated aging experiments and extrapolating these results using Arrhenius laws to predict the behavior for later periods [9,11,12,47,51,61]. However, as mentioned in [26], Arrhenius law can be applied only to an elementary process, whereas polymer oxidation in general involves at least six elementary reactions. Therefore, it might not be sufficient to use a generalized Arrhenius law for such complex cases. The deviation from Arrhenius law has been experimentally proved by Celina and coworkers in many of their research [14,[38][39][40][41]63].
Verdu and coworkers have extensively studied the reaction mechanism of polymer oxidation in the past decade [6-8, 20-26, 52, 57, 64], and developed a kinetic model by incorporating coupled oxygen diffusion and consumption kinetics. An extension of this kinetic model was proposed by Tandon and coworkers [59], [60], [54]; where oxidation process in polymer was described by a three zone model. In the three-zone model, the surface oxidative layer (zone I) is separated from the unoxidized polymer (zone III) with an active reaction zone (zone II).
This heterogeneity in oxidation results in non-uniform mechanical properties and stress-strain distribution in an oxidized polymer.
All these models mentioned above use semi-empirical approach to determine the mechanical properties of oxidized polymers. One such popular approach to correlate mechanical stiffness with the degree of oxidation is to measure the local variation of the indentation elastic moduli in the oxidized layer of a bulk polymer sample [37], [49], [48], [52], [53]. To study the coupled chemo-mechanical behavior of polymer matrix exposed to thermo-oxidative environments, a theoretical model considering the thermodynamics of irreversible process has been first proposed by Gigliotti and coworkers [37]. The method has been further extended to incorporate oxidation induced change in visco-elastic properties of epoxy polymers and experimentally validated with different aging time, temperature and oxygen partial pressure d r a f t in a recent work [36] by the same group. However, the model has developed in small deformation setting, and no diffusion of oxygen has been considered. Recently, diffusion limited oxidation (DLO) profile has been modeled by Celina and coworkers [56] using finite element method where the authors have demonstrated the effect of partial pressure of oxygen on the degree of oxidation in a polymer. However, the proposed methodology has not been used to predict the mechanical response of polymers undergoing progressive oxidation.
It is important to note, that the coupling of multiple physics in the constitutive modeling of materials is not new in the mechanics community. For instance, a coupled theory for high temperature oxidation in metallic alloys has been developed by [46]. The coupled stress-diffusion response of polymer gels has been studied by [17-19, 29, 43] in a great detail.
A very recent work by [66] on the reaction-diffusion coupling in photoresponsive gels is also noteworthy. Recent works on coupled photochemical reaction-large deformation theory has also been published by several research groups [28,58,65].
Similar to that trend in this work, a thermodynamically consistent continuum-level theory is developed that incorporates the coupling between diffusion and oxidation reactions together with the large deformation behavior of polymers. The model is implemented in a commercial FE package ABAQUS/Standard by writing a user element subroutine (UEL), to describe the coupled thermo-chemo-mechanical behavior that arises due to high temperature oxidation in polymers. The paper is organized as follows-the coupled theory of oxidation has been explained followed by the numerical implementation of the specific constitutive model. Various numerical simulations have been performed and the results have been analyzed in detail to understand the capability of the model. At the end, the concluding remarks based on the present work has been presented.

Chemistry of oxidation
With some minor variation at the initiation stage, the mechanistic reaction schemes proposed to describe the thermal oxidation of polymers had a general consensus [6-8, 20-26, 57, 64] in the literature. According to this standard mechanistic scheme, during oxidation, d r a f t polymer undergoes six set of closed loop chain reactions involving three major steps identified as initiation, propagation and termination, respectively, as explained below: where PH represents the polymer substrate; P * is the alkyl radicals; POOH is the hydroperoxide, P O * 2 is the peroxy radical and O 2 represents oxygen. The set of reactions as explained above was originally given by Verdu and coworkers [26]. For a given macromolecular structure, which reactions would be dominating during the oxidation process are determined by comparing the activation energy and reaction rates of each individual reaction [20].
Among the two initiation reactions, activation energy associated with the hydroperoxide decomposition (Ia) is noticeably lower than the substrate decomposition or substrate oxidation reaction (Ib), and is called the initiator of the oxidation reaction [26]. Therefore (Ia) tends to dominate the initiation rate of the oxidation process. On the other hand, at high temperature and low POOH concentration, initiation is predominantly due to the substrate consumption (Ib) and the rate of initiation is almost constant or slowly decreases. This causes significant structural change in the polymer as a result of substrate consumption [26].
In case of propagation reaction, the rate constant k 2 for reaction (II) is much faster than the rate constant k 3 for reaction (III). Therefore, in case of excess oxygen, reaction(III) becomes dominant rate-controlling process, as the propagation will depend on the availability of polymer, [26]. Further, due to the presence of excess oxygen, alkyl radicals (P * ) quickly d r a f t transforms into peroxyl radical (P O * 2 ); thus the termination reaction will occur mostly by reaction (VI). It was reported in the literature that under steady-state assumption, the overall rate of oxidation in case of excess oxygen is linked to the constant value k 2 3 /k 6 [6]. On the other hand, in the case of limited oxygen availability, the overall rate can be linked to the constant k 2 k 6 2k 5 k 3 [P H] [23], [24]. Therefore, the chemistry of oxidation reaction is influenced by several factors and requires a thorough understanding at the macromolecular level in order to couple its effect with the macro-scale deformation kinematics.
3. Thermodynamically consistent chemo-mechanically coupled theory for polymer oxidation

Extent of reaction
Based on the work of high-temperature oxidation in metallic alloys [46] and a more recent work in curing of glassy polymers [58], we define the extent of reaction for each reaction (n) in a local dimensionless form as, where (n) is the number of reaction and can vary from 1 to 6 depending on the oxidation scheme for a specific polymer. Since polymer oxidation is a multi-component system involv- ing multiple reactions (as described earlier), changes in the chemical concentration for each component β are related to the extent of all the reactions it is participating as, where r (β) is the rate of production or consumption of chemical species β and R (n) is the stoichiometric coefficient in reaction n (measured in moles per volume), that gives the consumption or production rate for any species β.
d r a f t

Kinematics of the deformation
The notation used in this work is standard of continuum mechanics [45]. The function χ defines a map which converts every material point in the reference configuration X ∈ B R to spatial or current one as, x ∈ B t . A smooth deformation is a one to one mapping, x = χ(X, t) with deformation gradient defined as, and from the standard requirement of the local condition of impenetrability we can assume that, Oxygen embrittlement causes chain scission and weight loss in the polymer, combination of which induces shrinkage strain in the oxide layer. Treating the shrinkage strain as irreversible strain (like plasticity), we consider a multiplicative decomposition of the total deformation gradient based on the large deformation theory of polymers [4], where F e is the elastic component of the deformation gradient, F p is the irreversible plastic component, which incorporates the deformation due to permanent chemical shrinkage and viscoplastic deformation that can occur in the glassy polymer. Referring to Eq.6 the determinant J can be decomposed into, where J e = detF e > 0 and J p = detF p > 0 so that F e and F p are both invertible. The right polar decomposition of F e is given by, where R e is the rotation and U e and V e are symmetric, positive-definite stretch tensor with, d r a f t The left and right Green-cauchy tensors are then defined as, The velocity gradient is defined as, Substituting Eq.6 into Eq.11, we get With, The elastic and plastic deformation tensor can be further defined from the standard continuum mechanics as, So that L e = D e +W e and L p = D p +W p . Assuming plastic flow is irrotational, (as,W p = 0), the plastic velocity gradient becomes L p ≡ D p . We also make another kinematical assumption that the plastic stretch D p can be additively decomposed as, where D s represents the inelastic strain rate due to oxidation induced shrinkage and D vp represents the inelastic strain rate due to bulk viscoplastic deformation in the polymer.
Assuming chemical shrinkage to be isotropic and depends on the rate of oxidation, we have considered the following simple form for D s as, d r a f t where γ ≤ 1 is a material parameter that determines the amount of volume shrinkage. The viscoplastic flow rule for a glassy polymer is widely available in the literature and thus not discussed in this paper [2-5, 10, 42, 50]. Finally, as a standard, we assume plastic flow is incompressible, so that, and then, by Eq.7 we have, It can be shown that, F p , L p and C e are invariant under any change in frame. Detailed proof of frame-indifference can be found in many literature, such as, [2-4, 17-19, 46].

Mass balance for the diffusing and reacting species
Let P t is any arbitrary part of the deformed body B surrounded by the boundary ∂P t .
Let c (β) 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 (β) R for any species β are caused by diffusion of that species across the boundary ∂P t , which is characterized by a flux j (β) (x, t), defined by the number of diffusing species measured per unit area per unit time and a chemical potential of the species µ (β) . The diffusion of the species β into the polymer depends on the difference of the chemical potential between the boundary ∂P t and within the body B t and − ∂Pt j (β) .nda represents the number of moles of diffusing species entering the body B t across ∂B t per unit time. Similarly r (β) (x, t) denotes the number of moles for 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, d r a f t for every part P t , n represents a unit vector in the current configuration. Bringing the time derivative in Eq.19 inside the integral and using the divergence theorem on the integral over Since P t is arbitrary, this leads to the local mass balance as,

Specific form of mass balance for polymer oxidation
Based on the kinetics of oxidation reaction provided in the literature and discussed in section 2, we make some reasonable assumptions for mass balance of the reacting and diffusing species as given below, • Oxygen is the only diffusing species-during oxidation, oxygen diffuses into the polymer and reacts with polymer macromolecules. As a consequence of the oxidation reactions some volatile compounds are also being produced which eventually leaves the polymer. Therefore, during oxidation both inward and outward diffusion occurs.
However, in this work, we consider oxygen is the only diffusing species and neglect the outward diffusion of the volatile compounds.
• Steady-state reaction-as it is well known in the chemistry literature, polymer radicals are highly reactive in nature and thus they react almost immediately after formation. Hence we do not need to track their concentrations during mass balance.
• Extent of reaction-for each elementary reaction, we have introduced a state variable, extent of reaction, ξ (n) , where n indicates the number of elementary reaction. We consider that, these extent of reactions are directly related to the local reaction rate for each elementary reaction, n, by a stoichiometric coefficient, R n , so that the local reaction rate can be expressed as R nξ(n) .
d r a f t • Global reaction rate for species-we employ a similar approach as proposed in the kinetic model [26] to calculate the global reaction rate for each species. Following this, the change in rate for any given species concentration is the algebraic sum of elementary rates. Thus for any species β, the rate can be expressed as, With the assumptions mentioned above and using the closed loop reactions given in Eq.1, we consider the following: i) Oxygen is the only diffusing species.
ii) POOH is the only product during the propagation step.
iii) PH substrate consumption is worth tracking, only in the cases where propagation reaction (III) is dominating (i.e. in the cases of excess or constant oxygen concentration).

Balance of oxygen concentration
Applying mass balance, to track the oxygen concentration during the oxidation reaction at any time, t, we can write, As discussed in section 2, the oxidation reactions in polymers involve a multi-component system consisting of peroxide (POOH) decomposition and polymer substrate (PH) consumption contributing in the formation of alkyl radicals (P * ); which further takes part in multiple reactions. Among them, oxygen is consumed or produced in two of the reactions. Therefore the change in oxygen concentration, due to the reactions can be written as: where r (O 2 ) is the rate of consumption or production of oxygen. From Eq.1, oxygen is consumed in reaction (II) and produced in reaction (VI). Hence we can re-write Eq.23 using d r a f t Eq.24 as, where R 2 and R 6 are the stoichiometric coefficients of oxygen in reaction (II) and (VI), respectively. Note that, the stoichiometric coefficients for the consumed and produced O 2 are considered as negative and positive, respectively.

Balance of hydroperoxide, POOH and the polymer substrate, PH
As discussed in section 2, POOH is the main initiator of the oxidation reactions. The reaction is auto-accelerated due to the formation of this initiator further in the propagation step. Therefore, the balance of POOH concentration comes from the initiation reaction (Ia) and propagation reaction (III); we can write, where R 1 and R 3 are the stoichiometric coefficients of hydroperoxide in reaction (Ia) and (III), respectively.
In the case of excess oxygen, the concentration of the substrate, PH remains no longer constant. As mentioned earlier in section 2, the initial reaction (Ib), where substrate directly reacts with oxygen, is highly unlikely because of its higher activation energy. Therefore, we can neglect any changes of PH-concentration change occurring form (Ib). Then the rate of change of PH-concentration can be tracked from only reaction (III) as, where R 3 is the stoichiometric coefficient of substrate, PH in reaction (III).

Balance of forces and moments
Since the time-scale associated with oxygen diffusion is much longer compared to the scale associated with deformation wave speed, we neglect all forms of inertial effects in the present d r a f t formulation and assume that body forces does not depend on time. Therefore, the force and moment balance in the deformed body B t can be expressed as, where T is Cauchy stress tensor and b is the non-inertial body force per unit volume of deformed body. Similarly, for the reference body, there exists a stress tensor, T R , called the Piola stress, such that the surface traction on the surface ∂B R of reference body B R , is given by, where, t R (n R ) is the surface traction and n R is a surface unit vector.
and T R satisfies the local force and moment balance as, where b R is non-inertial body force per unit volume of the reference body. Finally, as standard, Piola stress is related to Cauchy stress in the deformed body by, and also as standard, T is frame-invariant.

Energy balance
Let us consider an isothermal polymer oxidation process. Neglecting the inertia effect, the internal energy in any arbitrary part P t of a deformed body B t can be written as less than or equal to the power expended on the part plus the energy increase in the system due to change in concentration in the species due to the diffusion and reactions. Here, we calculate the free energy contribution due to coupled diffusion-reaction based on the works in [34], [44], [17] and [46]. Denoting ψ as the internal energy density per unit volume of the deformed body, we can write the energy imbalance as, is the chemical potential for any species β and µ (β) j (β) is the energy carried by species β into part P t by the flux j (β) . Similarly, µ (β) r (β) indicates energy change in part P t due to chemical reactions. Applying divergence theorem to the terms on the integrals over ∂P t and bringing the time derivative inside the integral, Eq.32 can be re-written as, Using mass balance Eq.21, 22 and 28 in Eq.33 and considering P t is arbitrary, yields the local form of energy balance as,

Stress-power
Assuming the constitutive behavior of a generic glassy thermoset as elasto-plastic, we can write the stress-power as, From Eq.35 the following stresses have been defined: • The elastic second Piola stress, which is symmetric, as Cauchy stress T is symmetric.
• The mandel stress, which in general is not symmetric.
Using the fact that T is invariant, we can conclude T e and M e are also invariant under the change of frame. Taking time derivative of Eq.10, we can write, Now, let us assume ψ I is free energy per unit intermediate volume such that So that the energy imbalance Eq.34 can be rewritten as,

Basic constitutive equation
Guided by the free energy imbalance Eq.42, a functional form is assumed for the free energy ψ I ; the second Piola stress T e and the chemical potential of species µ (β) are then determined by the constitutive equations of the form, Further, the species flux j (β) can be assumed to have the following constitutive form: Then, by invoking the thermodynamic restrictions it is possible to generate the specific forms for the constitutive equations.
d r a f t

Thermodynamic restrictions
By Eq.44ψ In view of Eq.45, the free energy imbalance Eq.42 is equivalent to the requirement that the following inequality must be satisfied for all constitutive processes: For any arbitraryĊ e andċ (β) R to hold the above equation true, their coefficients must vanish. Therefore, the thermodynamic restrictions can be obtained through the 'state relations' as given below: i. The free energy determines the second Piola stress as, and the chemical potential as, ii. The species flux satisfies the species-transport inequality as, iii. Finally, the following dissipation inequality must be satisfied by the chemical reactions: d r a f t Thus from Eq.50, it is possible to define a force of dissipative nature F (n) which is a conjugate toξ (n) for each reaction n. Thus, where A (n) is defined as affinity in chemistry literature and defined by, This reflects the fact that every chemical reaction relates to a dissipative mechanism. The extent of reaction for any reaction n, ξ (n) is assumed to evolve according to a constitutive with Fξ (n) > 0 wheneverξ (n) > 0.

Fluid flux and Fick's law
In this work we assume that the fluid flux obeys Fick's law, i.e. flux depends linearly on the gradient of the chemical potential and thus can be expressed by, with M (β) (C e , c R , ξ (n) ) is the mobility tensor. The consequence of species transportation inequality 49 is that, the mobility tensor is positive-definite for any non-zero concentration

Specialization of the constitutive equations
We only consider the materials that are initially and continually isotropic. Under such assumption, it can be shown that the response of the polymer is invariant under any arbitrary rotation in the reference configuration or intermediate configuration and the constitutive responses do not change [2-4, 17-19, 46].
In order to apply the coupled theory presented so far, we develop a specific constitutive form in this section. For this purpose, we consider a simple NeoHookean elastomer for the d r a f t mechanical part of the free energy. Our model has multiple species reacting among themselves in multiple reactions during oxidation process. Hence, we consider that the specific form of the chemical part of the free energy takes a simple quadratic form of the extent of reaction ξ (n) , based on the literature [62]. Finally, the free energy for diffusion is considered as the energy of mixing between the species β and the polymer and depends on the chemical potential µ (β) .
The energy for diffusion or mixing is taken as to follow the Flory-Huggin's theory based on earlier works of [17][18][19]33].
The following section presents the details of free energy and other constitutive relations for a simple NeoHookean elastomer, undergoing an isothermal oxidative diffusion-reaction process.

Free energy
Motivated by the free energy representation in [46], a separable form of the free energy has been considered aŝ (i) The mechanical energy for a NeoHookean elastomer has the form where C dis = J −2/3 C e and G(ξ (n) ) is the oxidation dependent shear modulus, K is the bulk modulus of the polymer.
(ii) ψ chem I is the chemical energy influencing the oxidation reaction and is given by where the parameter H (n) is the chemistry modulus of reaction n. This term in the free energy favors the local state (ξ (n) = 1) [62].
(iii) ψ dif f I is the part of free energy comes from diffusion of species β into the polymer.
As established in section 3.2.1 , oxygen is the only diffusing species. Here we have d r a f t incorporated Flory-Huggins theory for mixing of polymer to define the energy due to diffusion [33]. It is given by the following form: where µ 0(O 2 ) is the reference chemical potential of oxygen, R is universal gas constant, ϑ is the oxidation temperature, χ is the dimensionless Flory-Huggins interaction parameter and Ω is the volume of a mole of oxygen.
Thus the complete free energy expression Eq.55 can be written aŝ

Oxidation dependent shear modulus
Following the work of [37], we consider that bulk modulus of the polymer does not depend on oxidation. However, shear modulus becomes higher because of the oxide layer formation [37]. As polymer radials P * reacts with O 2 to create P O * 2 in the propagation reaction (II), it is logical to assume that the change in shear modulus happens mostly due to the propagation reaction (II). Thus, we can write, where G un and G ox corresponds to the shear modulus of unoxidized and completely oxidized polymer, and ξ 2 is the extent of reaction for propagation reaction (II).

Stress, chemical potential and affinity
From the free energy Eq.59 and using guidelines for thermodynamic restriction, it is possible to get the specific constitutive equations for Cauchy stress T, chemical potential of oxygen, µ (O 2 ) and affinity of each reaction n, A (n) .
d r a f t Using Eq.60 , 47 and 36, the Cauchy stress can be expressed in the following form: where (B dis ) 0 is the deviatoric part of elastic Left Cauchy-Green tensor B e . Also, using Eq.48 into Eq.59, the chemical potential of oxygen can be written as, And the affinity of any reaction n can be calculated from Eq.52, 60 and 59 as, where n is the reaction number.

Reactive force and evolution of extent of reaction, ξ (n)
Applying thermodynamic restriction, we get a dissipative natured positive force F (n) , which is a force conjugate to ξ (n) . This simply means that for any reaction n to occur, we need a positive driving force, F (n) . Using Eq.50, 59 and 62, we can define the driving force for reaction as We assume a thermally activated relation for the evolution of the rate of extent of reaction, where, k n is the rate constant for reaction n, Q act is the activation energy, C ξ is a preexponential factor and has a unit of 1/(s − MP a).
d r a f t It is evident that the reaction will occur only when F (n) > 0. Also, it is logical to assume that there will be a threshold value of species concentration below which no reaction will occur. Denoting the threshold as c thr R , we can write that that the reactions occur only when available c (β) R > c thr R at a material point. Following this, we can rewrite Eq.65 as,

Oxygen flux and diffusivity
From Eq.54, oxygen flux is given by, with m (O 2 ) being the mobility of oxygen. Ignoring the effect of temperature gradient we can write from Eq.62, Further, diffusivity of oxygen can be defined as, and using Eq.69 in 67 the constitutive relation for species mobility can be obtained as, The diffusivity of oxygen in virgin or unoxidized polymer is different than the diffusivity in the active oxidatione zone or completely oxidized zone [59,60]. Also, diffusivity of the oxide layer is the reaction controlling parameter, as it determines the availability of oxygen in the core to take part into further oxidation [59]. Thus to comply with the literature, a simple form of diffusivity is assumed in this work that varies linearly with the extent of reaction as, where the parameters D un and D ox indicate the diffusivity of oxygen in the unoxidized polymer and in the oxide layer, respectively. Further, diffusivity follows simple Arrhenius rule, to incorporate the temperature dependencies as, where D 0,un and D 0,ox are the reference diffusivity values and Q d,un and Q d,ox are the activation energy of diffusion for the unoxidized and the oxidized material, respectively.

Governing differential equations and the boundary conditions
There are two governing differential equations required to be solved in this case: 1. The local force balance in the current configuration as, 2. Considering oxygen consumption in propagation reaction (II) is much more prominent compared to the termination reaction (VI), the local balance equation for the concentration of oxygen from Eq.25 in the current configuration can be written as, The initial and boundary conditions are required to complete the solutions of these differential equations. Let, S 1 and S 2 are complementary subsurfaces of the boundary ∂B of a spatial body B such that S 1 ∪ S 2 = ∂B and S 1 ∩ S 2 = ∅. Similarly, let S c R and S j are complementary subsurfaces of the boundary ∂B = S c R ∪ S j with S c R ∩ S j = ∅. For a time interval t ∈ [0, T ], two boundary conditions can be described such that displacement is known on S 1 and traction on S 2 . Thus, we can write, Similarly, another pair of boundary condition can be considered for a time interval t ∈ [0, T ] such that oxygen concentration is known on S c R and oxygen flux on S j and thus, The initial conditions are, u(X, 0) = u 0 (X), and c Thus the coupled set of equations, Eq.73 and 74 together with 75, 76 and 77 pose an initial boundary value problem to be solved for the displacement u(x, t) and concentration c R (x, t) simultaneously.

Numerical Implementation of the coupled multiphysics problem
In the absence of body force, the coupled differential equations reduce to divT = 0 in B; with u =ǔ on S 1 ; Tn =ť on S 2 Next, the weak forms of these initial boundary value problems described above are obtained by multiplying each of the equations by two weight functions w 1 and w 2 , respectively; w 1 vanishes over the boundary S 2 and similarly, w 2 vanishes over S j . The corresponding weak forms are given by Next, the body is approximated using finite elements as B = ∪B e . Assuming the nodal solutions are the displacement vector and the oxygen concentration, the trial solutions are d r a f t interpolated inside each element by, where A denotes the node of the element and N A is the shape function, u A and c A R are the nodal displacement and oxygen concentration respectively. The weight functions are interpolated using same the shape function as, Using Eq.80 and 81 into 79 yields the following system of equations: Following which the system of equations are solved using an iterative Newton solver, with the element level residuals are given in the form as, for displacement residual and for oxygen concentration, respectively.
The tangents for Newton solver are defined as, Considering no traction force working over the boundary S 2 , the index form of K AB uu can be written as, where u B k indicates the nodal displacement in direction k at node B. Using the identity F mn = δ mn + ∂N A ∂Xn u A m , such that ∂Fmn ∂u B k = ∂N A ∂Xn δ mk , it can be obtained as, Then Eq.86 can be written as, Next, the tangent related to change in concentration c R is calculated as, where ∆t denotes the difference between two consecutive time steps. Similarly, the last two tangents are calculated as, and The system of equations are solved numerically for each element by writing a user element subroutine in ABAQUS/Standard (2017) [1].

Representative numerical simulations
There are numerous material parameters involved in our theory because of the complex coupling of the diffusion-reaction phenomena. Following a detailed review of the literature, we have chosen some typical values for these constants which is realistic to most of the d r a f t and extent of oxidation based on that. Under such assumption, affinity and reactive force simply come from reaction(II) and have the following form: and where c R is the oxygen concentration, R and ξ are the stoichiometric coefficient and extent of oxidation in reaction (II), respectively. As mentioned before, we take stoichiometric coefficient as negative for consumption of the reacting species. Hence, for this single reaction case, mass balance equation for the concentration of oxygen has the form, d r a f t The polymer specimen is a 1mm by 1 mm 2D-block. Using symmetry, we simulate only one quarter of the sample. In our numerical simulation, we consider a 2D plane strain geometry as shown in Fig.1. We mesh the geometry using 2D plane strain UPE4 user elements. For the mechanical boundary condition, we consider symmetry along x and y axis.
For the chemical boundary condition, we apply zero oxygen concentration along AD and AB and a constant oxygen pressure ofp O 2 along BC and CD.
Finally we oxidize the specimen isothermally at temperatureθ. Oxygen sorpsion occurs at the boundary by Henry's law, as, c R = Sp O 2 , where S is the solubility of oxygen into the polymer and p O 2 is the partial pressure of oxygen; followed by oxygen diffusion which is controlled by the diffusivity parameters, D un and D ox . We simulate the specimen to study the effect of temperature by comparing the extent of oxidation on the element E located at the top surface as shown in Fig.1.

Effect of diffusion and kinetic parameters on oxidation
In this section, we study the effect of diffusion and other kinetic parameters on the extent of oxidation. To accomplish this, we use the same specimen described in Fig.1 with the same boundary condition as mentioned in section 7.1.
To check the effect of diffusivity on the extent of oxidation, first, we keep D un unchanged at 2.1 × 10 −12 m/s 2 and study the oxidation in two cases with D ox equal to 2.1 × 10 −12 m/s 2 and 6.3 × 10 −12 m/s 2 , respectively. Later, we increase D un to 6.3 × 10 −12 m/s 2 keeping D ox at 6.3×10 −12 m/s 2 . For each cases, we measure the oxygen concentration at a node near the core of the specimen (marked as node N in Fig.1). The comparison between the concentration of oxygen at node N for the three different diffusivity values is demonstrated in Fig.3. As the diffusivity of the polymer increases, the oxygen can penetrate quicker into the specimen; thus increases the oxygen concentration at node N at any given time.
The chemistry modulus, H, is a measure of the energy required for a reaction to happen.
The affinity for oxidation reaction is given by, A = H(1−ξ). Hence higher chemistry modulus generates higher affinity, and higher affinity leads to a higher reactive force, F . Further, since  Table 1.
The comparison between the extent of oxidation as a function of time for these two different H-values is shown in Fig.4(a). As we can see from Fig.4(a), ξ reaches to 1 approximately ten times faster due to a 10 times increase in H value.
Activation energy of a reaction has a completely opposite effect on the extent of oxidation.
Higher activation energy of a material indicates high barrier against the extent of oxidation, which means material is less prone to oxidation. We demonstrate this by increase the activation energy of of the polymer specimen from 100kJ/mol to 110kJ/mol for propagation reaction (II). From Fig.4(b), we can see a significant drop in the extent of oxidation because of 10kJ/mol increase in the activation energy.

FE simulation for the extent of oxidation in a bulk polymer sample
To demonstrate how the extent of oxidation evolves with time and oxide layer forms, in this section, we demonstrate a case of a uni-directional diffusion of oxygen into a block of polymer as shown in Fig.5. The specimen is a 1x1 mm block with a thickness of 0.5mm, which is exposed to oxygen on the top surface and all other surfaces are impermeable as Finally we let the sample oxidize for 200 hours at 150 o C and study the oxide volume fraction or extent of oxidation evolution at different times. In Fig.6 a), b) and c), we present the contour of the extent of oxidation in the bulk specimen after 38 hours, 100 hours and 200 hours of oxidation, respectively. We consider ξ = 0.99 representing complete oxidation of the polymer at a material point.
The visible identification of oxidation appears after 38 hours of exposure with a thin oxide layer formation on the top surface as shown in Fig.6a). A 220 µm thick oxide layer is formed on the top surface after 100 hours of exposure which continue to increase and reach a thickness of 520 µm after 200 hours. In Fig.6, the "gray" zone indicates that the material In this section, we model a randomly distributed carbon fiber reinforced composite specimen to predict the shrinkage strain developed due to accelerated oxidation and corresponding residual stress following a similar work reported in [35]. The geometry of the model is shown in Fig.7. The dimension of the composite sample is 1 by 1 mm with the fiber radius of 3µm.
We model the fibers as linear elastic and has the properties as E = 230 GP a, ν = 0.3. We model the polymer matrix as the user defined material properties given in Table 1. The composite specimen is oxidized at 150 o C under 2 bar O 2 pressure for 38 hours. For the chemical boundary condition, we apply a constant oxygen pressure of 2 bar on the left side of the specimen. We consider fibers are impermeable and do not take part in the oxidation process. For the mechanical boundary condition, we consider symmetry along x and y-axis.   We can see the maximum negative strain is developed in the less-dense fiber zone, which are shown in "blue" in the contour plot. These zones are mostly effected by oxidation and fiber-matrix debonding can occur in these zones because of the inhomogeneous shrinkage strain. Fig.8(b) shows contour plot of the corresponding Von-Mises stress generated in the matrix of the composite specimen. Note that Von-Mises stress is quite non-homogeneous in these less-dense fiber zones. This inhomogeneous stress distribution within the composite specimen can eventually lead to initiation and propagation of micro-cracks into the material, leading to ultimate failure of the specimen.    The model shows the capability of predicting the progressive oxidation font in the material for given ambient conditions. The numerical simulations show the usefulness of the model to simulate oxidation process in a randomly distributed fiber composites. In future work, we plan to validate our model with real experimental data for specific polymers. Once fully validated, the present model would be capable of reducing the empiricism in long-term life prediction for polymers due to oxidative aging.

Acknowledgment
The authors gratefully acknowledge several valuable discussions with Prof.