Subsections

• 3.8.1 Discretization of the Simple Diffusion Model
• 3.8.2 Discretization of the Simple Extrinsic Diffusion Model
• 3.8.3 Discretization of the Three-Stream Mulvaney-Richardson Model
• 3.8.4 Analytical Solution of the Segregation Problem for One Dimension
• 3.8.5 Numerical Handling of the Segregation Model

3.8 Numerical Handling of the Diffusion Models

In following section we will present discretization approaches for several diffusion models described in the preceding sections. Starting with the partial differential equations which model the given diffusion phenomena, we construct the nucleus matrix and residuum vector which are needed to carry out the procedure presented in Section 2.8. Throughout this section the number of nodes of discretization $T_h(\Omega)$ is denoted as $N$.

3.8.1 Discretization of the Simple Diffusion Model

The discretization of the simple diffusion models follows directly from the discussion presented in Section 2.3. In order to construct the global matrix of the system, we write nucleus matrix of the simple diffusion model for each $T\in T_h(\Omega)$, $rank(\mathbf{\Pi}(T))= 4$,

$$\mathbf{\Pi}(T) = \mathbf{K}(T) + D \Delta t \mathbf{M}(T).$$ (3.76)

$\Delta t$ is the time step of the discretisized time and $\mathbf{K}(T)$ and $\mathbf{M}(T)$ the stiffness and mass matrices defined on single tetrahedra $T$ from $T_h(\Omega)$,

$$K_{pq}(T)=\int\limits_T \nabla \varphi_p(\mathbf{x}) \nabla \varp... ... M_{pq}(T)=\int\limits_T \varphi_p(\mathbf{x}) \varphi_q(\mathbf{x}) d\Omega,$$

for $p,q\in\{1,2,3,4\}$. This problem is linear and the calculation of the nucleus matrix $\mathbf{\Pi}(T)$ as the Jacobi matrix of the operator $\ell^e : \Bbb{R}^{4}\rightarrow \Bbb{R}^{4}$ (introduced previously in Section 2.5) is trivial.
the global matrix $\mathbf{G}$ is assembled out of the $\mathbf{\Pi}(T)$ matrices and has obviosly a dimension equal to the number of all points of discretization.

3.8.2 Discretization of the Simple Extrinsic Diffusion Model

We take the functions $w(\mathbf{x})$, $v(\mathbf{x})$ and $u(\mathbf{x})$ defined on bounded open domen $\Omega$. Starting from Green's theorem, a simple relationship can be proved,

$$\int\limits_\Omega \nabla v \nabla w u d\Omega + \int\limits_\O... ...w v d\Omega=\int\limits_\Gamma v u \frac{\partial w}{\partial n} d\Gamma,$$ (3.77)

where $\Gamma$ is the boundary of the domain $\Omega$. We write (3.33) as

$$\nabla (h(C_A)D) \nabla C_A + D h(C_A) \triangle C_A = \frac{\partial C_A}{\partial t}.$$ (3.78)

By multiplying (3.78) with the basis function $\varphi_i = \varphi_i(\mathbf{x})$ and integrating over the domen $\Omega$ we have,

$$\int\limits_\Omega \nabla (h(C_A)D) \nabla C_A \varphi_i d\O... ...mega = \int\limits_\Omega \frac{\partial C_A}{\partial t} \varphi_i d\Omega.$$ (3.79)

Applying (3.77) on (3.79) gives,

