Modeling the oxidation of a polymer-derived ceramic with chemo-mechanical coupling and large deformations

To get a better insight into the coating behavior of a polymer-derived ceramic material, we model and simulate the diffusion, oxidation and reaction-induced volume expansion of a specimen without outer mechanical loads. In this macroscale approach, we use an oxidation state variable which determines the composition of the starting material and the oxide material. The model contains a reaction rate which is based on the change of the free energy due to a change of the concentrations of the starting material, the oxide material and a diffusing gaseous material. Using this, we model a growing oxide layer in a perhydropolysilazane (PHPS)-based polymer-derived ceramic (PDC), containing silicon filler particles. Within the mechanical part of the modeling, we use the Neo-Hookean material law which allows for the consideration of volume expansion and the diffusion kinematics in terms of finite deformations. We derive this continuum formulation in 3D and reduce it later to 1D, as we show that a 1D formulation is sufficient for thin oxide layers in our consideration. In such a case, the reaction-induced volume expansion is mostly limited to strains orthogonal to the oxide layer, as the bulk material hinders transversal deformation. Both formulations, i.e., 1D and 3D, are implemented in the finite element software FEAP. We perform a parameter study and fit the results with experimental data. We investigate the diffusion kinematics in the presence of volume expansion. Additionally, we discuss the influence of the elastic energy on the reaction rate.


Introduction
In high-temperature applications in mechanical engineering, the knowledge about the oxidation behavior is of great interest. The intrusion of oxygen into the surface of a machine component often enables chemical reactions which can change the material properties or the geometrical shape of the component. This can lead to a decrease of the performance or to a complete failure of the machine. Sometimes, materials which have a good resistance against oxidation do not have sufficiently good mechanical properties to withstand the loads. Also, the weight should be as small as possible. To overcome this problem, coatings are often used. Materials with a high oxidation resistance are applied on the surface of the machine components to shield it from admitting reactive gases.
One promising candidate for coatings is a polymer-derived ceramic (PDC) based on perhydropolysilazane (PHPS). A lot of research in characterizing PHPS has been done by Günthner and co-workers. In their early works, they investigate different material combinations, e.g., PHPS with boron nitride (BN) as filler particles [1] or the unfilled PHPS compared to a different PDC precursor which is a specific polycarbosilazane abbreviated as ABSE [2]. ABSE is synthesized for the article while PHPS is commercially available. In their experiments, the latter shows better protection. However, filler-loaded ABSE is investigated in an additional article [3]. The group kept on investigating both precursors and focused on the conversion behavior and the mechanical properties like the Young's modulus and the hardness of the unfilled material [4]. They also tested a barrier coating system based on PHPS as a bond coat and a further filler-loaded PDC (namely KiON HTT1800) as a top coat [5]. One advantage of such systems is that it can increase the critical coating thickness of PDC coatings. It measures the maximum thickness which can be used without cracking and delamination.
Recent research on this topic has been done by our group. In most cases, PHPS was loaded with boron (B) and silicon (Si) fillers. Oxidation tests for the coating were performed at 800 • C [6] and in another work at 800 • C and at 1100 • C [7]. In the latter, the effect of different filler content with the focus on boron was examined in nitrogen (N 2 ) and argon (Ar) atmospheres with promising results. In another work, Mo 5 SiB 2 (with molybdenum) powders as fillers were examined [8]. Even more recently, the effect of the Si and the B particle size and of the protection layer thickness was investigated [9]. In another article, our group studied the rheological behavior of the coating slurries and the thermal expansion of the coating which can be matched to the one of the protected molybdenum alloy-based substrate by varying the content and composition of the fillers [10].
The topic is also part of our review [11] on barrier coatings for molybdenum-based alloys. In general, besides PHPS, other precursors for PDCs are used for coatings; one example is poly(hydridomethylsiloxane), cf. Torrey et al. on its phase evolution [12] and on its processing [13]. In one of the mentioned articles, we performed oxidation tests for filler-loaded PHPS using silicon and boron particles [7]. The experiments are the basis for the paper at hand. Here, we introduce a reaction-diffusion model which is used to reproduce specific oxidation test results.
Often, the results of oxidation tests can be quantified by simple algebraic equations which relate the mass gain or the width of an oxide layer with the time. To model the oxidation in a way that also allows the simulation of possibly complex geometries, usually reaction-diffusion equations can be used. The oxidation behavior then can be obtained, e.g., by numerical integration of these partial differential equations. While the reaction rate can be expressed pragmatically by formulating dependencies on the educt concentrations, also energetically motivated reaction rates can be used. Including energy changes allows a simple way of incorporating different phenomena which influence each other in the oxidation process like reaction swelling which then can affect the diffusion paths. In the paper at hand, we mostly focus on energetically motivated models. Some examples for such models are Guyer et al. [14], Liang et al. [15] and Bazant [16] which all approach oxidation of metals from an electro-chemical point of view. The latter one also considers small elastic deformations. Other examples which include small deformation frameworks in their model are on reactions in cement [17] and on thermal oxidation [18].
When dealing, e.g., with polymer materials, large deformations are often considered. One work including large reaction swelling is done by Drozdov [19]. The article is about a model for the behavior of a polymeric gel. However, there are also works on high temperature coating materials using a large deformation framework, cf. the articles of Loeffel and Anand on the theory of their model [20] and its application [21]. They approach a well-established oxide/metal bilayer consisting of a bond-coat for oxidation protection and a top-coat for thermal insulation. They model the oxidation of the bond-coat and include large volume changes due to chemical reactions. However, often, the oxidation of high-temperature coating materials like ceramics or metals is treated within small deformation frameworks, e.g., one article on the investigation of the oxide layer growth in undulated interfaces [22], or one on the infiltration of calcium-magnesium-alumino-silicates [23]. Nakajima et al. also only consider small deformations when studying the stresses beside undulated interfaces [24].
A typical element in coatings is silicon. The reaction from pure silicon to silicon-dioxide (Si O 2 ) causes a volume change. For silicon, this volume change can be determined by the ratio of the molar volumes, V M i , i.e., the Pilling-Bedworth ratio, (1) Fig. 1 Schematic illustration of the oxidation process for a hexahedral specimen The molar masses, M i , differ by approximately factor two, while the densities, ρ 0 i , do not differ enough to compensate this. The ratio Φ is still nearly factor two which motivates a large deformation framework. The superscript zero indicates, that here, we include standard values, e.g., the density without any additional mechanical compression. For the material, Φ is the favored volume change; however, influences like the lack of space for a volume expansion can affect the actual change. This can be included in the mechanical material law which is shown later in this contribution.
When dealing with diffusion and large deformations, one has to choose if the diffusion is deformation dependent or deformation independent. Here, the latter means that the current shape of the body determines the diffusion behavior, regardless of any present deformations. This diffusion assumption is widely spread in the field of polymer modeling. Hong and colleagues use this variant in their theory for the modeling of elastomeric gels [25]. An FEM implementation of the same theory is done by Zhang et al. [26]. Duda et al. derive similar governing equations, being able to cover more possible boundary conditions [27]. But they use the same concept for diffusion at large deformations. Also with this diffusion concept, Bouklas and Huang compare linear and nonlinear theories on polymeric gels [28]. A more recent work is done by Ganser and colleagues on polymers within an electro-chemo-mechanical model [29]. Chester and Anand, in contrast, based the diffusion on the initial shape of the body for elastomeric materials [30] and also within thermo-mechanical coupling [31]. This neglects the current shape of the body for the diffusion. However, later, still for elastomeric materials, they changed their approach for the diffusion [32] following Hong et al. [25] and Duda et al [27]. Before this change, Anand together with Loeffel also based different works in the field of high-temperature coatings on the diffusion in the initial shape of the body. We already mentioned their works earlier. First, they developed the theory [20] and then applied it within numerical modeling of experimental results [21]. Di Leo, Rejovitzky and Anand additionally applied the diffusion in the initial shape of the body to lithium-ion electrodes [33], also before the mentioned change. A different group applied it to silicon carbide fibers [34]. Neff et al. use some kind of hybrid in their model on rubber by basing the diffusion flux on the initial shape, but consider influences of volume changes rather than both, volumetric and isochoric changes [35]. There is no general agreement for the choice of the diffusion approach, and there are not many works on diffusion at large deformations in high-temperature coating materials. For our model, we choose the deformation-independent diffusion flux, i.e., we consider the current shape of the body.
In this paper, we approach the modeling of oxidation at 1100 • C for PHPS with silicon fillers. For this, we formulate the diffusion of a gas and we formulate the reaction with the starting material to oxide. In this diffusion-limited process, an oxide layer forms which grows over time due to further oxidation, cf. Fig. 1. The reaction triggers volume expansion which leads to an additional increase of the oxide layer width. Thus, the distance which the gas, i.e., oxygen has to travel to reach the reaction zone increases over time. The volume expansion is incorporated in a Neo-Hookean material law. Besides this, also the diffusion and the reaction rate are based on a free energy function.
The outline of this article is as follows. We first introduce the basics of the model in Sect. 2 including the balance equations. Then, in Sect. 3, we provide the constitutive relations based on specific free energy formulations. We summarize the final system of equations and prepare them for a finite element method (FEM) implementation, cf. Section 4. In Sect. 5, we first motivate a model reduction from a 3D formulation to a 1D formulation, and then we approach a fit of the model to experimental data. This is followed by the investigation of the influence of the reaction swelling on the oxidation process and finally a parameter study. Section 6 is the summary and the outlook. In this section, we provide the basics for the mechanical part of the modeling, including basic equations like the balance of linear momentum. We also provide the framework for the chemical reactions and the evolution equations of the chemical species. For this, we also introduce an oxidation state variable.

