1. Chemical Gibbs function variation

For the fully coupled chemomechanical process, the following chemical Gibbs function variation is obtained: ^{Reference Hu and Shen68}

where g _c is the chemical Gibbs function, and D _gc is the complement dissipative energy. ${Q^ * }$ and ${W^ * }$ are the complement heart absorbed from the external heart source and work of the external mechanical force. In the case of ZrB₂–SiC oxidation at high temperature, the oxidation reaction occurs at the interface between the oxide and substrate. By the same derivational process in Ref. Reference Hu and Shen68, one gets

where ${\sigma _{ij}}$ , $\mathop {{T_i}}\nolimits^ *$ , u _i , f _i , and ρ are the stress, surface traction, displacement, body force, and density, respectively. ${v_k}$ is the particle velocity. ${c_N}$ , $\mathop {{J_i}}\nolimits^N$ , and ${\mu _N}$ are the molar fraction, diffusion flux, and chemical potential of species N, respectively. ${v_{Nr}}$ and ${\dot w_r}$ are the stoichiometric coefficient and reaction rate of the rth reaction, respectively. $\mathop {{J_N}}\nolimits^ *$ is the diffusion flux of species N on the boundary. V _k is the velocity of the moving interface Γ. n _i is the component of the normal vector of the boundary and interface. $\left\| \cdot \right\|$ denotes the jump.

Due to the arbitrariness of $\delta {u_i}$ and $\delta {\mu _N}$ , one can get from Eq. (A2) the governing equations

and the boundary conditions

On the moving interface, one can get the jump conditions

2. One-dimensional model

At the temperature range of 1200–1800 °C, the chemical reaction can be schematically expressed as

According to the law of mass action, the chemical reaction rate of Eq. (A7) is ^{Reference Krishnamurthy and Srolovitz52,Reference Zhou, Qu and Cherkaoui54}

where ${k_1}$ is the chemical reaction rate constant of Eq. (A7). $C_{{{\rm{O}}_2}}^{i2}$ is the molar concentration of oxygen at the interface 2 in Fig. 1. The molar fraction and molar concentration are related by ${c_{{{\rm{O}}_2}}} = {V_{\rm{m}}}{C_{{{\rm{O}}_2}}}$ , where V _m is the average molar volume of the oxide per atom.

Based on the assumption that the amount of oxygen consumed by the SiC is proportional to the volume fraction of SiC, the following expression is obtained.

where ${\dot w_2}$ is the chemical reaction rate of Eq. (A8).

The generated ZrO₂ causes the growth of the oxide layer. According to the conservation of zirconium atom, one can get

where ${h_{{\rm{ox}}}}$ and ${V_{{\rm{Zr}}{{\rm{O}}_2}}}$ are the thickness of the oxide layer and molar volume of the ZrO₂, respectively.

The net increase of liquid products (SiO₂ and B₂O₃) first fills up the pore in the oxide layer, then forms the outer glassy layer. The outer glassy layer thickness h _g is expressed as

where ${J_{{{\rm{B}}_2}{{\rm{O}}_3} - {\rm{vap}}}}$ , ${J_{{\rm{Si}}{{\rm{O}}_2} - {\rm{vap}}}}$ , and ${J_{{\rm{SiO}} - {\rm{vap}}}}$ are the evaporation flux of B₂O₃, SiO₂, and SiO at the surface, respectively. ${V_{{{\rm{B}}_2}{{\rm{O}}_3}}}$ and ${V_{{\rm{Si}}{{\rm{O}}_2}}}$ are the molar volume of B₂O₃ and SiO₂, respectively.

The recession rate of ZrB₂ and SiC is related to the change rate of h _ox through molar volume and volume fraction.

where ${V_{{\rm{Zr}}{{\rm{B}}_2}}}$ and ${V_{{\rm{SiC}}}}$ are the molar volume of ZrB₂ and SiC, respectively. ${R_{{\rm{Zr}}{{\rm{B}}_2}}}$ and R _SiC are the recession thickness of ZrB₂ and SiC, respectively.

Then the SiC-depleted layer thickness h _SiC is easily expressed as

The weight gain due to the oxidation is expressed as

where ${\rho _{{\rm{Zr}}{{\rm{O}}_2}}}$ , ${\rho _{{\rm{SiC}}}}$ , ${\rho _{{\rm{Zr}}{{\rm{B}}_2}}}$ , and ${\rho _{\rm{g}}}$ are the density of ZrO₂, SiC, ZrB₂, and borosilicate, respectively.

Ignoring the convection item, Eq. (A4) is written as

where $J_3^{{{\rm{O}}_2}}$ is the flux of oxygen in x ₃ direction.

The stress-dependent chemical potential of oxygen is ^{Reference Li69}

where R and T are the gas constant and absolute temperature, respectively. ${\eta _{{{\rm{O}}_2}}}$ and ${\sigma _{kk}}$ are the coefficient of compositional expansion and hydrostatic stress in the oxide layer. $\mu _{{{\rm{O}}_2}}^0$ is the chemical potential at the reference state.

The oxygen flux is expressed as

where ${D_{{{\rm{O}}_2} - {\rm{g}}}}$ is the diffusivity of oxygen in the glassy borosilicate.

Then, Eq. (A17) can be rewritten as

Since the chemical potential is stress-dependent, the chemical reaction is stress-dependent too. The stress-dependent chemical reaction rate constant is expressed as ^{Reference Zhou, Qu and Cherkaoui54,Reference Swaminathan, Qu and Sun70}