$$D\int\limits_\Gamma h(C_A) \varphi_i \frac{\partial C_A}{\part... ...Omega = \int\limits_\Omega \frac{\partial C_A}{\partial t} \varphi_i d\Omega.$$ (3.80)

Assuming zero-Neumann boundary conditions on $\Gamma$ we obtain the weak formulation of the equation (3.33)

$$\int\limits_\Omega \frac{\partial C_A}{\partial t} \varphi_i d\Omega+D\int\limits_\Omega h(C_A) \nabla C_A \nabla \varphi_i d\Omega = 0.$$ (3.81)

By introducing time discretization with time step $\Delta t$ and writing the last equation for the single element $T\in T_h(\Omega)$ we obtain,

$$\int\limits_T (C^n_A- C^{n-1}_A) \varphi_i d\Omega+ D\Delta t \int\limits_T h(C_A^n) \nabla C_A^n \nabla \varphi_i d\Omega = 0.$$ (3.82)

Functions $C_A(\mathbf{x},t_n)$, $\nabla C_A(\mathbf{x},t_n)$, and $h(C_A(\mathbf{x},t_n))$, for $t_n=n\cdot \Delta t$, are linearly approximated on the element $T\in T_h(\Omega)$,

$$C_A^n=C_A(\mathbf{x},t_n) = \sum_{i=1}^{4} \varphi_i C^{n}_{A,i},$$ (3.83)

$$\nabla C_A^n = \nabla C_A(\mathbf{x},t_n) =\sum_{i=1}^{4} \nabla \varphi_i C^{n}_{A,i},$$ (3.84)

$$h(C_A^n)= h(C_A(\mathbf{x},t_n))= h(C^{n}_{A,m}).$$ (3.85)

The $C_1^n, C_2^n, C_3^n, C_4^n$, are the values of the concentration at the nodes of the element $T$, $C^{n}_{m}$ is the concentration value at some point inside the $T$. Normally, for $C^{n}_{A,m}$ we use the following simple approximation,

$$C^{n}_{A,m} = \frac{1}{4} \sum_{i=1}^{4} C^{n}_{A,i}.$$ (3.86)

We can define discrete operator $\ell^e$ for each $T\in T_h(\Omega)$. Since we have to deal only with a single partial differential equation and not with a system of equations, we can omit the second index in (2.27) by defining the operator,

$$\ell_{i}^e(C_{A,1}^n,C_{A,2}^n,C_{A,3}^n,C_{A,4}^n)=\frac{1}{\Del... ...{ik} C_{A,i}^n+D h(C_{A,m}^n)\sum_{i=1}^4 K_{ik} C_{A,i}^n,\quad k=1,2,3,4.$$ (3.87)

The nucleus matrix is than defined as

$$\mathbf{\Pi}(T)=\frac{\partial (\ell_{1}^e,\ell_{2}^e,\ell_{3}^e,\ell_{4}^e)}{\partial (C_{A,1}^n,C_{A,2}^n,C_{A,3}^n,C_{A,4}^n)}.$$ (3.88)

and the residuum vector as,

$$R_k^e=\ell_{k}^e(C_{A,1}^n,C_{A,2}^n,C_{A,3}^n,C_{A,4}^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik} C_{A,i}^{n-1},\quad k=1,2,3,4.$$ (3.89)

The nucleus matrix for this simple case can be also calculated analytically. In this case is also $rank(\mathbf{G})=N$.

3.8.3 Discretization of the Three-Stream Mulvaney-Richardson Model

First we consider dopant continuity equation of Mulvaney-Richardson model (3.45). The local change of dopants is due to the divergence of fluxes of dopant-interstitial and dopand vacancy pairs,

$$\frac{\partial C_A}{\partial t}=\nabla\cdot(\mathbf{J}_{AI}+\mathbf{J}_{AV}).$$ (3.90)

In this model the diffusion coefficient $D_A$ of the point defect-dopant pairs is considered constant and we can write (based on (3.45)),

$$\nabla\cdot\mathbf{J}_{AI}= \frac{f_I D_A}{C_{I}^{eq}}\Delta (C_I C_A) + \frac{f_I D_A}{C_{I}^{eq}}\nabla(C_I C_A \nabla ( n)),$$ (3.91)

$$\nabla\cdot\mathbf{J}_{AV}=\frac{(1-f_I) D_A}{C_{V}^{eq}}\Delta (C_V C_A) + \frac{(1-f_I) D_A}{C_{V}^{eq}}\nabla(C_V C_A \nabla ( n)).$$ (3.92)

In order to obtain a weak formulation of (3.91) on the single element $T\in T_h(\Omega)$ we multiply both sides of the equation with the basis function $\varphi_k$, and integrate over $T$,

$$\int_T\nabla\cdot\mathbf{J}_{AI} \varphi_k d\Omega = -\frac{f_... ...a \varphi_k d\Omega +\frac{f_I D_A}{C_{I}^{eq}}\int_T \nabla(C_I C_A \nabla ( n)) \varphi_k d\Omega=$$ (3.93)

$$=-\frac{f_I D_A}{C_{I}^{eq}}\int_T \nabla(C_I C_A)\nabla \varphi_k d\Omega-\frac{f_I D_A}{C_{I}^{eq}}\int_T C_I C_A \nabla ( n) \nabla \varphi_k d\Omega.$$ (3.94)

The nonlinear terms of (3.94), $C_I(\mathbf{x}) C_A(\mathbf{x})$, $\nabla(C_I(\mathbf{x}) C_A(\mathbf{x}))$, and ln$n$ are expressed by the basis nodal functions $\varphi_1,\varphi_2,\varphi_3,\varphi_4$ on the following way,

$$C_I(\mathbf{x}) C_A(\mathbf{x}) = \sum_{i,j=1}^4 C_{I,i}C_{A,j}\varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}),$$

$$\nabla(C_I(\mathbf{x}) C_A(\mathbf{x})) = \sum_{i,j=1}^4 C_{I,i}C_{A,j}\nabla(\varphi_i(\mathbf{x}) \varphi_j(\mathbf{x})),$$ (3.95)

$$\nabla( n)=\frac{1}{\sqrt{4n_{i}^{2}+C_{A,m}^{2}}} \sum_{i=1}^4 C_{A,i}\nabla\varphi_i(\mathbf{x}).$$

In the last expression we used the relationship (3.13) for the electron concentration $n$ in the case of single charged dopant. We substitute the relationships (3.95) into (3.94) and obtain,

$$\int_T \nabla\cdot\mathbf{J}_{AI} \varphi_k d\Omega=-\frac{f_I ... ...^4C_{I,i}C_{A,j}\int_T \nabla (\varphi_i \varphi_j)\nabla \varphi_k d\Omega -$$

$$-\frac{f_I D_A}{C_{I}^{eq}\sqrt{4n_{i}^{2}+C_{A,m}^{2}}}\sum_{i,j... ...}C_{A,l}\int_T \varphi_i \varphi_j \nabla \varphi_l \nabla \varphi_k d\Omega.$$ (3.96)

To complete the formulation of the weak equation form we have to evaluate following two integrals,

$$\int_T \nabla (\varphi_i \varphi_j)\nabla \varphi_k d\Omega=det... ...}_i \varphi^{e}_j)\mathbf{\Lambda}\nabla_e \varphi^{e}_k d\zeta d\eta d\xi,$$ (3.97)

