« PreviousUpNext »Contents
Previous: 4.2 Implementation in FEDOS    Top: 4 Numerical Implementation    Next: 5 Results and Discussion

4.3  Implementation in COMSOL

COMSOL [30] is a commercial program widely used in industry and research for the simulation of FEM based problems. It offers a complete solution for various physics and engineering applications starting from the geometry generation over the modeling to the simulation and evaluation with various plotting tools and data exporting capabilities. In addition it offers a material data base from where material parameters can be imported.

4.3.1  Electro-Mechanical Model

The electro-mechanical model is a preconfigured model including the electric potential

(4.30–4.32) \{begin}{gather} \mathbf {E}=-\nabla V \mathref {(4.30)}\\ \displaystyle \mathbf {J}=\Bigl (\sigma +\epsilon _{0}\epsilon _{\mathrm {r}}\frac {\partial }{\partial t}\Bigr )\mathbf {E}+\mathbf {J}_{\mathrm {e}} \mathref {(4.31)} \\ \nabla \cdot \mathbf {J}=Q_{\mathrm {e}}, \mathref {(4.32)} \{end}{gather}

where \( V \) stands for the electric potential, \( \mathbf {J} \) for the current density, \( \mathbf {J}_{\mathrm {e}} \) for a user predefined current density, \( Q_{\mathrm {e}} \) for the free charge distribution and \( \sigma \) for the conductivity. The Joule heating modeling is formulated by the equation system

(4.33) \{begin}{gather} Q=\mathbf {E}\cdot \mathbf {J} \mathref {(4.33)} \{end}{gather}

(4.34) \{begin}{gather} \displaystyle \rho C_{\mathrm {p}}\frac {\partial T}{\partial t}+\rho C_{\mathrm {p}}\mathbf {u}\cdot \nabla T=\nabla \cdot (k\nabla T)+Q. \mathref {(4.34)} \{end}{gather}

where \( Q \) stand for the thermal energy dissipated per time, \( C_{\mathrm {p}} \) for the thermal expansion coefficient under constant pressure, \( \rho \) for the volumetric mass density, \( k \) for the thermal conductivity, and \( \mathbf {u} \) for the velocity field. The mechanical stress modeling is implemented

Figure 4.5.:  Flowchart of the EM simulation in FEDOS.

for small deformations by

(4.35–4.38) \{begin}{align} -\nabla \cdot \pmb \sigma &=\b {F}_{\mathrm {v}} \\ \pmb \sigma -\pmb \sigma _{0}&=C: (\pmb \epsilon -\pmb \epsilon _{0}-\pmb \epsilon _{\mathrm {i}\mathrm {n}\mathrm {e}1}) \\ \pmb \epsilon _{\mathrm {i}\mathrm {n}\mathrm {el}}&=\alpha (T-T_{\mathrm {r}\mathrm {e}\mathrm {f}})\b {1} \\ \pmb \epsilon &=\frac {1}{2}\Bigl ((\nabla \mathbf {u})^{T}+\nabla \mathbf {u}\Bigr ). \{end}{align}

where \( \b {F}_{\mathrm {V}} \) stands for the body force density, \( T_{\mathrm {r}\mathrm {e}\mathrm {f}} \) for the reference temperature of the thermal expansion, and \( \pmb \epsilon _{\mathrm {i}\mathrm {n}\mathrm {e}1} \) for the strain originating from the thermal expansion. This model has the ability to include strains by \( \pmb \epsilon _{0} \) or stresses by \( \pmb \sigma _{0} \) from another model and is used to includes the strain arising from the vacancy flux and vacancy generation/annihilation in the EM model.

4.3.2  PDE as Electromigration Model

Since not for all physical problems models are included in COMSOL, it provides the “Coefficient Form PDE”, which allows to express a huge set of technical relevant linear and non-linear PDEs. This general PDE has the following form:

(4.39) \{begin}{gather} e_{a}\displaystyle \frac {\partial ^{2}u}{\partial t^{2}}+d_{a}\frac {\partial u}{\partial t}-\nabla \cdot (c\nabla u+\pmb \alpha u-\pmb \gamma )+\pmb \beta \cdot \nabla u+au=f \mathref {(4.39)} \{end}{gather}

EM is represented by two differential equations, one for the vacancy flux and one for the strain due to the movement of the vacancies as well as their generation/annihilation.
The vacancy dynamics is modeled by a conservation law with a generation/annihilation term given by (3.9). By comparing this equation against (4.39) it follows that \( e_{a}, \pmb \gamma \) and \( \pmb \beta \) have to be set to zero and \( d_{a} \) has to be set to one. Taking this into account and rearranging the equation leads to

(4.40) \{begin}{gather} \displaystyle \frac {\partial u}{\partial t}=\nabla \cdot (c\nabla u+\pmb \alpha u)+(f-au). \mathref {(4.40)} \{end}{gather}

The equation parameters have to be set as follows to model vacancy dynamics.

(4.41–4.44) \{begin}{align} c&=-D_{\mathrm {v}} \\ \displaystyle \pmb \alpha &=D_{\mathrm {v}}\Bigl (\frac {|Z^{*}|e}{k_{\mathrm {B}}T}\mathbf {E}-\frac {f\Omega }{k_{\mathrm {B}}T}\nabla \sigma +\frac {Q^{*}}{k_{\mathrm {B}}T^{2}}\nabla T\Bigr ) \\ f&=\displaystyle \ \frac {C_{\mathrm {e}\mathrm {q}}}{\tau _{\mathrm {v},\mathrm {e}\mathrm {q}}} \\ a&=\frac {1}{\tau _{\mathrm {e}\mathrm {q}}} \{end}{align}