Deformation
We base the modeling of the chemical swelling on the widely used multiplicative split of the deformation gradient with the elastic part F el and the inelastic part F in . Especially in plasticity, the split is often termed after Kröner [36] and Lee [37]. We use it to describe the chemical reaction-induced volume change via where λ in is an inelastic stretch. Such a spherical shape of F in is a common choice when describing swelling or volume changes in general. The origin of multiplicative splitting the stretch for such applications is sometimes attributed to Flory in his work on swelling polymer networks [38]. For the determinants, the split subsequently gives J = J el J in and J in = (λ in ) 3 . We illustrate this split in Fig. 2 for a reaction-induced swelling in combination with an elastic deformation. As another example, imagine a chemical reaction which forces the body to swell. When there are no outer mechanical restraints, the body actually swells; however, when there is no space for swelling, F = 1 is constant, and we can calculate F el = F(F in ) −1 = 1. The imposed reaction swelling F in increases F el which then increases the potential elastic energy of the system.
Other deformation gradient-based quantities that we use are the left and the right Cauchy-Green tensors, the Green-Lagrange strain tensor and the isochoric part of the right Cauchy-Green tensor These quantities, in combination with a superscript index el or in, yield such quantities for the elastic part or for the chemical swelling part, respectively. The quantities will be used as strain measurements and in our constitutive law. Note that for the isochoric part of the right Cauchy-Green tensor, we haveĈ =Ĉ el , since F in solely contributes volume changes. The local balance of linear momentum for the absence of body forces and when neglecting inertia is div(σ ) = 0 with the Cauchy stress tensor σ . In this paper, it is the basic differential equation for the mechanical part of the model. It is one of the three differential equations which need to be solved to describe the whole coupled model for oxidation in this article.