$$\int_T \varphi_i \varphi_j \nabla \varphi_l \nabla \varphi_k d\... ...a_e \varphi^{e}_l \mathbf{\Lambda}\nabla_e \varphi^{e}_k d\zeta d\eta d\xi,$$ (3.98)

here we transfer the calculation to the normed $(\xi,\zeta,\eta)$ coordinate system. Transformation of the gradient operators is done by multiplication with the matrix $\mathbf{\Lambda}=\mathbf{(J^{-1})^T}$.

$$\mathbf{\Lambda}\nabla_e (\varphi^{e}_i \varphi^{e}_j)\mathbf{\La... ...partial\eta} \frac{\partial(\varphi^{e}_k)}{\partial\zeta} \end{bmatrix}=$$

$$=\Bigl(\Lambda_{0,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{... ...\partial\eta}+\Lambda_{0,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr)+$$

$$+ \Bigl(\Lambda_{1,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}... ...tial\eta}\!+\!\Lambda_{1,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr)+$$

$$+ \Bigl(\Lambda_{2,0}\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}... ...l\eta}+\Lambda_{2,2}\frac{\partial(\varphi^{e}_k)}{\partial\zeta}\Bigr).\quad$$ (3.99)

Furthermore we rearrange (3.99),

$$\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\xi}\frac{\p... ...bda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lambda_{1,1}+\Lambda_{2,0}\Lambda_{2,1}) +$$

$$+\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\xi}\frac{\... ...mbda_{0,1}\Lambda_{0,0}+\Lambda_{1,1}\Lambda_{1,0}+\Lambda_{2,1}\Lambda_{2,0})+$$

$$+ \frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\eta}\frac... ...mbda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\Lambda_{1,2}+\Lambda_{2,1}\Lambda_{2,2})+$$

$$+\frac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\zeta}\frac... ...mbda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lambda_{1,1}+\Lambda_{2,2}\Lambda_{2,1})+$$

$$+\frac{\partial(\varphi^{e}_i\varphi^{e}_j)}{\partial\zeta}\frac{... ...^{e}_k)}{\partial\zeta}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}).$$ (3.100)

and introduce the matrix $\mathbf{K}_{\xi}$, $\mathbf{K}_{\eta}$, $\mathbf{K}_{\zeta}$.

$${K}_{\xi}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zeta=... ...frac{\partial(\varphi^{e}_i \varphi ^{e}_j)}{\partial\xi} d\zeta d\eta d\xi,$$

$${K}_{\eta}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zeta... ...rac{\partial(\varphi^{e}_i \varphi ^{e}_j)}{\partial\eta} d\zeta d\eta d\xi,$$ (3.101)

$${K}_{\zeta}^{i,j}=\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi}\int_{\zet... ...rac{\partial(\varphi^{e}_i \varphi^{e}_j)}{\partial\zeta} d\zeta d\eta d\xi.$$

The values of the integrals (3.97) and (3.98) can be expressed as functions of local node indexes,

$$\alpha_{T}(i,j,k)=\int_T \nabla (\varphi_i \varphi_j)\nabla \varp... ...i,j}( \varphi_{\xi}^{k}(\Lambda_{0,0}^{2}+\Lambda_{1,0}^{2}+\Lambda_{2,0}^{2})+$$

$$+\varphi_{\eta}^{k}(\Lambda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lamb... ...bda_{0,0}\Lambda_{0,2}+\Lambda_{1,0}\Lambda_{1,2}+\Lambda_{2,0}\Lambda_{2,2}))+$$

$$+ {K}_{\eta}^{i,j}(\varphi_{\xi}^{k}(\Lambda_{0,1}\Lambda_{0,0}+\... ...,0})+\varphi_{\eta}^{k}(\Lambda_{0,1}^{2}+\Lambda_{1,1}^{2}+\Lambda_{2,1}^{2})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\La... ...mbda_{0,2}\Lambda_{0,0}+\Lambda_{1,2}\Lambda_{1,0}+\Lambda_{2,2}\Lambda_{2,0})+$$