The electric field \( \mathbf {E} \) and the hydraulic pressure \( \sigma \) are taken from the electro-mechanical model (cf. Section 4.3.1).

For the strain build-up by the vacancy dynamics and generation/annihilation the same general equation (4.39) can be used with the parameters set as following:

(4.45–4.48) \{begin}{align} c&=-\Omega (1-f) \\ \displaystyle \pmb \alpha &=\Omega (1-f)\Bigl (\frac {|Z^{*}|e}{k_{\mathrm {B}}T}\mathbf {E}-\frac {f\Omega }{k_{\mathrm {B}}T}\nabla \sigma +\frac {Q^{*}}{k_{\mathrm {B}}T^{2}}\nabla T\Bigr ) \\ f&=\displaystyle \Omega f\frac {C_{\mathrm {e}\mathrm {q}}}{\tau _{\mathrm {v},\mathrm {e}\mathrm {q}}} \\ a&=\displaystyle \Omega f\frac {1}{\tau _{\mathrm {e}\mathrm {q}}} \{end}{align}

Since also the strain build-up is described by a first derivative in time analogously the parameters \( e_{a}, \pmb \gamma \), and \( \pmb \beta \) are set to zero and \( d_{a} \) to one. It should be emphasized that the vacancy dynamics model is connected through the stress gradient to the solid mechanics model, which by itself depends on the strain build-up by the vacancies as an input value making those models highly nonlinear.

4.3.3  Phase Field Model

In COMSOL a phase field model is implemented by the equation

(4.49) \{begin}{gather} \displaystyle \frac {\partial \phi }{\partial t}+\mathbf {u}\cdot \nabla \phi =\nabla \ \frac {\gamma \lambda }{\epsilon _{\mathrm {p}\mathrm {f}}^{2}}\nabla \psi (\phi ), \mathref {(4.49)} \{end}{gather}

where \( \mathbf {u} \) is a velocity field. Taking into account that the velocity field is divergence free \( (\nabla \cdot \mathbf {u}=0) \) this equation can be rearranged to

(4.50) \{begin}{gather} \displaystyle \frac {\partial \phi }{\partial t}=\nabla \cdot \left (-\mathbf {u}\phi +\frac {\gamma \lambda }{\epsilon _{\mathrm {p}\mathrm {f}}^{2}}\nabla \psi (\psi )\right ). \mathref {(4.50)} \{end}{gather}

Independent of the form of the function \( \psi (\phi ) \) , this equation has the form of a conservation.

(4.51) \{begin}{gather} \displaystyle \frac {\partial \phi }{\partial t}=\nabla \cdot \mathbf {J} \mathref {(4.51)} \{end}{gather}

Therefore, the phase field model of COMSOL is not suitable to implement the vacancy exchange between the bulk and the interface of the void evolution model, described in Section 3.6, represented by the last term of

(4.52) \{begin}{gather} \displaystyle \frac {\partial \phi }{\partial t}=\frac {2}{\epsilon _{pf}\pi }\nabla \cdot \mathcal {D}_{s}(\phi )(\nabla \mu _{s}-eZ\mathbf {E})-\frac {4A}{\epsilon _{pf}\pi }(\mu _{s}-\mu _{v}). \mathref {(4.52)} \{end}{gather}

Therefore, for the implementation of the phase field model the “Coefficient Form PDE” (4.39) has to be chosen with the parameters set as follows:

(4.53–4.59) \{begin}{align} e_{a}&=0, \\ \displaystyle \pmb \gamma &=\frac {2}{\epsilon _{pf}\pi }\mathcal {D}_{s}(\phi )(\nabla \mu _{s}-eZ\mathbf {E}), \\ \pmb \beta &=\b {0}, \\ c&=0, \\ \pmb \alpha &=\b {0}, \\ f&=-\displaystyle \frac {4A}{\epsilon _{pf}\pi }(\mu _{s}-\mu _{v}), \\ a&=0. \mathref {(4.59)} \{end}{align}

Figure 4.6.:  Flowchart of the EM simulation in COMSOL.

Figure 4.7.:  Flowchart of the EM induced void simulation in COMSOL.

4.3.4  Time Stepping

COMSOL offers two methods to compute the solution. The first method called ”Fully Coupled” generates the matrices \( \b {A} \) and \( \mathbf {C} \) of (4.27) where all differential equations are included at once. These matrices depend on the solution vector (e.g. the coupling of the vacancy flux and the stress equations) and therefore, first (4.27) is solved where \( \pmb \beta \) and \( \pmb \gamma \) is retrieved from the old solution vectors. This step is followed by a regeneration of the matricies using \( \pmb \beta \) and \( \pmb \gamma \) calculated from the new solution vector and solving again.
These steps are repeated untill convergence is reached. The simulation is terminated under user chosen conditions (e.g. maximum step number), if the problem does not converge. After the solution has converged, the next time step is computed. The time step size depends on various conditions (e.g. step size to the next chosen output time).

The second method called “Segregated Step” generates the matrices \( \b {A} \) and \( \mathbf {C} \) in (4.27) for every system separately and calculates the solutions sequentially. The number and the included differential equations of the systems can by chosen by the user. For EM without void modeling the sequence is depicted in Figure 4.6 and for EM modeling with void modeling in Figure 4.7. This method is the method of choice as the Fully Coupled” method often needs too much memory, since the size of the matrices grows with second order of the number of mesh points.

« PreviousUpNext »Contents
Previous: 4.2 Implementation in FEDOS    Top: 4 Numerical Implementation    Next: 5 Results and Discussion