Concentrations and chemical reactions
The matrix of the investigated material consists of a polymer-derived ceramic of the polysilazane type. The bulk samples were produced using a doctor blade procedure applied to a mixture of 45 vol.% PHPS and 55 vol.% silicon fillers. The obtained samples were subjected to pyrolysis at 1000 • C in N 2 for 1 hour. The further oxidation of the samples was conducted in air at 1100 • C. The details of the production procedure and the cyclic oxidation test are described in Smokovych et al. [7]. As a result of the polymer-to-ceramic conversion, the obtained composites were formed on the base of Si 2 N 2 O (silicon oxynitride), Si O 2 and elemental silicon according to the following reaction: Additional free elemental silicon results directly from the PHPS after pyrolysis. Its interaction with oxygen during the cyclic oxidation procedure causes the constantly growth of the Si O 2 oxidation barrier layer: For our model, we assume that the oxygen determines the rate in the oxidation process. Whenever new oxygen reaches the reaction zone, the oxidation of the specimen progresses to form Si O 2 . We consider the Si 2 N 2 O as an intermediate state. We use the representative reaction equation starting material + gas → oxide.
The stoichiometric coefficients for the starting material (index s), the gas (index g) and the oxide (index o) are s s = −1, s g = −1 and s o = 1, respectively. For now, we neglect any diffusion and use the stoichiometric coefficients to writeċ Here, c i are current placement concentrations of the species i, i.e., the amounts of moles related to the current placement volume. Equation (9) states that a chemical reaction with the reaction rate r increases the concentration of the oxide and decreases the concentrations of the gas and the starting material by the same number of moles per volume. Later in this article, we model r based on the free energy. This then determines the evolution of the species concentrations, c i . In the following section, we include diffusion into the system.

Conservation of mass
The reaction-diffusion equation is derived from the global mass balance. The total mass of the domain B, occupied by the current placement densities, ρ i , for the number of n s species is The change of mass in this domain is equal to the mass fluxes, f i , over the boundary ∂B. The temporal change of the mass,ṁ = dm dt , then iṡ The normal vector n of the surface increment points outward which yields the negative algebraic sign in the flux term. The divergence part has its origin in temporal changes of the volume, i.e., reaction swelling. Chemical reactions can also influence the system by consuming specific constituents and by producing other constituents. For this, we introduce a source term which can change the local mass concentrations for every species. A chemical reaction does not change the total mass of a system; thus, the sum of all mass changes should be zero, with an arbitrary reaction rate r for a single chemical reaction. Each coefficient a i is introduced to construct the expression a i r which quantifies each mass change due to r . To ensure conservation of mass in the reaction, the sum of these changes has to be zero. Since it is zero, we can add this sum to Eq. (11). The global balance of mass then is Using the divergence theorem and some rearrangements, the local formulation for each species becomeṡ When this expression is satisfied, the global formulation is also satisfied. As stated in Eq. (12) and due to the arbitrariness of r , the sum for all involved a i has to be zero, Dividing Eq. (14) by the molar mass M i gives the formulation in terms of molar concentrationṡ Here, s i = a i M i is the stoichiometric coefficient of species i in the reaction, cf. Sect. 2.2. The concentration and the density are related by c i = ρ i M i . Here, we call the quantity j i = f i M i the diffusion flux of species i.

Constitutive modeling
In the following section, we provide the connection between the balance equations and the Helmholtz free energy which we use to model the process. Then, we derive the expressions for the free energy. This yields the formulations for the stress, the reaction rate and the diffusion flux.

Thermodynamic restrictions and energy correlations
We base our argumentation on a standard dissipation inequality, Here, ψ R is the Helmholtz free energy in reference placement and S is the second Piola-Kirchhoff stress tensor. The quantities μ i are chemical potentials, C i = J c i are concentrations in reference placement, and J i are reference placement diffusion flux vectors. The sum over i represents three species with the indices v, g and o.
In such a dissipation inequality, typically, all terms except the ones which contain the fluxes are considered to be zero in sum. Then, to satisfy the remaining part, we can formulate for each species The quantities M i are positive semi-definite mobility tensors. We assume that mass transport which causes the oxide and the starting material to rearrange, underlies the relation J o = −J s . Increasing the oxide concentration decreases the starting material concentration by the same amount. From this and with Eq. (18), we can formulate In our model, we use μ o = −μ s which satisfies this relation for equal mobilities. While with J o = −J s the solid diffusion fluxes do not change the total amount of the solid substances, we can formulate something similar for the chemical reaction: In Eq. (8), 1 mol of the starting material reacts to 1 mol of the oxide, i.e., the chemical reaction also does not affect the amount. For a constant reference placement volume, we can write Compared to the gas diffusion, the solid diffusion does not play a major role, i.e., for the rest of the paper, J g is the only flux which we consider. The reference placement flux can be pushed forward to current placement, We use Eq. (21) with an isotropic m g =m g 1 which has the scalar mobilitym g . This formulation gives deformation independent diffusion behavior, i.e., the diffusion depends on the current shape of the body, regardless of the state of deformation.
With the relation for the concentrations, C s = C max − C o , and using standard arguments, from Eq. (17) we can extract where the factor 1 2 is due to μ s = −μ o andĊ s = −Ċ o since this allows to use μ sĊs = μ oĊo as a replacement in Eq. (17). For convenience, we choose the Cauchy stress tensor for our model, Another important quantity which was mentioned in the previous Sect. 2.3 is the reaction rate r . For now, we simply describe it as a function of the chemical potentials of the species which are involved in the reaction. It is involved in the evolution of the oxidation state variable, i.e., Eq. (28), and in the reaction-diffusion equation for the gas, i.e., Eq. (16).