$$+\varphi_{\eta}^{k}(\Lambda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lamb... ...+ \varphi_{\zeta}^{k}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}))),$$ (3.102)

$$\beta_{T}(i,j,l,k)=\int_T \varphi_i \varphi_j \nabla \varphi_l \n... ...^{l}( \varphi_{\xi}^{k}(\Lambda_{0,0}^{2}+\Lambda_{1,0}^{2}+\Lambda_{2,0}^{2})+$$

$$+\varphi_{\eta}^{k}(\Lambda_{0,0}\Lambda_{0,1}+\Lambda_{1,0}\Lamb... ...bda_{0,0}\Lambda_{0,2}+\Lambda_{1,0}\Lambda_{1,2}+\Lambda_{2,0}\Lambda_{2,2}))+$$

$$+{\varphi}_{\eta}^{i}(\varphi_{\xi}^{k}(\Lambda_{0,1}\Lambda_{0,0... ...,0})+\varphi_{\eta}^{k}(\Lambda_{0,1}^{2}+\Lambda_{1,1}^{2}+\Lambda_{2,1}^{2})+$$

$$+ \varphi_{\zeta}^{k}(\Lambda_{0,1}\Lambda_{0,2}+\Lambda_{1,1}\La... ...mbda_{0,2}\Lambda_{0,0}+\Lambda_{1,2}\Lambda_{1,0}+\Lambda_{2,2}\Lambda_{2,0})+$$

$$+\varphi _{\eta}^{k}(\Lambda_{0,2}\Lambda_{0,1}+\Lambda_{1,2}\Lam... ...)+\varphi_{\zeta}^{k}(\Lambda_{0,2}^{2}+\Lambda_{1,2}^{2}+\Lambda_{2,2}^{2}))).$$ (3.103)

In (3.102) and (3.103), for the sake of simplicity, we used the abbervations, $\varphi_{X}^k=\frac{\partial\varphi^{e}_k}{\partial X}$ for $X=\xi,\eta,\zeta$. It can easily been shown that the functions $\alpha_T(i,j,k)$ and $\beta_T(i,j,l,k)$ satisfy following relations,

$$\sum_{k=1}^4\alpha_T(i,j,k)=0, \sum_{i=1}^4\alpha_T(i,j,k)=K_{jk}(T), \sum_{j=1}^4\alpha_T(i,j,k)=K_{ik}(T),$$ (3.104)

$$\sum_{k=1}^4\beta_T(i,j,l,k)=0, \sum_{l=1}^4\beta_T(i,j,l,k)=0, \sum_{i,j=1}^4\beta_T(i,j,l,k)=K_{lk}(T).$$ (3.105)

We introduce now the backward Euler time discretization scheme for (3.90) and by applying the functions $\alpha_{T}(i,j,k)$ and $\beta_{T}(i,j,l,k)$ we write the $\ell_{k,1}$ functions (for $k=1,2,3,4$),

$$\ell_{k,1}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\frac{... ...I\frac{C_{I,i}^n}{C_{I}^{eq}}+(1-f_I)\frac{C_{V,i}^n}{C_{V}^{eq}}\Bigr)\Bigr\}+$$ (3.106)

$$+\frac{D_A }{\sqrt{4n_{i}^{2}+(C_{A,m}^n)^{2}}}\sum_{i,j,l=1}^4\B... ...I\frac{C_{I,i}^n}{C_{I}^{eq}}+(1-f_I)\frac{C_{V,i}^n}{C_{V}^{eq}}\Bigr)\Bigr\},$$ (3.107)

and residuum vector $R_{k,1}^e$ as

$$R_{k,1}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,1... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{A,i}^{n-1}.$$ (3.108)

Equations (3.74) and (3.75) of the Mulvaney-Richardson model have the same structure, it is sufficient, if we consider in more detail only (3.74),

$$\frac{\partial C_I}{\partial t}=D_I\Delta C_I - k_{f} ( C_I C_V - C_{I}^{eq} C_{V}^{eq}) +\nabla \cdot \mathbf{J}_{AI}.$$ (3.109)

We need to consider only $D_I\Delta C_I - k_{f} ( C_I C_V - C_{I}^{eq} C_{V}^{eq})$, the rest of the equation has already been treated above. Again, we use linear basis nodal functions $\varphi_k$ and apply Green's theorem,

$$-D_I\sum_{i,=1}^{4}C_{I,i}\int_T \nabla \varphi_i \nabla Nk d\Om... ... \varphi_j \varphi_kd\Omega + k_fC_{I}^{eq} C_{V}^{eq}\int_T \varphi_k d\Omega.$$ (3.110)

we set $\gamma_T(i,j,k)=\int_T \varphi_i \varphi_j \varphi_k d\Omega$, and proceed,

$$-D_I\sum_{i=1}^{4}C_{I,i} K_{i,k} - k_f\sum_{i,j=1}^{4}C_{I,i}C_{V,j}\gamma_T(i,j,k) + k_f C_{I}^{eq} C_{V}^{eq}V_k.$$ (3.111)