where p represents the oxide. $k_1^0$ is the chemical reaction rate constant when there is no stress.

Stress-dependent chemical reaction rate is

Substituting Eqs. (A22) and (A23) into Eqs. (A11) and (A12), they are rewritten as

The evaporation flux at the surface is express as follows

where ${D_{{\rm{species}}}}$ is the diffusivity in mixed gas. ${p_{{\rm{species - vap}}}}$ is the evaporation pressure. Species refers to SiO₂, B₂O₃, and SiO. ${\delta _{{\rm{bdry}}}}$ is the boundary layer thickness. ${l_{{\rm{specimen}}}}$ is the length of the specimen. ${v_{{\rm{fluid}}}}$ , ${\eta _{{\rm{fluid}}}}$ , and ${\rho _{{\rm{fluid}}}}$ are the velocity, viscosity, and density of the ambient fluid, respectively.

According to Eq. (A6), the moving boundary condition of the interface 2 is

where $J_3^{{{\rm{O}}_2},i2}$ is the oxygen flux at the interface 2. $V$ is the velocity of the interface 2 and is expressed as

At the interface 3, the flux of the oxygen is equal.

Since the quasi-steady state is assumed, the oxygen flux can be expressed as

where ${D_{23}}$ and ${D_{3a}}$ are the oxygen diffusivity in the region (2–3) and region (3–a), respectively. Since B₂O₃ evaporates much faster than SiO₂ at the surface, the glassy layer contains more SiO₂ than borosilicate filled in the porous channels. As a result, the diffusivity for oxygen in region (2–3) is different from that in region (3–a). $C_{{{\rm{O}}_2}}^{i3}$ and $C_{{{\rm{O}}_2}}^a$ are the oxygen concentration at the interface 3 and a, respectively. $C_{{{\rm{O}}_2}}^a$ is expressed as

where $p_{{{\rm{O}}_2}}^a$ and $S$ are the oxygen partial pressure at the air interface and solubility of oxygen in the glassy borosilicate.

Substituting Eqs. (A30) and (A31) into Eq. (A29), we obtain

and then we get the oxygen flux at the moving boundary

Solving Eqs. (A27), (A33), and (A34) simultaneously, we obtain the oxygen concentration at the interface 2, $C_{{{\rm{O}}_2}}^{i2}$ .

Accompanying the oxide layer growth, the stress is produced. Based on Clarke’s assumption, ^{Reference Clarke44,Reference Suo and Shen71} the lateral growth strain is proportional to the oxide layer thickness.

where ${D_{{\rm{ox}}}}$ is the lateral growth coefficient.

The elastic strain is

where i can be replaced by ox for oxide and sub for substrate. ${\nu _i}$ , ${E_i}$ , and ${\sigma _i}$ are the Poisson ratio, Young’s modulus and plane stress, respectively.

For the porous oxide layer, the relationship between Yong’s modulus and porosity is ^{Reference Li, Li, Zhang and Fang61}

where $E_{{\rm{ox}}}^0$ is the Yong’s modulus of oxide when there is no pore.

For the high-temperature oxidation, the creep effect should be taken into account. Here, Norton’s creep power law ^{Reference Golan, Arbel, Eliezer and Moreno72} is adopted in the work.

where B, Q, and n are the material constant, activation energy for creep process and the stress exponent, respectively. Since the stress is compressive in the oxide layer and tensile in the substrate, ${\rm{sign}}({\sigma _{{\rm{ox}}}}) = - 1$ and ${\rm{sign}}({\sigma _{{\rm{sub}}}}) = 1$ .

The total strain in the oxide and substrate is expressed, respectively, as

For the cracking, bulking and spalling are not taken into account in the work, the strain in the oxide and substrate should be compatible.

The rate of Eq. (A41) is

Substituting Eqs. (A35), (A36) and (A38) into Eq. (A42), we obtain

The force equilibrium equation is

where ${h_{{\rm{sub}}}}$ is the thickness of substrate. In the work, we deal with the average stress in the oxide and substrate. From Eq. (A44), ${\sigma _{{\rm{sub}}}}$ can be expressed as

The rate of Eq. (A45) is

Substituting Eqs. (A45) and (A46) into Eq. (A43), we obtain

Rearranging Eq. (A47), we get

In the same way, we get

The equations can be calculated by the evolutionary algorithm by assuming an initial thickness of ${h_{{\rm{ox}}}}$ and ${h_{\rm{g}}}$ .

The variation of pore radius r and borosilicate volume fraction f _g is

where ${r^0}$ and $f_{\rm{g}}^0$ are the pore radius and glassy borosilicate volume fraction when there is no stress.

The relationship between the macroscopic strain and the microscopic strain in the pore is expressed as ^{Reference Dormieux, Kondo and Ulm73}

where I is the forth-order identity tensor. S is the Eshelby tensor. E is the macroscopic strain. If the effect of lateral growth strain on the pore strain is neglected, the pore strain is

To implement the model, the equations can be rewritten in incremental forms with respect to incremental time, $\Delta t$ . Eqs. (A24) and (A25) are rewritten as

The thickness at ${t^{(k + 1)}}$ can be easily obtained as

In the same way, Eqs. (A48) and (A49) are rewritten as

The stress at ${t^{(k + 1)}}$ can be easily obtained as