Oxidation state variable
Based on the maximum concentration C max from Eq. (20), we can formulate a dimensionless quantity which we call the oxidation state variable, It ranges from zero to one where η = 0 indicates the complete absence of oxide in a material point and η = 1 indicates pure oxide in this point. We use η to track the evolution of the oxidation state variable which then can be used to calculate the concentrations C o and C s or c o and c s in the post-processing. Note that introducing η is not necessary, the simulation could also be based directly on the concentrations. However, the dimensionless η leads to slightly more compact expressions, and for tracking the progress of oxidation, this relative variable is more meaningful compared to an absolute concentration. Using Eq. (16) for the oxide and neglecting its diffusion, we can writeċ where we performed a pull-back in reference placement for the second one. The evolution of the oxidation state variable then isη Note that we prefer the deformation independent C max over c max since the reference placement version can be used as a constant input parameter.

Elastic energy, Cauchy stress and elastic constitutive tensor
For the elastic energy density, we choose a Neo-Hookean material with the common isochoric-volumetric split. The energy density reads Here, G is the shear modulus. We write it as a function of the material quantityĈ el instead of the spatial quantityb el as we derive the stresses and the fourth order material tensor in reference placement and transform them to current placement afterward. For scaling the volume expansion due to oxidation, we use the oxidation state variable from Eq. (26). We model the reaction induced swelling as follows. A reaction from starting material to oxide causes a volume change of J in = Φ. When there is no reaction at all, i.e., the composition at a point solely consists of starting material and J in = 1 remains. Using the oxidation state variable η, we can interpolate this via In the volumetric part of the elastic free energy density function, f (J el ), the equilibrium state is demanded to be at J el = 1. With J el = J J in , we can rewrite this as with the abbreviation β. Any deviation from the equilibrium state, i.e., from Eq. (31), increases the elastic energy. For this, we write and use the bulk modulus K to scale the increase of the energy. Furthermore, we insertedĈ el =Ĉ. A similar equation can be found in the work of Hong and Wang [39]. As mentioned before, the partial derivative of the total free energy with respect to C gives the second Piola-Kirchhoff stress tensor. Note that in this model, the elastic energy is the only contribution of the free energy which depends on C, i.e., other terms which are introduced later in Sect. 3.6 do not need to be included in the stress calculation. The Cauchy stress tensor, which we are more interested in, can be calculated based on Eqs. (22) and (25). Additional calculation steps are shown in Appendix A. The Cauchy stress then is and the hydrostatic pressure is Note that for a deformation-free state with constant b = 1, i.e., when there is no space for swelling, a chemical reaction causes a pressure of Kβη. The higher the concentration of the oxide, the higher the pressure in the case of lack of space. An increase of the pressure or the energy of the system due to a chemical reaction can influence the reaction kinetics which we can observe in the chemical potential. Here, the elastic contribution for the oxide is, For the FEM implementation which utilizes the Newton-Raphson method, we need the stiffness-tetrad related to the Cauchy stress tensor. For this, we apply a push-forward on the second derivative of the free energy with respect to the right Cauchy-Green tensor, cf. Appendix A. The current placement stiffness tetrad for Eq. (33) then is Here, we make use of the fourth-order identity tensor I = δ ik δ jl e i ⊗ e j ⊗ e k ⊗ e l .

A 1D representation for blocked transversal strains
Besides the general 3D formulation, we are also interested in the case that the swelling mostly results in a strain normal to the oxide layer. For this, we can approximate the behavior by assuming zero transversal strains. For this approximation, we use the e 1 direction as the only direction which has nonzero strains. The corresponding deformation gradient then is This gives the expressions For the 1D FEM implementation, we only need the σ 11 and the C 1111 component of Eqs. (33) and (36). Using Eqs. (38)(39)(40), we get and

Parameters for the elastic law and the reaction induced volume change
In experiments, Günthner et al. [4] determined the Young's modulus of their PHPS specimens, pyrolyzed at 1000 • C in nitrogen, to be about E PHPS = 150 GPa. The Young's modulus of pure silicon, which is anisotropic, usually has the values from 130 to 188 GPa, depending on its orientation [40]. However, the cubic diamond crystalline silicon particles in the PHPS-derived matrix are distributed and orientated quite randomly. Thus in our model, the overall Young's modulus of the material is assumed to be isotropic. For this we use E s = E PHPS = 150 GPa. Additionally, in our work, we use the same Young's modulus for the oxide phase. The Poisson's ratio also varies strongly in the anisotropic silicon [41]. Similar like Konegger et al. [42] who approach a different polysilazane-derived ceramic, we use a Poisson's ratio of 0.2. Using it together with the Young's modulus, we can calculate the shear modulus, G = 62.5 GPa, and the bulk modulus, K = 83.3 GPa for the equations provided in this section. The desired volume change Φ can be calculated based on the Pilling-Bedworth ratio. For starting material and oxide, it is The molar masses, which we involve here, are M s = M Si = 28 g mol and M o = M Si O 2 = 60 g mol . In the equation, we use a zero as a superscript which, here, indicates that we use the numeric value for mostly pure and macroscopically undeformed materials, i.e., standard values from textbooks or research articles. These densities are no evolving field variables in our oxidation model, in contrast to the densities without the superscript zero. We treat the densities ρ 0 s and ρ 0 o as material parameters. To calculate the mechanically unconstrained volume change from pure starting material to pure oxide, we use the values of the density of the silicon filler particles. The density of pure solid crystalline silicon at 1100 • C is ρ 0 Si = 2.32 g cm 3 [43]. For amorphous silicon-dioxide at 1100 • C, it is ρ 0 Si O2 = 2.21 g cm 3 [44]. While the ratio of the molar masses has approximately a value of two, the densities, in contrast, are quite similar, i.e., its ratio is approximately one. The reaction of one mole of starting material to one mole of oxide approximately doubles the mass. However, this would also double the density, as long as there is no volume change. In our implementation, we thus use