Finally, we have the discretizied form of equation (3.74),

$$\ell_{k,2}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\sum_{... ...j=1}^{4}C^{n}_{I,i}C_{V,j}^{n}\gamma_T(i,j,k) - C_{I}^{eq} C_{V}^{eq}V_k\Bigr)+$$

$$+\frac{f_I D_A}{C_{I}^{eq}}\sum_{i,j=1}^{4}\Bigl\{C^{n}_{A,i}\Big... ...{2}+(C^n_{A,m})^2}}\sum_{l=1}^{4}C_{I,l}^{n}\beta_T(i,j,l,k)\Bigr)\Bigr\},\quad$$ (3.112)

and the residuum vector $R_{k,2}^e$ as,

$$R_{k,2}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,2... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{I,i}^{n-1}.$$ (3.113)

and analogously for the vacancy balance equation,

$$\ell_{k,3}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\sum_{... ...j=1}^{4}C^{n}_{I,i}C_{V,j}^{n}\gamma_T(i,j,k) - C_{I}^{eq} C_{V}^{eq}V_k\Bigr)+$$

$$+\frac{(1-f_I) D_A}{C_{V}^{eq}}\sum_{i,j=1}^{4}\Bigl\{C^{n}_{A,i}... ...{2}+(C^n_{A,m})^2}}\sum_{l=1}^{4}C_{V,l}^{n}\beta_T(i,j,l,k)\Bigr)\Bigr\},\quad$$ (3.114)

and the corresponding residuum vector $R_{k,3}^e$ as,

$$R_{k,3}^e(\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)=\ell_{k,3... ...thbf{c}_I^n,\mathbf{c}_V^n)-\frac{1}{\Delta t}\sum_{i=1}^4 M_{ik}C_{V,i}^{n-1}.$$ (3.115)

Using the functions $\ell^e_{k,1}$, $\ell^e_{k,2}$, and $\ell^e_{k,3}$ for $k=1,2,3,4$ we construct the nucleus matrix $\mathbf{\Pi}(T)$ of the diffusion problem,

$$\mathbf{\Pi}(T)=\frac{\partial (\ell_{1,1}^e,\ell_{1,2}^e,\ell_{1... ...3,3}^e,\ell_{3,4}^e)}{\partial (\mathbf{c}_A^n,\mathbf{c}_I^n,\mathbf{c}_V^n)}.$$ (3.116)

Note that the global matrix $\mathbf{G}$ in this case has rank $3\cdot N$.

3.8.4 Analytical Solution of the Segregation Problem for One Dimension

In order to assess the numerical scheme for the problem described in Section 3.5.3 it is useful to construct an analytical solution for a special one-dimensional case. For the sake of simplicity, we consider in both segments intrinsic diffusion (Section 3.3.1). As segment $S_0$ the region $x>0$ is assumed and as segment $S_1$ region $x<0$. The one-dimensional segregation problem fullfils the following initial and interface conditions,

$$C_0(x,0) = C_{init} x>0 \quad C_1(x,0) = 0 x<0,$$ (3.117)

$$D_0 \frac{\partial C_0}{\partial x} = - h \Bigl(C_1-\frac{C_0}{m}\Bigr)$$ (3.118)

$$D_0 \frac{\partial C_0}{\partial x} = D_1 \frac{\partial C_1}{\partial x}, x=0.$$ (3.119)

Note that condition (3.117) also means that the media from both sides of the interface has an infinite length. We are searching for the solution of the problem given by (3.30), (3.117), (3.118) and (3.119) in the form,

$$x > 0, \quad C_0(x,t) = A_0 + B_0 C(x,t,\alpha_0,D_0)$$ (3.120)

$$x < 0, \quad C_1(x,t) = A_1 + B_1 C(-x,t,\alpha_1,D_1)$$ (3.121)

where $A_0, A_1, B_0, B_1, \alpha_0, \alpha_1$ are constants to be determined and

$$C(x,t,\alpha,D)= \Bigl (\frac{x}{2 \sqrt{D t}}\Bigr) - \exp\Bigl(\frac{h x \alpha + h^2 t \alpha^{2}}{D}\Bigr) \Bigl (\frac{x}{2 \sqrt{D t}} + \frac{h\sqrt{D t} \alpha}{D}\Bigr )$$ (3.122)

is a solution of the diffusion equation ( $\ref{eq:intrinsic-diffusion:1}$) for the case of the surface evaporation condition already studied in [38].
We determine constants $A_0, A_1, B_0, B_1, \alpha_0, \alpha_1$ from the initial and interface conditions (3.117), (3.118), (3.119) as follows. From the initial conditions we have

$$A_0 = C_{init} \quad A_1 = 0.$$ (3.123)

The interface condition (3.119) yields

$$\frac{D_0}{D_1}\frac{\frac{\partial C_0}{\partial x}}{\frac{\part... ...{D_0}\Bigr )}{\text{erfc} \Bigl (\frac{h\sqrt{D_1 t} \alpha_1}{D_1}\Bigr )}=-1.$$ (3.124)

