3.1.2 Visco-Elastic model using FEM

A visco-elastic oxide growth model has been developed at our institute using FEM as the simulation environment [80]. The model is three-dimensional and does not rely on a simple expression for the interface motion, but rather relies on calculations of oxidant diffusion and oxide volume expansion in three dimensions to move the FEM mesh appropriately. A major limitation of FEM analysis is related to the mesh used for simulations. As either the time or space increments are decreased, or the mesh is made finer, the number of calculations necessary to simulate oxide growth can increase by a factor dependent on the model itself [34], [159]. Some models also do not account for orientation dependence adequately or the error can be a non-linear function which can amplify significantly with further iterations [229].

The model presented in [80] is based on a few main equations. The first equation is meant to describe oxidant diffusion

$$D\left(T\right)\Delta C\left(\vec{x},t\right)=k\left(\eta\right)C\left(\vec{x},t\right),$$ (72)

where $\Delta=\cfrac{\partial^{2}}{\partial x^{2}}\,+\,\cfrac{\partial^{2}}{\partial y^{2}}\,+\,\cfrac{\partial^{2}}{\partial z^{2}}$ is the Laplace operator, $C\left(\vec{x},t\right)$ is the oxidant concentration in the material, and $D\left(T\right)$ is the temperature dependent low stress diffusion coefficient. $k\left(\eta\right)$ is the strength of a spatial sink and not a simple reaction coefficient at a sharp interface

$$k\left(\eta\right)=\eta\left(\vec{x},t\right)\cdot k_{max},$$ (73)

where $k_{max}$ is the maximal possible strength of the sink.

The next equation is meant to deal with the dynamics of $\eta$

$$\cfrac{\partial\eta\left(\vec{x,t}\right)}{\partial t}=-\cfrac{1}{\lambda}\cdot\cfrac{k\left(\eta\right)C\left(\vec{x},t\right)}{N},$$ (74)

where $\lambda$ is the volume expansion factor (=2.25) for the reaction from Si to SiO$_2$, and $N$ is the number of oxidant molecules incorporated into a unit of SiO$_2$ volume.

The volume increase of the generating oxide occurs successively and not abruptly because of the diffuse interface concept. A volume increase of the oxidized material is calculated using the $\eta$ and $C$ values. After a given time $\Delta t$, the normalized additional volume is given by

$$V_{rel}^{add}=\cfrac{\lambda-1}{\lambda}\cdot\cfrac{\Delta t\,k\left(\eta\right)C\left(\vec{x},t\right)}{N}.$$ (75)

An important aspect of (3.4) is that the sum of $V^{add}_{rel}$ over all time steps must not exceed 125%, which is the maximum volume expansion of the material during oxidation.

Since the principal axis components of the residual strain tensor $\tilde{\epsilon}_{0}$ are linearly proportional to $V^{add}_{rel}$ in the form

$$\epsilon_{0,xx}=\epsilon_{0,yy}=\epsilon_{0,zz}=\cfrac{1}{3}\,V_{rel}^{add},$$ (76)

the normalized additional volume directly loads the mechanical problem. Using Hook's law, the stress tensor $\tilde{\sigma}$ is given by

$$\tilde{\sigma}=\mathbf{D}\left(\tilde{\epsilon}-\tilde{\epsilon_{0}}\right)+\tilde{\sigma_{0}},$$ (77)

where $\mathbf{D}$ is the so-called material matrix, constructed in Lame's form [80], $\tilde{\epsilon}$ is the strain tensor, $\tilde{\epsilon}_{0}$ is the residual strain tensor, and $\tilde{\sigma}_{0}$ is the residual stress tensor. The material matrix is rebuilt and (3.6) becomes

$${\scriptstyle \left(\begin{array}{c} \sigma_{xx}\\ \sigma_{yy}\\ ... ...y}\\ \gamma_{yz}-\gamma_{0,yz}\\ \gamma_{zx}-\gamma_{0,zx} \end{array}\right)},$$ (78)

where $E$ is the Young's modulus, $\nu$ is the Poisson ratio, and $\gamma$ are the shear strain components.

The effective shear modulus $G_{eff}$ can be added to $\mathbf{D}$ to handle elastic and visco-elastic materials. For elastic materials

$$G_{eff}=G=\cfrac{E}{2\left(1+\nu\right)},$$ (79)

while for visco-elastic materials

$$G_{eff}=G\cfrac{\tau}{\Delta T}\left(1-e^{\left(-\Delta T\slash\tau\right)}\right),$$ (80)

where $\tau$ is the Maxwellian relaxation time constant and $G$ is the effective shear modulus.