Chemical energy and chemical potentials
In our model, the chemical free energy is or in terms of the oxidation state variable Here, g f i is the molar energy of formation of the species i. Multiplied with the concentration, the term represents the energy of formation per volume with the contributions of each species, i.e., of the starting material, the oxide and the gas. The gas constant is denoted as R, while T is the absolute temperature. We model the system as a binary mixture of the two solid materials, while the dissolved diffusing gas is treated differently. For this, the second term in Eq. (44) is based on the entropy of mixing of the starting material and the oxide. Expressions like that can be derived from Boltzmann's equation from statistical thermodynamics [45]. We use the last term to describe the free energy of the diffusing gas. The shape of the term is motivated by the work of Loeffel and Anand [20]. In our model, the parameter C b g is the boundary gas concentration of the system. The chemical free energy in combination with the elastic free energy, given in Eq. (32), yields the total free energy ψ R = ψ chem R + ψ el R . The chemical potentials are

Modeling and investigating the reaction rate
We use the following reaction rate which is based on the derivation in the article of Pekař [46], We define it in the current placement. Note that the expression in the square brackets is dimensionless, while r and k have the dimension of a concentration per time. With the chemical potentials, the rate is The first term represents the forward reaction, the second term the backward reaction. In both terms, there are dependencies on the oxidation state. A state with η = 0 or η = 1, i.e., pure starting material or pure oxide causes an infinite reaction rate driving the system to η = 0 or η = 1, respectively. The ratio scales the forward reaction. It contributes values from zero to one and represents the gas concentration. A higher concentration of the gas yields a higher forward reaction rate. Further scaling of both rates is induced by the additional exponential factors. For illustrative purposes, in combination with Fig. 3, we introduce which, here, we call the trimmed reaction rate. Note that in Eq. (51), we assume a constant maximum presence of gas concentration, i.e., In the following paragraphs, we discuss the behavior of the reaction rate in Eq. (50) by using the simplified one in Eq. (51) with different definitions of α 1 and α 2 . The hydrostatic pressure from Eq. (34) is p = −K (J − 1 − βη). As mentioned before, the special case of having a lack of space for volume change, i.e., J = 1 gives p = Kβη. In the reaction rate equation, the presence of positive hydrostatic pressure decreases the forward reaction rate. In turn, the backward reaction rate is increased. With a lack of space for reaction swelling, a state with much oxide is energetically less favorable.
In the plots on the left-hand side of Fig. 3, we plot Eq. (51) with α 1 = exp(−zη) mol cm 3 s and α 2 = exp(zη) mol cm 3 s using z = {0, 1, 2, . . . , 10}. We choose these specific shapes of α 1 and α 2 because they then have the same kind of dependency on η as the factors in Eq. (50) when inserting p for p. The quantity z substitutes the combination of constants multiplied by η. Here, z = 0 gives α 1 = α 2 = 1 mol cm 3 s which results in the same evenly large forward and backward reactions. Then, the equilibrium state, i.e., when the reaction rate becomes zero, is at η = 0.5. Increasing z which represents increasing the elastic contribution shifts the equilibrium in the direction of η = 0. Here, oxidation increases the elastic energy und thus hinders the forward reaction and accelerates the backward reaction.
In our actual reaction rate in Eq. (50), the hydrostatic pressure becomes very large, since no dissipative effects are implemented for the deformation. To incorporate the elastic contribution in a reasonable way in the model, one needs to include such dissipative effects first. However, this would be beyond the scope of this article. We decide to drop the elastic contribution from the reaction rate for the whole paper by setting p = 0 in Eq. (50).
We now discuss the influence of the molar energies of formation in our formulation with g , the enthalpy of formation, the absolute temperature and the entropy of formation. We take values from a reference book [47]. Furthermore, besides a temperature of ϑ = 1100 • , for this, we use the values for solid silicon, gaseous oxygen molecules and solid silicon dioxide. The molar enthalpies of formation of silicon and Based on this, we can neglect the back reaction term in Eq. (50). In the diagram on the right-hand side of Fig. 3, we plot Eq. (51) for α 1 = 1 mol cm 3 s and α 2 = exp(−z) mol cm 3 s with z = {1, 2, . . . , 6}. The higher z, the less significant the contribution of the α 2 term. Neglecting the back reaction term gives the thick red line with α 2 = 0. This qualitatively represents the reaction rate of Eq. (50) excluding back reactions. When redefining k to include the constant energies of formation, it reduces to Note that the reaction rate in Eq. (53) drives the system to pure oxide, i.e., η = 1 instead of just extremely close to η = 1.
The reaction rate in Eq. (53) becomes infinity for pure starting material, i.e., η = 0. In the numerical implementation, we approximate this by replacing the rate for small η with a finite value. For values smaller than η = 0.05, we use the constant value r = r (η). Furthermore, for values larger than η = 0.95, we use a linear interpolation between the corresponding rate r = r (η) and r = 0 at η = 1, since the slope of Eq. (53) is infinity for η = 1.

Modeling the diffusion flux
We only consider the diffusion of the gas. We drop the index for the gas diffusion flux, as no other flux is used. In the considered oxidation process, an oxide layer forms at the surface and keeps growing over time. For this, starting material is consumed in the chemical reaction. The gas molecules have to travel a longer distance to reach the reaction zone. Additional to this growth of the oxide layer, the oxide layer swells mechanically. This increases the distance furthermore. The latter effect is included by defining the diffusion flux with an isotropic formulation in the current placement, cf. Sect. 3.1, The corresponding chemical potential only depends on the gas concentration, i.e., the flux in terms of the gradient of the concentration is Here, D is a diffusion coefficient. Often, the mobility depends on the material composition. In the article Nauman and He [48], the diffusion of a binary mixture is modeled including one diffusion flux for each of the two components. Each mobility depends linearly on the mole fraction of the diffusing component. We assume a mobility which depends linearly on the gas concentration and linearly on the oxidation state. The mobility and the diffusion coefficient arê The diffusion flux drives the system to develop a homogeneous μ g field in the whole domain. In our specific formulation, in which the gas chemical potential solely depends on the gas concentration, this is equivalent to developing a homogeneous gas concentration field. At the same time, the chemical reaction consumes all of the gas as long as there is starting material available.