This equation is fullfiled if

$$\frac{\alpha_0}{\sqrt{D_0}}=\frac{\alpha_1}{\sqrt{D_1}} \quad B_0\alpha_0 = -B_1\alpha_1.$$ (3.125)

From (3.118) and (3.123) follows,

$$\frac{h}{D_0}\frac{\Bigl(C_{1} - \frac{C_{0}} {m}\Bigr)}{\frac{\p... ...}{D_1}) \text{erfc} \Bigl (\frac{h\sqrt{D_1 t} \alpha_1}{D_1}\Bigr)\Bigr )-$$

$$-\frac{B_0}{m}\Bigl(1-\exp(\frac{h^2 t \alpha_0^{2}}{D_0}) \Bigl(\frac{h\sqrt{D_0 t} \alpha_0}{D_0}\Bigr )\Bigr)-\frac{C_{init}}{m}\Bigr)=1.$$ (3.126)

The last equality is ensured for the condition,

$$B_1-\frac{B_0}{m} = \frac{C_{init}}{m} \quad -B_1+\frac{B_0}{m} = B_0\alpha_0.$$ (3.127)

By solving the equation system given by (3.125) and (3.127) we have,

$$B_0 = -\frac{C_{init}}{1+m\sqrt{\frac{D_0}{D_1}}},\quad B_1 = \fr... ...+ \sqrt{\frac{D_1}{D_0}}},\quad\alpha_0 = \frac{1}{m} + \sqrt{\frac{D_0}{D_1}},$$

$$\alpha_1 = 1 + \frac{1}{m}\sqrt{\frac{D_1}{D_0}}.$$ (3.128)

So we can write a solution for the problem posed by (3.30), (3.117), (3.118) and (3.119). For $x>0$,

$$C_0(x,t) = C_{init} \Bigl ( 1- \Bigl ( \frac{1}{1+m\sqrt{\frac{D_0}{D_1}}}\Bigr)\Bigl ( \Bigl (\frac{x}{2 \sqrt{D_0 t}}\Bigr) - \exp\Bigl(\frac{h x \alpha_0 + h^2 t \alpha_{0}^{2}}{D_0}\Bigr) \Bigl (\frac{x}{2 \sqrt{D_0 t}} + \frac{h\sqrt{D_0 t} \alpha_1}{D_0}\Bigr ) \Bigr ) \Bigr),$$ (3.129)

and for $x<0$,

$$C_1(x,t) = \frac{C_{init}}{ m + \sqrt{\frac{D_1}{D_0}}}\Bigl ( \Bigl (-\frac{x}{2 \sqrt{D_1 t}} \Bigr) -\exp\Bigl(\frac{-h x \alpha_1 + h^2 t \alpha_{1}^{2}}{D_1} \Bigr ) \Bigl (-\frac{x}{2 \sqrt{D_1 t}} + \frac{h\sqrt{D_1 t} \alpha_1}{D_1}\Bigr ) \Bigr ).$$ (3.130)

3.8.5 Numerical Handling of the Segregation Model

Let us assume that the segments $S_0$ and $S_1$ are comprised of three-dimensional areas $\Omega_0$ and $\Omega_1$ and the connecting interface $I$ with two-dimensional area $\Theta$, respectively. The tetrahedralization of areas $\Omega_0$, $\Omega_1$ and the triangulation of area $\Theta$ are denoted as $T_h(\Omega_0)$, $T_h(\Omega_1)$ and $T_h(\Theta)$. We discretize (3.30) on the element $T\in T_h(\Omega_0)$ using a linear basis function $\varphi_k$. After introducing the weak formulation of the equation and subsequently applying Green's theorem we have

$$\int\limits_T\frac{\partial C} {\partial t} \varphi_k d\Omega =... ...l n} \varphi_k d\Gamma -D_0\int\limits_T \nabla C \nabla \varphi_k d\Omega,$$ (3.131)

where $\Gamma$ is the boundary of the element $T$. Assuming that $T \cap T_h(\Theta) = \Gamma_{\Theta}\neq\emptyset$ and marking all inside faces of $T$ as $\Gamma_{in}$ we have $\Gamma = \Gamma_{\Theta} \cup \Gamma_{in}$ and we can write

$$\int\limits_\Gamma\frac{\partial C} {\partial n} \varphi_k d\Ga... ...+ \int\limits_{\Gamma_{in}}\frac{\partial C} {\partial n} \varphi_k d\Gamma.$$ (3.132)

In the standard finite element assembling procedure [5], we take into account only the terms $\int_T \nabla C \nabla \varphi_k d\Omega$ when building up the stiffness matrix. Thereby the terms $\int_{\Gamma_{in}}\frac{\partial C} {\partial n} \varphi_k d\Gamma$ do not need to be considered because of their annihilation on the inside faces.
The term $\int_{\Gamma_{\Theta}}\frac{\partial C} {\partial n} \varphi_k d\Gamma$ makes sense only on the interface area $\Theta$ and there it can be used to introduce the influence of the species flux from the neighboring segment area $\Omega_1$ by applying the segregation flux formula (3.65)

