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,$$ (194)

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.$$ (195)

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.$$ (196)

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.133) 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)$$ (197)

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\}.$$ (198)

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}$$ (199)

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.$$ (200)

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,$$ (201)

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,$$ (202)

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},$$ (203)

where $J(\Gamma_{\Theta})$ is the Jacobi matrix evaluated on element $\Gamma_{\Theta}$. Using (3.134) 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}.$$ (204)

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},$$ (205)

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}.$$ (206)

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.

$N=20$, $m=3$

$N=80$, $m=3$

Figure: 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=0.2$

$N=80$, $m=0.2$