Summary of equations and comments on the solution procedure
We summarize the primary variables and the important equations in Table 1. Three differential equations have to be solved, namely, the balance of linear momentum, Eq. (5), the evolution of the oxidation state variable, Eq. (28), and the evolution of the gas concentration, Eq. (16).
Besides the obvious couplings of the oxidation state variable, the stresses and the gas concentration, another coupling is more hidden: The balance of linear momentum generates the displacement field and thus determines the spatial coordinates which are used in the spatial partial derivatives of the concentrations, i.e., the deformation changes the kinematics of the diffusion.
As we approach an implementation in a finite element method software, we are furthermore interested in the weak form of the problem, Here, grad S (w) = sym(grad(w)) is the symmetric part of the gradient. The scalar and the vectorial testfunctions are w and w, respectively. The right-hand sides of the equations contain quantities to implement the traction boundary conditions with t = σ n and the diffusion flux boundary conditions. We implemented this formulation as a user element in the finite element software FEAP. We use linear shape functions, i.e., we created an 8-node element for the 3D-case and a 2-node element for the 1D-case. To solve the discretized nonlinear equation system, the Newton-Raphson method is used. For this, based on the weak forms, Eqs. (58-60), and the shape functions, we can construct the elemental residual vector and tangent matrix which both are given for the 3D-case in Appendix B. For all simulations, we chose the Newmark method for time integration with equidistant time discretization.

Discussion
In this chapter, we first discuss a 1D simplification for the 3D-model. Then, we approach a fit to experimental data. We discuss the influence of the reaction swelling on the oxidation process and we perform a parameter study for different parameters. Gas concentration Differential equations for the primary fields div(σ ) = 0 Local balance of linear momentuṁ η = Jr Cmax Evolution of η c g + c g divẋ = div D grad c g − r Evolution of the gas concentration Additional equations  In this section, we want to legitimate a simplification that we use in our FE model for an experimental fit and a parameter study. Numerically, 3D simulations are extremely costly, so it is usually a good idea, to try to reduce the model to 2D or even to 1D. By doing so, of course, one has to keep in mind, that the reduced models should approximate the 3D simulations sufficiently well. When reducing to 1D, there are two central aspects which have to be considered. The first one addresses the reaction-induced volume expansion which results in large strains normal to the surface and quite small strains in transversal directions. The smaller these transversal strains are, the better is an approximation assuming zero transversal strains. This can allow the neglection of two of the three dimensions for the mechanical modeling. The second aspect relates to the diffusion. The diffusion behavior near to the surface is different when the surface is plane compared to regions next to edges where diffusion fluxes superimpose. When the influence of the edges is sufficiently small, the model can be reduced to 1D from the diffusion point of view. Both considerations do not require the exact specimen dimensions which we use later in this article. Furthermore, some parameters can be changed which we do for simplification. First, we investigate the transversal strains due to the volume expansion. For this, the correct Poisson's ratio and the correct Pilling-Bedworth ratio are important. The Poisson's ratio is implied by the shear modulus and the bulk modulus; therefore, we stick to the values for Φ, β, G and K which we described in the previous sections. Additionally, we use the already introduced η to clip the reaction rate to prevent infinite values, cf. Sect. 3.7. All other parameters are chosen to control the size of an oxide layer keeping the numerical effort small and allowing a good illustration of the effects. The parameters are listed in Table 2.

Starting material concentration
In Fig. 4, we illustrate a cube of the size 4mm× 4mm × 4mm. Three of the six faces are exposed to the gas and the other three faces are mechanically constrained in a way that forms symmetry planes, i.e., displacements  Fig. 4. Shown are the oxidation state variable and different Green-Lagrange strains of a 32×32×32 elements simulation (top) and of a 64×64×64 elements simulation (bottom). The oxide layers are nearly identical, while the percentage of the bulk material differs in normal direction are inhibited. The actual FE model thus represents 1 / 8 of the larger cube and provides insight into the reaction-diffusion process. In the later used 1D representation we aim to only investigate a line along the x 3 -axis which we also use for line plotting here in the 3D investigations.
Later in this paper, we will use small oxide layers and a big amount of bulk material. When we have a look at the detail which shows the mesh in Fig. 4, we observe that the initially perfectly cube-shaped elements mostly deform in the orthogonal direction of the oxide layer. This motivates the following consideration. We compare two cubes with different dimensions, i.e., 2 mm×2 mm×2 mm and 4 mm×4 mm×4 mm, but with nearly the same oxide layer thickness. Thus, effectively, the ratio of the oxide layer and the bulk material is different in the two cases, cf. Figure 5. To get the same quality of discretization, especially the reaction zone between oxide layer and starting material should be discretized equally in both simulations. This is easily achieved by using the same element edge lengths of 0.03125mm for both simulations. For the same size of oxide layer and reaction zone, this gives the same discretization. Using the same element sizes for different cube sizes gives in our case 32×32×32 elements in the small cube and 64×64×64 elements in the bigger one.
In Fig. 5, different profile plots along the x 3 -axis in Fig. 4 are plotted. The plot is extended in the negative x 3direction, undoing the reduction which uses the symmetry planes. This results in a profile plot along the full cube which ranges from oxide layer to bulk material to oxide layer on the other side. The curves in Fig. 5 show us different effects. The E 33 -strain is the normal strain coaxial to the x 3 -axis in Fig. 4. The E 11 and E 22 are normal strains orthogonal to it. As mentioned in the introduction of this section, we want to validate a 1D Fig. 6 Illustration of oxide layers of different thicknesses (h j ) and the segments (with length l j ) which are expected to be not significantly influenced by corner effects, such as the superimposition of diffusion fluxes model which should represent the 3D model. The 1D model assumes that the E 33 strain is much larger than the transversal strains, E 11 and E 22 . They are set to zero, similarly like in the plane strain concept.
In the graphic we can observe that a smaller oxide layer relative to the bulk material size gives smaller transversal strains. Additionally, the shape of the transversal strain curve becomes flatter along the whole axis. In the middle, at x 3 = 0, we have E