$$\int\limits_{\Gamma_{\Theta}}\frac{\partial C} {\partial n} \var... ...s_{\Gamma_{\Theta}} h\Bigl(C^{1} - \frac{C^{0}} {m}\Bigr) \varphi_k d\Gamma.$$ (3.133)

After a usual assembling procedure on the tetrahedralization $T_h(\Omega_0)$ and $T_h(\Omega_1)$ has been carried out and the global stiffness matrix for both segment areas of the segregation problem has been built, the interface inputs (3.131) for the segregation fluxes $J_{01}$ and $J_{10}$ are evaluated on the triangulation $T_h(\Theta)$ and assembled into the global stiffness matrix according to the particular assembling algorithm developed in this work.

The numerical algorithm based on the concept described above is carried out in two steps.
Step 1.
We assemble the general matrix $\mathbf{G}$ of the problem for both segments, i.e., for both diffusion processes. The number of points in segments $S_0$ and $S_1$ is denoted as $s_0$ and $s_1$, respectively. The general matrix has dimensions $(s_0+s_1) \times (s_0+s_1)$ and the inputs are correspondingly indexed. The matrix is assembled by distributing the inputs from matrix $\mathbf{\Pi_i}(T_i)$, $dim(\mathbf{\Pi_i}(T_i))= 4 \times 4$, defined for each $T_i\in T_h(\Omega_i), i=0,1$

$$\mathbf{\Pi_i}(T_i) = \mathbf{K}(T) + D_i \Delta t\mathbf{M}(T_i)$$ (3.134)

The $\Delta t$ is the time step of the discretisized time and $\mathbf{K}(T_i)$ and $\mathbf{M}(T_i)$ are stiffness and mass matrix defined on single tetrahedra $T_i$ from $T_h(\Omega_i)$

$$K_{pq}(T_i)=\int\limits_T \nabla N_p \nabla N_q dx dy dz M_{pq}(T_i)=\int\limits_T N_p N_q dx dy dz\quad p,q\in\{0,1,2,3\}.$$

Let us denote vertices of the element $T_i$ from the tetrahedralization $T_h(\Omega_i)$ by $P_0, P_1, P_2, P_3$ and their indexes in the segment $S_i$ by $0\leqslant k^{i}_{P_{0}}, k^{i}_{P_{1}}, k^{i}_{P_{2}}, k^{i}_{P_{3}}<s_i$. Assembling means, for each $T_i \in T_h(\Omega_i)$, and $r,q\in\{0,1,2,3\}$ adding the term $\mathbf{\Pi_i}(r,q)$ to $\mathbf{G}(i s_0+k^{i}_{P_{r}}, i s_0+k^{i}_{P_{q}})$.
After this assembling procedure is carried out the general matrix $\mathbf{G}$ has the following structure,

$$\mathbf{G}= \begin{bmatrix}\mathbf{S_{00}} & \mathbf{0} \mathbf{0} & \mathbf{S_{11}} \end{bmatrix}$$ (3.135)

Where $dim(\mathbf{S_{00}}) = s_0 \times s_0$ and $dim(\mathbf{S_{11}}) = s_1 \times s_1$. The matrix $\mathbf{S_{00}}$ and $\mathbf{S_{11}}$ are the finite element discretizations of equation (3.30) for the segments $S_0$ and $S_1$, respectively.
Step 2.
For the element $T_0 \in T_h(\Omega_0)$ with one of it's faces ( $\Gamma_{\Theta}$) laying on the interface $I$ according to the idea presented above, weak formulation is,

$$\int\limits_{T_{0}}\frac{\partial C} {\partial t} N_k d\Omega =... ...m}\Bigr) N_k d\Gamma - D_0\int\limits_{T_{0}} \nabla C \nabla N_k d\Omega.$$ (3.136)

The segregation term on the right side of ( $\ref{SegrWeak}$) is evaluated on the two-dimensional element $\Gamma_{\Theta}$. If we discretisize equation ( $\ref{SegrWeak}$) and take a backward Euler time scheme with time step $\Delta t$, we obtain for segment $S_0$,

$$\mathbf{M}(T_0) (C_0^n-C_0^{n-1}) =h\Delta t \mathcal{M} (C_1^n-\frac{1}{m}C_0^n) - D_0 \Delta t \mathbf{K}(T_0) C_0^n,$$ (3.137)

where

$$C_0^n=(C_{0,P_0}^n, C_{0,P_1}^n, C_{0,P_2}^n, C_{0,P_3}^n)^T, C_0^{n-1}=(C_{0,P_0}^{n-1}, C_{0,P_1}^{n-1}, C_{0,P_2}^{n-1}, C_{0,P_3}^{n-1})^T,$$ (3.138)