= E
(2) 22 = 0.0899. While here in the small cube, the transversal strain is 1 / 4 of the strain in normal direction, in the large cube, it is just 1 / 10 of the strain in normal direction. We expect that even smaller oxide layers relative to the bulk material size would lead to even smaller transversal strains.
For plotting along another axis, e.g., x 1 , the labels of the transversal strains would be E 22 and E 33 . A consistent consequence is E 11 = E 22 = E 33 in the middle of the cube. That there are negative strains right behind the oxide layer is logical, because the swelling layer which encloses the whole cube naturally tends to expand in both directions, x 3 and negative x 3 , compressing the inclusion.
In our following simulations, we have oxide layers from 10 μm to 300 μm. The smallest dimension of the actual specimen is 2.5 mm= 2500 μm. We expect that in these cases, the transversal strains are sufficiently small. Also note that an oxide layer of 300 μm thickness after 100 h, most of the simulation time, has a smaller size, since the layer initially has 0 μm and grows over time.
Near the edges, the oxidation behaves differently as the oxygen supply is bigger compared to the rest of the surface. When using a 1D representation in the simulations, this effect is not considered. The smaller the effect which is neglected, the better the simulation via 1D. In Fig. 6, we illustrate different thicknesses h j of oxide layers and the lengths l j of the parts which are not much influenced by the corner geometry. According to Fig. 6, a thinner oxide layer usually gives a bigger length l j . The thinner the oxide layer in comparison to the edge lengths, the better the 1D approximation. Additionally, in Fig. 6 (right), the radius is changed for the different layer thicknesses. Since usually, thinner layers have smaller radii, the influence is reduced even more. In our work, the proportion which is influenced by the corners and edges is expected to be neglectable. Due to numerical efficiency, we use the 1D model in all following simulations.

Fit to experimental data
In Fig. 7, we show a fit to experimental data. The data are taken from oxidation tests at 1100 • C for 100 h. The specimen has the dimensions 25 mm × 5 mm × 2.5 mm. The plotted curve on the left-hand side is the mass gain per exposed surface over time. In the performed oxidation test, the exposed surface can be approximated quite well by the whole surface of the specimen. More details on the experimental setup can be found in Fig. 7 Oxidation test with simulation results and fit to experimental data. Shown are the mass gain over time (left), the profile of oxidation state, η, and the normalized gas concentration, c g /c b g , after 100h (upper right) and after 50h (lower right) in the simulation Table 3 Parameters which we used for the fit and which are different from Table 2 Parameter Value Unit Comment There are not much data for diffusion coefficients in PDCs in literature. There is a work which uses diffusion coefficients of the magnitude of around 10 −4 cm 2 s in a PDC near room temperature [49]. Using polysiloxane as a precursor, their material is not the same as ours which is based on polysilazane. Another difference is their higher porosity. We use their material parameters as rough orientation while keeping in mind that the actual diffusion coefficient at the elevated temperature is much higher, i.e., several powers of ten. Silicon is known to have a really small diffusion coefficient for oxygen (D Si O 2 = 10 −10 cm 2 s at 1100 • C [50]). We expect that most of the diffusion is restricted to the PDC matrix which surrounds the filler particles and makes up 45 vol. % of our composite material.
We base the fit on another rough orientation which we use for the boundary gas concentration. Itoh and Nozaki [50] investigated the solubility of oxygen in silicon at different temperatures. At 1100 • C, an oxygen concentration of (3.5 ± 0.5) × 10 17 at cm 3 was determined at zero depth of the material. Using the Avogadro constant, we get approximately 5.8 × 10 −7 mol cm 3 as the boundary gas concentration. However, we do not model pure silicon or silicon dioxide, i.e., our c b g differs in the fit. For the maximum concentration of the solids, we use C max =  Table 3. On the right-hand side in the figure, we show a profile plot corresponding to the mass gain curve at 50h and at 100h. We can observe the oxide layer, a reaction zone and the unreacted core of the specimen. The gas concentration is a linear function in the oxide layer, followed by a smooth nonlinear section in the reaction zone where it tends to zero as it is consumed.
The initial simulated length is 0.04 cm which is not the full thickness of the specimen. In this 1D simulation, the influence of the size of the bulk material is already included and does not need to be modeled. The length is increased by swelling of the oxide phase and then it is no longer 0.04 cm but 0.054 cm. Furthermore, a larger

Influence of the reaction swelling on the oxidation process
In Fig. 8, we illustrate the influence of including the volume expansion. In our considered representative chemical reaction, one mole of starting material always reacts to one mole of oxide material. Thus, the concentration, i.e., moles per volume, is the same for the oxide layer and the starting material. However, one mole of oxide is heavier than one mole of starting material. When there is no volume change, the density increases. The total density is which we plot in the bottom right corner. We can observe that the density in the oxide layer is unrealistically high for the case excluding volume expansion. In the Pilling-Bedworth ratio, cf. Sect. 3.5, for the oxide, we used ρ 0 o = 2.21 g cm 3 . Including volume expansion, in the simulation, the actual density in the oxide layer still is 25% higher than ρ 0 o , since the swelling is mechanically restricted to some extent. In the same figure, we can observe that the plot range of the deformed version is longer which is due to the expansion. Again note, that we do not have to model the full size of the specimen, as the influence of the bulk material are included in the 1D simplification. Furthermore, the mass gain in the case without volume expansion is higher, which we expect, since without swelling, the distances for the diffusing gas are shorter and the concentrations are higher. Fitting the model parameters to the experimental values is also possible when neglecting the volume expansion. However, obviously, the parameters then are different.

Parameter study
In this section, we vary the fitted parameters shown in Table 3. We always vary one single parameter while keeping the rest unchanged.
In the first row in Fig. 9, we illustrate the influence of the reaction rate parameter k. In the plot on the right-hand side, we can observe that a higher reaction rate parameter causes a shorter reaction zone with a steeper change of η. In the plot on the left-hand side, the corresponding mass gain curves are illustrated. In the beginning of the process, the reaction zones are formed and as they are fully developed, the oxide layer grows and the reaction zones move deeper and deeper into the material, keeping their shape. The simulation with the blue curves completes the development of the reaction zone quite early, while the simulation with the green curves finishes late. The initial slope of the mass gain strongly depends on the reaction rate. In the later Fig. 9 Results of a parameter study for the reaction rate parameter (first line), the diffusion coefficient for all phases (second line), the diffusion coefficient for the oxide phase (third line) and the boundary gas concentration (bottom line). Shown are for each the mass gain over time (left) and a profile plot (right) stage, i.e., when the fully developed reaction zone is being shifted, the mass gain curves stay close to each other, as then the process is mostly controlled by the diffusion which behaves the same for all three simulations.
In the second row, we illustrate the influence of the diffusion coefficient, having the same value for the oxide and the starting material. Higher diffusion coefficients cause higher oxygen provision and thus a higher mass gain over time and a faster movement of the reaction zone. Furthermore, higher diffusion coefficients, keeping the reaction rate constant, lead to wider reaction zones. This is similar to the previous case, using smaller reaction zones while keeping the diffusion coefficient constant.
The third investigation addresses the reduction of the diffusion coefficient for the oxide, i.e., formulating D as a function of η. For this, we use the third row in Fig. 9. The diffusion coefficient of the oxide controls the shift of the reaction zone after it is developed. Setting it to zero prevents such a shift completely. The process stagnates. Furthermore, the shape of the reaction zones is slightly different, since especially in the reaction zones, we have 0 < η < 1 which already have a reduced diffusion coefficient.
The last investigation of this parameter study focuses on the boundary gas concentration. The behavior is very similar to the second one, i.e., the one which varies the diffusion coefficient of all phases. However, here, in contrast, the slope of η in the reaction zone is unaffected. Due to this, also the mass gain is marginally different, when compared to the second investigation.

Summary and outlook
In this paper, we derived a model for the oxidation process of a PDC which is based on PHPS with silicon filler particles. The model was created within a large deformation framework. One important quantity is the oxidation state variable which determines the solid material composition at each point in the domain. The stresses, the reaction rate and the diffusion depend on free energy functions. We introduced several couplings. For instance, the oxidation induces large reaction swelling which increases the distance that the gas particles have to travel through the oxide layer. In the model, there is the possibility to include the reaction swelling into the reaction rate, i.e., slowing down the reaction rate when lack of space would lead to a higher elastic energy when reacting. However, later we neglected this influence, as for a realistic description, further mechanical modeling, i.e., considering inelastic effects, is needed. Another coupling is the diffusion where a diffusion coefficient is derived which depends on the oxidation state. We introduced the model for 3D simulations and furthermore reduced the model to 1D. We showed that the 1D representation is sufficient for specific cases, i.e., when the growing oxide layer is thin. We also investigated the influence of the reaction swelling on the diffusion kinematics in the oxidation process. We found significant differences for simulation results excluding the reaction swelling and for the ones including the reaction swelling. We performed a curve fit for some of our model parameters; however, the resulting parameter combination is not unique, as other combinations can produce fits of a similar quality. In the last part of the paper, we investigated the behavior of the system within a parameter study. Varying the reaction rate and the diffusion coefficients changed the mass gain curves and the size of the reaction zone in between the oxide phase and the starting material phase. We furthermore were able to create plateaus in the mass gain curve which means a stagnation of the process after a specific time. This was achieved by using a zero diffusion coefficient for the oxide phase. Here, after the reaction zone is developed, no further oxidation proceeds. Behavior like this can be found in other experiments of our group, using boron particles as additional fillers in the PHPS matrix, cf. Smokovych et al. [7]. The motivation for our investigations is to understand the oxidation protection behavior as a coating material for Mo-Si-B alloys. With this paper, we created a basic model which we plan to improve in the future. This includes the consideration of the microscale, i.e., the distribution and the orientation of the filler particles besides gaps and cracks.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Derivation of the stress tensor and the constitutive tensor
The elastic energy density is With the product rule and chain rule, we can write the second Piola-Kirchhoff stress tensor as We use the identities with the fourth-order identity tensor I = ∂C ∂C which is based on the Kronecker delta in index notation I i jkl = δ ik δ jl . Then, the stresses are S = G J − 2 The stiffness tetrad is derived by with E which we give in index notation as E i jkl = −C −1 ik C −1 jl . The stiffness tetrad in current configuration is For clarity, here, we used capital letters in the index for reference placement and lowercase ones for current placement. The spatial tetrad then is

B Residual and tangent for the Newton-Raphson scheme
The numerical implementation within the finite element method is based on the elemental residual and the elemental tangent matrix. The elemental domain is B e . In the following, N I and N J are shape functions for the nodes I or J , respectively, while, e.g., N I ,k represents a spatial derivative with respect to the kth room direction. The residual is [0] T N I N J 0 Note that here, the stress matrix is a different object than the stress vector from Eq. (70), [σ ] = [σ ]. However, it uses the same entries. Furthermore, the tangent moduli matrix consists of components of the tetrad in Eq. (36) with respect to our chosen orthonormal basis. The stress matrix and the tangent moduli matrix are