are the values of the species concentration for the $n^{th}$ and $n-1^{st}$ time step at the vertices of element $T_0$ and analogously $C_1^n=(C_{1,P_0}^n, C_{1,P_1}^n, C_{1,P_2}^n, C_{1,P_3}^n)^T$ for $T_1\in T_h(\Omega_1)$, $T_0 \cap T_1 = \Gamma_{\Theta}$. Without loosing generality we assume that verticex $P_3$ of the tetrahedras $T_0$ and $T_1$ is the point which doesn't belong to the interface $\Theta$. In that case the matrix $\mathcal{M}$ from ( $\ref{FirstDisc}$) has a simple structure,

$$\mathcal{M} = \frac{det(J(\Gamma_{\Theta}))}{24} \begin{bmatrix}2 & 1 & 1 & 0 1 & 2 & 1 & 0 1 & 1 & 2 & 0 0 & 0 & 0 & 0 \end{bmatrix},$$ (3.139)

where $J(\Gamma_{\Theta})$ is the Jacobi matrix evaluated on element $\Gamma_{\Theta}$. Using (3.132) we have,

$$\mathbf{\Pi_0}(T_0) C_0^{n} - h\Delta t \mathcal{M} (C_1^n-\frac{1}{m}C_0^n) = \mathbf{M}(T_0) C_0^{n-1}.$$ (3.140)

We introduce now $\mathbf{H_0}(\Gamma_{\Theta})= -\frac{h}{m}\Delta t \mathcal{M}$ and $\mathbf{H_1}(\Gamma_{\Theta})= h\Delta t \mathcal{M}$ and write for the element $T_0$,

$$\mathbf{\Pi_0}(T_0) C_0^n - \mathbf{H_0}(\Gamma_{\Theta}) C_0^n - \mathbf{H_1}(\Gamma_{\Theta}) C_1^n = \mathbf{M}(T_0) C_0^{n-1},$$ (3.141)

and analogously for the element $T_1$,

$$\mathbf{\Pi_1}(T_1) C_1^n + \mathbf{H_0}(\Gamma_{\Theta}) C_0^n + \mathbf{H_1}(\Gamma_{\Theta}) C_1^n = \mathbf{M}(T_1) C_1^{n-1}.$$ (3.142)

In the following, for the sake of simplicity, we omit $\Gamma_{\Theta}$ from $\mathbf{H_i}(\Gamma_{\Theta})$ and write $\mathbf{H_i}$.

The contributions of $\mathbf{\Pi_0}(T_0)$ and $\mathbf{\Pi_1}(T_1)$ are already included in the general matrix $\mathbf{G}$ by the assembling procedure made in the first step, the build up of $\mathbf{G}$ can now be completed by adding the inputs from matrix $\mathbf{H_0}$ and $\mathbf{H_1}$ in order to take into account the segregation effect on interface $I$.
We denote the vertices of the element $\Gamma_{\Theta} \in T_h(\Theta)$ as $P_0, P_1, P_2$. In the tetrahedralization $T_h(\Omega_0)$ these points have indices $0\leqslant k^{0}_{P_{0}}, k^{0}_{P_{1}}, k^{0}_{P_{2}} < s_0$ and indices $0\leqslant k^{1}_{P_{0}}, k^{1}_{P_{1}}, k^{1}_{P_{2}} <s_1$ in the tetrahedralization $T_h(\Omega_1)$. The actual implementation of the scheme ( $\ref{SegrBasis}$) is for each element $\Gamma_{\Theta}$ of the interface triangulation $T_h(\Theta)$ and for $r,q\in\{0,1,2,3\}$:

• adding the term $-\mathbf{H_0}(r,q)$ to the input $\mathbf{G}(k^{0}_{P_{r}}, k^{0}_{P_{q}})$
• adding the term $-\mathbf{H_1}(r,q)$ to the $\mathbf{G}(k^{0}_{P_{r}}, s_0 + k^{1}_{P_{q}})$
• adding the term $\mathbf{H_0}(r,q)$ to $\mathbf{G}(s_0+k^{1}_{P_{r}}, k^{0}_{P_{q}})$
• adding the term $\mathbf{H_1}(r,q)$ to the $\mathbf{G}(s_0 + k^{1}_{P_{r}}, s_0+k^{1}_{P_{q}})$.

Now the assembling of the general matrix $\mathbf{G}$ is completed. These procedure is carried out for each time step of the simulation. The numerical solution of the segregation problem is compared with an analytical solution on Figure 3.2.

Figure 3.2: The comparison between numerical (points) and analytical solution (full line) for three different time steps in each figure above. Two different spatial discretizations (N=20 and N=80 points) are used. Because of the assumption of the infinite media by the derivation of analytical solution there is a deviation of this solution from the numerical one in the proximity of the point $-10$ on the abscissa. Numerical solution considers finite segment with zero Neumann boundary condition at the $-10$ point.

$N=20$, $m=3$

$N=80$, $m=3$

$N=20$, $m=0.2$

$N=80$, $m=0.2$