Interdiffusion in many dimensions: mathematical models,numerical simulations and experiment

Abstract

Over the last two decades, there have been tremendous advances in the computation of diffusion and today many key properties of materials can be accurately predicted by modelling and simulations. In this paper, we present, for the first time, comprehensive studies of interdiffusion in three dimensions, a model, simulations and experiment. The model follows from the local mass conservation with Vegard’s rule and is combined with Darken’s bi-velocity method. The approach is expressed using the nonlinear parabolic–elliptic system of strongly coupled differential equations with initial and nonlinear coupled boundary conditions. Implicit finite difference methods, preserving Vegard’s rule, are generated by some linearization and splitting ideas, in one- and two-dimensional cases. The theorems on the existence and uniqueness of solutions of the implicit difference schemes and the consistency of the difference methods are studied. The numerical results are compared with experimental data for a ternary Fe-Co-Ni system. A good agreement of both sets is revealed, which confirms the strength of the method.

1. Introduction

Diffusion controls many chemical and physical processes; therefore, a knowledge of diffusion and the ability to study it theoretically and experimentally are very valuable in materials science. Common examples of diffusion-controlled processes are: sintering of powders, grain growth, nitriding and carburizing by vapour deposition, oxidation and diffusional creep, as well as processes applied in microelectronics devices, such as the addition of n or p impurities to semiconducting silicon, germanium or other crystals, welding, crystallization and many other alloying and metallurgical processes. Among the given examples of engineering applications, there are situations in which the solids under study do not present one-dimensional geometry. These are mainly examples of metallic joints in welding and soldiering, and of the stability of coating–substrate systems. Therefore, a knowledge of the mass transport in bodies with an arbitrary two- or three-dimensional geometry is important.

A suitable formalism for an $s$ -component interdiffusion, appearing as a movement of solute species when two solids are in contact, is based on a fundamental physical law expressing mass conservation

$${\partial }_t{c}_i=-\mathrm{div}{J}_i,i=1,\dots ,s,$$

(1)

where $c_i$ and $J_i$ denote the concentration and flux of the $i$ th component. Very often, a suitable formalism is based on the mass conservation law with the reaction term (sources or sinks) $A_i=A_i(c_1,\dots ,c_s)$ ,

$${\partial }_t{c}_i=-\mathrm{div}{J}_i+A_i,i=1,\dots ,s.$$

(2)

Mass conservation needs to be combined with a constitutive equation for the flux. A theory that is considered the most fundamental to describe interdiffusion in multicomponent systems is provided by linear phenomenological laws [1] combined with Onsager reciprocal relations [2, 3]. Linear phenomenological laws predict that when the system is close to the state of equilibrium and local equilibrium holds, the flows $J_k$ , $k=1,\dots ,s$ are expected to be linear combinations of all forces $X_j$ , $j=1,\dots ,s$ acting in the system

$$J_k=\sum _{j=1}^s{L}_{\mathrm{kj}}{X}_j,k=1,\dots ,s.$$

(3)

In this equation, $L_{\mathrm{kj}}$ , $k,j=1,\dots ,s$ are referred to as phenomenological transport coefficients. They form a two-dimensional matrix and obey the Onsager reciprocal relations $L_{\mathrm{kj}}=L_{\mathrm{jk}}$ , which imply so-called cross-effects. Reciprocal relations occur between different pairs of the forces and flows [2, 4–7]. The knowledge of the symmetric matrix of phenomenological coefficients is often somewhat ambiguous. In computations, one often uses the bi-velocity method, which combines Darken’s and Brenner’s concepts that the volume velocity is essential in defining the local material velocity at nonequilibrium, $v^{\mathrm{d}}$ , the so-called drift or Darken velocity. In the bi-velocity (Darken) method, two fluxes and two velocities are applied for each component to quantify the Kirkendall shift in solid solutions [8–13]

$$J_i=J_i^{\mathrm{diff}}+J_i^{\mathrm{d}}$$

(4)

and

$$v_i=v_i^{\mathrm{diff}}+v^{\mathrm{d}},i=1,\dots ,s,$$

(5)

where the superscripts “ $\mathrm{diff}$ ” and “ $\mathrm{d}$ ” relate to the diffusion and drift (Darken) fluxes or velocities, respectively. We postulate the incompressible transport (Vegard’s rule [14])

$$\sum _{i=1}^s{\mathrm{\Omega }}_i{c}_i=1,$$

(6)

where ${\mathrm{\Omega }}_i$ means the partial molar volume of the $i$ th component. The drift velocity is generated by the volume continuity equation [15]

$$\mathrm{div}(v^{\mathrm{d}}+\sum _{i=1}^s{\mathrm{\Omega }}_i{J}_i^{\mathrm{diff}})=0,$$

(7)

which is implied by equations (1), (4) and (6). In the one-dimensional case, if the fluxes on the boundaries of a domain equal zero, then the explicit formula on $v^{\mathrm{d}}$ follows immediately from equation (7),

$$v^{\mathrm{d}}=-\sum _{i=1}^s{\mathrm{\Omega }}_i{J}_i^{\mathrm{diff}}.$$

(8)

However, in the multidimensional cases, this is obviously not true. Nonetheless, the three-dimensional case has been studied, for example in [16, 17] with the use of equation (8). Taking into account equation (7), this is not a mathematical method and, physically, a part of information from the fluxes is neglected. We omit this assumption by a postulate that the drift field is a potential, i.e. $v^{\mathrm{d}}=-\nabla F$ . An attempt to construct a model with the potential drift field has been made in [18]. Unfortunately, in that study, the idea of deriving formulae is not correct because the rotation of a scalar potential is considered. Moreover, the uniqueness condition is not given. There are also no proofs of existence and uniqueness of numerical solutions, nor is there a comparison with any real experiment.

Much research has been conducted to understand the interdiffusion behaviour and its after-effects in multicomponent systems. For an extensive study, see, for example, the book by Gusak [19]. If a process is described by the diffusive model, the diffusion equation must be solved and initial and boundary conditions should be established. For two- and three-dimensional diffusion, a geometry of the solids must be assumed. Owing to many difficulties, numerical solution of the diffusion equation in three dimensions remains a challenging mathematical problem.

There have been only a few attempts to study multicomponent diffusion in the systems of an arbitrary geometry. Verbrugge et al. [20] simulated Fickian diffusion ( $s=1$ ) and stress formation in a spherical core-shell system, as in lithium-ion batteries. Experimentally, the interdiffusion and migration of associated voids, as manifested in bimetallic core-shell nanoparticles, were studied by Chee et al. [21]. The diffusion with drift in a multicomponent mixture was studied by Wierzba [22]. As a result, a good agreement of the experiment and calculated evolution of voids in two-dimensional space was achieved.

To make progress, we combine, in the present studies, computations and experiments of the interdiffusion in two-dimensional geometry. As an example, a ternary Fe-Co-Ni system is considered. A model based on Darken’s method with the use of Vegard’s rule is formulated. It is expressed by the strongly coupled nonlinear parabolic–elliptic system

$$\{\begin{array}{c}{\partial }_t{c}_i+\mathrm{div}(-D_i(c_1,\dots ,c_s)\nabla c_i-c_i\nabla F)=0,\hfill \\ \multicolumn{1}{c}{-\mathrm{\Delta }F=\mathrm{div}(\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(c_1,\dots ,c_s)\nabla c_k),}\\ \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}F\mathrm{d}x=0,}\end{array}$$

(9)

for $i=1,\dots ,s$ , with the initial and the nonlinear coupled boundary conditions. Such a model was introduced previously [23, 24]. Implicit finite difference methods (FDM), preserving Vegard’s rule, generated by some linearization and splitting ideas are used to present numerical solutions. Numerical methods that preserve some futures of a continuous model are very important and interesting from a physical and mathematical point of view. For example, boundedness-, positivity- or monotonicity-preserving numerical methods were studied in [25–35] and in the references therein. We prove that the implicit difference schemes given have a unique solution and that the difference methods are consistent. This paper extends the results in [24], where only the case of constant diffusion coefficients is studied. In our paper, these are nonlinear; more precisely, they are functions of mole fractions. A real experiment is made and the results are compared with our numerical datasets. The comparison presented in Section 4 shows that our difference methods are convergent and stable. Since the differential system (equation (9)) is a parabolic–elliptic strongly nonlinear system, theoretical convergence and stability analysis for the proposed numerical methods is a challenging task. In our previous paper [24], where the diffusion coefficients were constant, we have proved that the convergence and stability with respect to the concentrations were equivalent to the convergence and stability with respect to the potential, for suitable difference maximum norms. In the proofs, we have used a technique studied in [33–35]. Unfortunately, we have not yet been able to use this technique in the case of the nonlinear diffusion coefficients considered in this paper. The same problem is open for the well-known parabolic–elliptic Nernst–Planck–Poisson model, even in the case of constant diffusion coefficients.

Let us stress that such strongly coupled systems as in our model (i.e., by the second derivatives) have seldom been investigated mathematically. Equation (9) is a little similar in the structure to the Nernst–Planck–Poisson system [36–39]. We suppose, after analysis of the papers mentioned, concerning the Nernst–Planck–Poisson model, that if the initial concentrations and the fluxes on the boundaries are sufficiently regular, then equation (9) has a unique regular solution. It is a completely new model and our analytical investigations are not finished.

The paper is organized in the following way. In Section 2, the initial–boundary differential problem is formulated. Section 3 deals with the construction of the implicit finite difference methods, the analysis of the existence and the uniqueness of solutions to the suitable difference schemes and consistency of the difference methods in the one- and two-dimensional cases, respectively. Experimental and simulated data are discussed in Section 4. We show preparation of the diffusion couple, diffusion experiment, concentration maps and chosen profiles. The diffusion fluxes and drift velocity are computed and the time evolution is simulated. In Section 5, the results are summarized.

2. Interdiffusion model

In [23, 24], the following model of interdiffusion is constructed. Let an open and bounded set $\mathrm{\Omega }\subset {\mathbb{R}}^n$ with the smooth boundary $\partial \mathrm{\Omega }$ be given. The processing time $t_{\mathrm{end}}>0$ and the number of components $s\in \mathbb{N}\backslash \{1\}$ are fixed, and denote ${\mathbb{R}}_{+}=(0,\infty )$ . We assume that the drift velocity $v^{\mathrm{d}}:[0,t_{\mathrm{end}}]\times \mathrm{\Omega }\to {\mathbb{R}}^n$ is potential, i.e. there exists $F:[0,t_{\mathrm{end}}]\times \mathrm{\Omega }\to \mathbb{R}$ such that

$$-\nabla F=v^{\mathrm{d}}.$$

(10)

The following data are given:

• (a) ${\mathrm{\Omega }}_i=\mathrm{const}\in {\mathbb{R}}_{+}$ – the partial molar volume of the $i$ th component of the mixture, $i=1,\dots ,s$ ;

• (b) $D_i:[0,1/{\mathrm{\Omega }}_1]\times \dots \times [0,1/{\mathrm{\Omega }}_s]\to {\mathbb{R}}_{+}$ – the diffusion coefficient of the $i$ th component of the mixture, $i=1,\dots ,s$ ;

• (c) $c_{0i}:\mathrm{\Omega }\to {\mathbb{R}}_{+}$ – the initial concentration of the $i$ th component of the mixture, $i=1,\dots ,s$ ;

• (d) $j_i:[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega }\to \mathbb{R}$ – the evolution of a mass flux of the $i$ th component of the mixture through the boundary $\partial \mathrm{\Omega }$ , $i=1,\dots ,s$ .

The following functions are unknown:

• (a) $c_i:[0,t_{\mathrm{end}}]\times \mathrm{\Omega }\to {\mathbb{R}}_{+}$ – the concentration of the $i$ th component of the mixture, $i=1,\dots ,s$ ;

• (b) $F:[0,t_{\mathrm{end}}]\times \mathrm{\Omega }\to \mathbb{R}$ – the potential of the drift velocity.

We assume that each component of the mixture is a continuous medium, i.e. it satisfies the local mass conservation law (continuity equation)

$${\partial }_t{c}_i+\mathrm{div}{J}_i=0\mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },i=1,\dots ,s,$$

(11)

where

$$J_i=-D_i(c_1,\dots ,c_s)\nabla c_i-c_i\nabla F$$

(12)

is a flux of the $i$ th component of the mixture, and it is a sum of the diffusional and Darken drift fluxes (see equation (4)). Moreover, we postulate the constant partial molar volumes ${\mathrm{\Omega }}_i$ , i.e. the concentrations invariant and the incompressible transport (Vegard’s rule)

$$\sum _{i=1}^s{\mathrm{\Omega }}_i{c}_i=1\mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega }.$$

(13)

The solutions of the interdiffusion in the $s$ -component solid are the functions $c_i=c_i(t,x)$ and the potential of the drift velocity $F=F(t,x)$ , which fulfil the strongly coupled nonlinear parabolic–elliptic system

$$(\begin{array}{cc}{\partial }_t{c}_i+\mathrm{div}(-D_i(c_1,\dots ,c_s)\nabla c_i-c_i\nabla F)=0\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{-\mathrm{\Delta }F=\mathrm{div}(\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(c_1,\dots ,c_s)\nabla c_k)}& \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}F\mathrm{d}x=0}& \mathrm{on}[0,t_{\mathrm{end}}],\hfill \end{array}$$

(14)

with the initial conditions

$$c_i(0,x)=c_{0i}(x)\mathrm{on}\mathrm{\Omega },$$

(15)

and the coupled nonlinear boundary conditions

$$(\begin{array}{cc}-D_i(c_1,\dots ,c_s)\frac{\partial c_i}{\partial \mathbf{n}}-c_i\frac{\partial F}{\partial \mathbf{n}}=j_i(t,x)\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{\frac{\partial F}{\partial \mathbf{n}}=-\sum _{k=1}^s{\mathrm{\Omega }}_k(D_k(c_1,\dots ,c_s)\frac{\partial c_k}{\partial \mathbf{n}}+j_k(t,x))}& \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \end{array}$$

(16)

where $\mathbf{n}$ is the outside normal to the boundary $\partial \mathrm{\Omega }$ , $i=1,\dots ,s$ .

We assume the constant overall volume; thus, from incompressible transport

$${\int }_{\partial \mathrm{\Omega }}\sum _{i=1}^s{\mathrm{\Omega }}_i{j}_i(t,x)\mathrm{d}S=0\mathrm{on}[0,t_{\mathrm{end}}].$$

(17)

This can be treated as the compatibility condition to the elliptic subsystem to be well posed. This condition also follows immediately from the Gauss theorem.

Moreover, we assume the Vegard rule on the initial concentrations

$$\sum _{i=1}^s{\mathrm{\Omega }}_i{c}_{0i}(x)=1\mathrm{on}\mathrm{\Omega }.$$

(18)

Remark 1. In the case $n=1$ , equation (14), with the boundary conditions (equation (16)), implies, by integration of the elliptic equation on $F$ , the strongly coupled nonlinear evolution system

$${\partial }_t{c}_i+{\partial }_x(-D_i(c_1,\dots ,c_s){\partial }_x{c}_i+c_i(\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(c_1,\dots ,c_s){\partial }_x{c}_k+K(t)))=0$$

(19)

on $[0,t_{\mathrm{end}}]\times \mathrm{\Omega }$ and the coupled nonlinear boundary conditions

$$-D_i(c_1,\dots ,c_s)\frac{\partial c_i}{\partial \mathbf{n}}+c_i(\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(c_1,\dots ,c_s)\frac{\partial c_k}{\partial \mathbf{n}}+K(t)\mathbf{n})=j_i(t,x)$$

(20)

on $[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega }$ , $i=1,\dots ,s$ , where

$$K(t)=\sum _{i=1}^s{\mathrm{\Omega }}_i{j}_i(t,\mathrm{\Lambda })=-\sum _{i=1}^s{\mathrm{\Omega }}_i{j}_i(t,-\mathrm{\Lambda })\mathrm{on}[0,t_{\mathrm{end}}],$$

$\mathrm{\Omega }=(-\mathrm{\Lambda },\mathrm{\Lambda })$ . The one-dimensional model (equations (19), (15) and (20)) has been well-known for many years [11, 23, 40–42].

Remark 2. The general situation

$$J_i=-D_i(c_1,\dots ,c_s)\nabla c_i+c_i{v}^{\mathrm{d}}$$

leads to a differential algebraic system like equations (11) and (13), with $s+1$ equations and $s+n$ unknowns. It is obviously badly posed in the multidimensional case, i.e. $n>1$ . If $n=2$ or $n=3$ , $\mathrm{\Omega }$ is a simply connected region and $\mathrm{rot}{v}^{\mathrm{d}}=0$ , then there exists $F$ in equation (10). In solids, viscosity is very high and unknown in most cases and the diffusion processes in solids show negligible turbulence. Thus, it can be neglected and the postulate about $\mathrm{rot}{v}^{\mathrm{d}}=0$ is appropriate. Under equation (10), a number of equations and unknowns in the differential algebraic system of equations (11) and (13) is the same; it is $s+1$ for each $n$ . However, from a numerical and analytical point of view, it is better to study the equivalent differential system (equation (14)). For details of its construction we refer the readers to [23, 24].

3. Implicit difference method

3.1. Case n = 1

Let $\mathrm{\Omega }=(-\mathrm{\Lambda },\mathrm{\Lambda })$ . We assume that each diffusion coefficient $D_i$ depends on the mole fraction $N_i=c_i/(\sum _{k=1}^s{c}_k)$ of the $i$ th component of the mixture and the physical system is closed, i.e. $j_i(t,\mathrm{\Lambda })=j_i(t,-\mathrm{\Lambda })\equiv 0$ for $t\in [0,t_{\mathrm{end}}]$ , $i=1,\dots ,s$ . We will write, for simplicity, $D_i(N_i)$ , where $D_i:[0,1]\to {\mathbb{R}}_{+}$ . The initial–boundary differential problem (equations (14) to (16)) now takes the form

$$(\begin{array}{cc}{\partial }_t{c}_i={\partial }_x\left(D_i\left(N_i\right){\partial }_x{c}_i\right)+{\partial }_x{c}_i{\partial }_{x}F+c_i{\partial }_{\mathrm{xx}}F\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{-{\partial }_{\mathrm{xx}}F=\sum _{k=1}^s{\mathrm{\Omega }}_k{\partial }_x\left(D_k\left(N_k\right){\partial }_x{c}_k\right)}& \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}F\mathrm{d}x=0}& \mathrm{on}[0,t_{\mathrm{end}}],\hfill \end{array}$$

(21)

$$c_i(0,x)=c_{0i}(x)\mathrm{on}\mathrm{\Omega },$$

(22)

$$(\begin{array}{cc}-D_i(N_i){\partial }_x{c}_i+c_i\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(N_k){\partial }_x{c}_k=0\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\partial }_{x}F=-\sum _{k=1}^s{\mathrm{\Omega }}_k{D}_k(N_k){\partial }_x{c}_k}& \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \end{array}$$

(23)

$i=1,\dots ,s$ . Taking into account the Vegard rule (equation (13)), we calculate

$$c_s=\frac{1}{{\mathrm{\Omega }}_s}(1-\sum _{i=1}^{s-1}{\mathrm{\Omega }}_i{c}_i),$$

(24)

and then equations (21) to (23) are reduced as follows:

$$(\begin{array}{cc}{\partial }_t{c}_i={\partial }_x(D_i(N_i){\partial }_x{c}_i)+{\partial }_x{c}_i{\partial }_{x}F+c_i{\partial }_{\mathrm{xx}}F\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{-{\partial }_{\mathrm{xx}}F=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k{\partial }_x\left((D_k(N_k)-D_s(N_s)){\partial }_x{c}_k\right)}& \mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}F\mathrm{d}x=0}& \mathrm{on}[0,t_{\mathrm{end}}],\hfill \end{array}$$

(25)

$$c_i(0,x)=c_{0i}(x)\mathrm{on}\mathrm{\Omega },$$

(26)

$$(\begin{array}{cc}-D_i(N_i){\partial }_x{c}_i+c_i\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_k(N_k)-D_s(N_s)){\partial }_x{c}_k=0\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\partial }_{x}F=-\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_k(N_k)-D_s(N_s)){\partial }_x{c}_k}& \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \end{array}$$

(27)

$i=1,\dots ,s-1$ . We assume additionally that

$$\sum _{i=1}^{s-1}{\mathrm{\Omega }}_i{c}_{0i}(x)\le 1\mathrm{on}\mathrm{\Omega }.$$

(28)

Define a mesh on $\mathrm{\Omega }$ in the following way. Let $h=2\mathrm{\Lambda }/(M+1)$ and $\tau =t_{\mathrm{end}}/K$ stand for the space and time steps of the mesh, respectively, where $M,K\in \mathbb{N}$ are given. Define nodal points $t^{\mu }$ , $x_m$ , $x_{m+\frac{1}{2}}$ as follows:

$$\begin{array}{c}{t}^{\mu }=\mu \tau ,\mu =0,1,\dots ,K,\hfill \\ x_m=-\mathrm{\Lambda }+\mathrm{mh},m=0,1,\dots ,M+1,\hfill \\ x_{m+\frac{1}{2}}=-\mathrm{\Lambda }+(m+\frac{1}{2})h,m=0,1,\dots ,M.\hfill \end{array}$$

The terms ${\partial }_x(D_i(N_i){\partial }_x{c}_i)$ in equation (25) will be approximated at a point $(t^{\mu },x_m)$ with the use of a modification of the method given in [43], namely by the difference quotients

$$\frac{1}{h}(D_{i,m+\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m+1}^{\mu +1}-c_{i,m}^{\mu +1})-D_{i,m-\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m}^{\mu +1}-c_{i,m-1}^{\mu +1})),$$

where

$$\begin{array}{c}{D}_{i,m-\frac{1}{2}}^{\mu }=D_i\left(\frac{1}{2}(N_{i,m-1}^{\mu }+N_{i,m}^{\mu })\right),\hfill \\ D_{i,m+\frac{1}{2}}^{\mu }=D_i\left(\frac{1}{2}(N_{i,m}^{\mu }+N_{i,m+1}^{\mu })\right),\hfill \end{array}$$

for $m=1,\dots ,M$ . Analogously, the terms ${\partial }_x((D_k(N_k)-D_s(N_s)){\partial }_x{c}_k)$ will be approximated, but the concentrations will be taken at points $t^{\mu }$ instead of $t^{\mu +1}$ . To approximate the remaining derivatives, the central, forward and back difference quotients will be used. Moreover, the terms $D_i(N_i)$ in the first $(s-1)$ equations in equation (27) will be approximated by

$$D_{i,m}^{\mu }=D_i\left(N_{i,m}^{\mu }\right)$$

for $m=0,M+1$ and similarly $D_k(N_k)$ , $D_s(N_s)$ . But the terms $D_k(N_k)$ , $D_s(N_s)$ in the last equation in equation (27) will be approximated by $D_{k,1/2}^{\mu }$ , $D_{s,1/2}^{\mu }$ , $D_{k,M+1/2}^{\mu }$ , $D_{s,M+1/2}^{\mu }$ . We will find an approximate solution of equations (25) to (27) using a splitting linearized implicit difference method.

We define an implicit difference scheme for the elliptic subsystem on the potential $F$ in equations (25) to (27). It is a system of linear algebraic equations of the form

$$(\begin{array}{c}{F}_0^{\mu +1}-F_1^{\mu +1}=P_0^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,\frac{1}{2}}^{\mu }-D_{s,\frac{1}{2}}^{\mu })(c_{k,1}^{\mu }-c_{k,0}^{\mu }),\hfill \\ \multicolumn{1}{c}{-F_{m-1}^{\mu +1}+2F_m^{\mu +1}-F_{m+1}^{\mu +1}=P_m^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k((D_{k,m-\frac{1}{2}}^{\mu }-D_{s,m-\frac{1}{2}}^{\mu })c_{k,m-1}^{\mu }}\\ \multicolumn{1}{c}{-(D_{k,m+\frac{1}{2}}^{\mu }-D_{s,m+\frac{1}{2}}^{\mu }+D_{k,m-\frac{1}{2}}^{\mu }-D_{s,m-\frac{1}{2}}^{\mu })c_{k,m}^{\mu }+(D_{k,m+\frac{1}{2}}^{\mu }-D_{s,m+\frac{1}{2}}^{\mu })c_{k,m+1}^{\mu }),}\\ \multicolumn{1}{c}{-F_M^{\mu +1}+F_{M+1}^{\mu +1}=P_{M+1}^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+\frac{1}{2}}^{\mu }-D_{s,M+\frac{1}{2}}^{\mu })(c_{k,M}^{\mu }-c_{k,M+1}^{\mu }),}\\ \multicolumn{1}{c}{F_0^{\mu +1}+2\sum _{m=1}^M{F}_m^{\mu +1}+F_{M+1}^{\mu +1}=0,}\end{array}$$

(29)

for $m=1,\dots ,M$ , $\mu =0,\dots ,K-1$ . For each $\mu$ , this system has $M+3$ equations and $M+2$ unknowns.

Then we define an implicit difference scheme for the parabolic subsystem on the concentrations $c_i$ , $i=1,\dots ,s-1$ in equations (25) to (27). It is a system of linear algebraic equations of the form

$$(\begin{array}{c}{q}_{i,0}^{\mu }{c}_{i,0}^{\mu +1}+q_{i,1}^{\mu }{c}_{i,1}^{\mu +1}=Q_{i,0}^{\mu }:=0,\hfill \\ \multicolumn{1}{c}{d_{i,m-1}^{\mu }{c}_{i,m-1}^{\mu +1}+v_{i,m}^{\mu }{c}_{i,m}^{\mu +1}+u_{i,m+1}^{\mu }{c}_{i,m+1}^{\mu +1}=Q_{i,m}^{\mu }:=e_{i,m}^{\mu }{c}_{i,m}^{\mu },}\\ \multicolumn{1}{c}{q_{i,M}^{\mu }{c}_{i,M}^{\mu +1}+q_{i,M+1}^{\mu }{c}_{i,M+1}^{\mu +1}=Q_{i,M+1}^{\mu }:=0,}\end{array}$$

(30)

where

$$\begin{array}{c}\kappa =\frac{\tau }{h^2},\hfill \\ \multicolumn{1}{c}{v_{i,m}^{\mu }=1+\kappa (D_{i,m+\frac{1}{2}}^{\mu }+D_{i,m-\frac{1}{2}}^{\mu }),}\\ \multicolumn{1}{c}{d_{i,m-1}^{\mu }=\kappa (-D_{i,m-\frac{1}{2}}^{\mu }+\frac{1}{4}(F_{m+1}^{\mu +1}-F_{m-1}^{\mu +1})),}\\ \multicolumn{1}{c}{u_{i,m+1}^{\mu }=\kappa (-D_{i,m+\frac{1}{2}}^{\mu }-\frac{1}{4}(F_{m+1}^{\mu +1}-F_{m-1}^{\mu +1})),}\\ \multicolumn{1}{c}{e_{i,m}^{\mu }=1+\kappa (F_{m+1}^{\mu +1}-2F_m^{\mu +1}+F_{m-1}^{\mu +1}),}\\ \multicolumn{1}{c}{q_{i,1}^{\mu }=-D_{i,0}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,0}^{\mu }-D_{s,0}^{\mu })(c_{k,1}^{\mu }-c_{k,0}^{\mu }),}\\ \multicolumn{1}{c}{q_{i,M}^{\mu }=-D_{i,M+1}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+1}^{\mu }-D_{s,M+1}^{\mu })(c_{k,M}^{\mu }-c_{k,M+1}^{\mu }),}\\ \multicolumn{1}{c}{q_{i,0}^{\mu }=D_{i,0}^{\mu },}\\ \multicolumn{1}{c}{q_{i,M+1}^{\mu }=D_{i,M+1}^{\mu },}\end{array}$$

for $m=1,\dots ,M$ , $i=1,\dots ,s-1$ and $\mu =0,\dots ,K-1$ . For each $\mu$ and $i$ , this system has $M+2$ equations and $M+2$ unknowns.

If the steps $h$ , $\tau$ are small enough that $v_{i,m}^{\mu (m)}\ne 0$ , then we make a sequence of the Gauss substitutions

$$\begin{array}{c}{v}_{i,k}^{\mu (0)}=q_{i,k}^{\mu },k=0,1,\hfill \\ \multicolumn{1}{c}{v_{i,m}^{\mu (m)}=v_{i,m}^{\mu }-d_{i,m-1}^{\mu }{(v_{i,m-1}^{\mu (m-1)})}^{-1}{v}_{i,m}^{\mu (m-1)},v_{i,m-1}^{\mu (m)}=0,v_{i,m+1}^{\mu (m)}=u_{i,m+1}^{\mu },}\\ \multicolumn{1}{c}{v_{i,M+1}^{\mu (M+1)}=q_{i,M+1}^{\mu }-q_{i,M}^{\mu }{(v_{i,M}^{\mu (M)})}^{-1}{v}_{i,M+1}^{\mu (M)},v_{i,M}^{\mu (M+1)}=0,}\\ \multicolumn{1}{c}{Q_{i,0}^{\mu (0)}=Q_{i,0}^{\mu },}\\ \multicolumn{1}{c}{Q_{i,m}^{\mu (m)}=Q_{i,m}^{\mu }-d_{i,m-1}^{\mu }{(v_{i,m-1}^{\mu (m-1)})}^{-1}{Q}_{i,m-1}^{\mu (m-1)},}\\ \multicolumn{1}{c}{Q_{i,M+1}^{\mu (M+1)}=Q_{i,M+1}^{\mu }-q_{i,M}^{\mu }{(v_{i,M}^{\mu (M)})}^{-1}{Q}_{i,M}^{\mu (M)},}\end{array}$$

(31)

with $m=1,\dots ,M$ , $i=1,\dots ,s-1$ .

Theorem 1.

(1) For all the steps $h$ , $\tau$ , equation (29) has exactly one solution $F^{\mu +1}$ for given concentrations $(c_1^{\mu },\dots ,c_{s-1}^{\mu })$ of the form

$$\begin{array}{c}{F}_{M+1}^{\mu +1}=\frac{\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k[(D_{k,\frac{1}{2}}^{\mu }-D_{s,\frac{1}{2}}^{\mu })c_{k,0}^{\mu }+\sum _{l=1}^M((2l+1)(D_{k,l+\frac{1}{2}}^{\mu }-D_{s,l+\frac{1}{2}}^{\mu })}{2(M+1)}\hfill \\ \frac{-(2l-1)(D_{k,l-\frac{1}{2}}^{\mu }-D_{s,l-\frac{1}{2}}^{\mu }))c_{k,l}^{\mu }-(1+2M)(D_{k,M+\frac{1}{2}}^{\mu }-D_{s,M+\frac{1}{2}}^{\mu })c_{k,M+1}^{\mu }]}{2(M+1)},\hfill \end{array}$$

(32)

$$F_m^{\mu +1}=F_{m+1}^{\mu +1}+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,m+\frac{1}{2}}^{\mu }-D_{s,m+\frac{1}{2}}^{\mu })(c_{k,m+1}^{\mu }-c_{k,m}^{\mu }),m=M,\dots ,0.$$

(ii) Equation (30) has exactly one solution $(c_1^{\mu +1},\dots ,c_{s-1}^{\mu +1})$ for given concentrations $(c_1^{\mu },\dots ,c_{s-1}^{\mu })$ and potential $F^{\mu +1}$ if and only if the steps $h$ , $\tau$ are small enough that $v_{i,m}^{\mu (m)}\ne 0$ , $m=0,\dots ,M+1$ , $i=1,\dots ,s-1$ . It has the formula

$$\begin{array}{c}{c}_{i,M+1}^{\mu +1}={\left(v_{i,M+1}^{\mu (M+1)}\right)}^{-1}{Q}_{i,M+1}^{\mu (M+1)},\hfill \\ c_{i,m}^{\mu +1}={\left(v_{i,m}^{\mu (m)}\right)}^{-1}(Q_{i,m}^{\mu (m)}-v_{i,m+1}^{\mu (m)}{c}_{i,m+1}^{\mu +1}),m=M,\dots ,0,\hfill \end{array}$$

(33)

for $i=1,\dots ,s-1$ .

Proof. We will first prove point (i). Let any $\mu$ be fixed and let concentrations $(c_1^{\mu },\dots ,c_{s-1}^{\mu })$ be known. Adding the first equation to the second one, the new second equation to the third one, and so on, leads to the system

$$\left(\begin{array}{ccccccccc}1\hfill & -1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & -1\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & -1\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 2\hfill & 2\hfill & 2\hfill & \cdots \hfill & 2\hfill & 2\hfill & 2\hfill & 1\hfill \end{array}\right)\left(\begin{array}{c}{F}_0^{\mu +1}\hfill \\ \multicolumn{1}{c}{F_1^{\mu +1}}\\ \multicolumn{1}{c}{F_2^{\mu +1}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F_M^{\mu +1}}\\ \multicolumn{1}{c}{F_{M+1}^{\mu +1}}\end{array}\right)=\left(\begin{array}{c}\sum _{j=0}^0{P}_j^{\mu }\hfill \\ \multicolumn{1}{c}{\sum _{j=0}^1{P}_j^{\mu }}\\ \multicolumn{1}{c}{\sum _{j=0}^2{P}_j^{\mu }}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{j=0}^M{P}_j^{\mu }}\\ \multicolumn{1}{c}{\sum _{j=0}^{M+1}{P}_j^{\mu }}\\ \multicolumn{1}{c}{0}\end{array}\right).$$

(34)

It follows from elementary calculations that $\sum _{j=0}^{M+1}{P}_j^{\mu }=0.$ Hence, the last equation but one can be omitted and, after the next eliminations, we get the equivalent system

$$\left(\begin{array}{ccccccccc}1\hfill & -1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & -1\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & -1\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & 2(M+1)\hfill \end{array}\right)\left(\begin{array}{c}{F}_0^{\mu +1}\hfill \\ \multicolumn{1}{c}{F_1^{\mu +1}}\\ \multicolumn{1}{c}{F_2^{\mu +1}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{F_M^{\mu +1}}\\ \multicolumn{1}{c}{F_{M+1}^{\mu +1}}\end{array}\right)=\left(\begin{array}{c}\sum _{j=0}^0{P}_j^{\mu }\hfill \\ \multicolumn{1}{c}{\sum _{j=0}^1{P}_j^{\mu }}\\ \multicolumn{1}{c}{\sum _{j=0}^2{P}_j^{\mu }}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{j=0}^M{P}_j^{\mu }}\\ \multicolumn{1}{c}{\sum _{l=0}^M(-1-2l)\sum _{j=0}^l{P}_j^{\mu }}\end{array}\right)$$

(35)

with the square $(M+2)\times (M+2)$ nonsingular matrix. In consequence, for all the steps $h$ , $\tau$ , equation (29) has exactly one solution $F^{\mu +1}=(F_0^{\mu +1},\dots ,F_{M+1}^{\mu +1})$ of the form of equation (32).

Now we will prove point (ii). Let any $\mu$ and $i$ be fixed and let concentrations $(c_1^{\mu },\dots ,c_{s-1}^{\mu })$ and the potential $F^{\mu +1}$ be known. We see that the matrix in equation (30) has a tridiagonal form. Assume that the steps $h$ , $\tau$ are small enough that $v_{i,m}^{\mu (m)}\ne 0$ , $m=0,\dots ,M+1$ in equation (31). In consequence, after the use of the Gauss elimination method (equation (31)), equation (30) has the equivalent form

$$\left(\begin{array}{cccccccc}{v}_{i,0}^{\mu (0)}\hfill & v_{i,1}^{\mu (0)}\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& v_{i,1}^{\mu (1)}\hfill & v_{i,2}^{\mu (1)}\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & v_{i,2}^{\mu (2)}\hfill & v_{i,3}^{\mu (2)}\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & v_{i,M}^{\mu (M)}\hfill & v_{i,M+1}^{\mu (M)}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & v_{i,M+1}^{\mu (M+1)}\hfill \end{array}\right)\left(\begin{array}{c}{c}_{i,0}^{\mu +1}\hfill \\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{c_{i,1}^{\mu +1}}\\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{c_{i,M+1}^{\mu +1}}\end{array}\right)=\left(\begin{array}{c}{Q}_{i,0}^{\mu (0)}\hfill \\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{Q_{i,1}^{\mu (1)}}\\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{Q_{i,M+1}^{\mu (M+1)}}\end{array}\right)$$

(36)

Hence, equation (30) has a unique solution $c_i^{\mu +1}=(c_{i,0}^{\mu +1},\dots ,c_{i,M+1}^{\mu +1})$ of the form of equation (33). □

Theorem 2. Assume that the initial concentrations $c_{0i}$ , $i=1,...,s$ are of such regularity that the solution $(c,F)$ of equations (25) to (27) belongs to $C^3([0,t_{\mathrm{end}}]\times \overline{\mathrm{\Omega }},{\mathbb{R}}^s)$ , $D_i\in C^2([0,1],{\mathbb{R}}_{+})$ , $i=1,...,s$ and $\sum _{k=1}^s{c}_k\ge \alpha >0$ for some $\alpha$ . Then the difference method (equations (29) and (30)) is consistent and the truncation errors $r,R=O(\tau +h)$ .

Proof. Taking into account the regularity assumptions, the following relations for the terms that generate the difference method (equations (29) and (30)) in the internal points hold:

$$\begin{array}{c}\frac{1}{h}((D_{k,m+\frac{1}{2}}^{\mu }-D_{s,m+\frac{1}{2}}^{\mu })\frac{1}{h}(c_{i,m+1}^{\mu }-c_{i,m}^{\mu })-(D_{k,m-\frac{1}{2}}^{\mu }-D_{s,m-\frac{1}{2}}^{\mu })\frac{1}{h}(c_{i,m}^{\mu }-c_{i,m-1}^{\mu }))-{\partial }_x((D_{k,m}^{\mu }-D_{s,m}^{\mu }){\partial }_x{c}_{i,m}^{\mu })=O(h),\\ \frac{1}{h^2}(F_{m+1}^{\mu +1}-2F_m^{\mu +1}+F_{m-1}^{\mu +1})-{\partial }_{\mathrm{xx}}{F}_m^{\mu }=O(\tau +h),\\ \sum _{m=0}^M\frac{1}{2}(F_m^{\mu +1}+F_{m+1}^{\mu +1})h-{\int }_{\mathrm{\Omega }}{F}^{\mu }\mathrm{d}x=O(\tau +h),\\ \frac{1}{\tau }(c_{i,m}^{\mu +1}-c_{i,m}^{\mu })-{\partial }_t{c}_{i,m}^{\mu }=O(\tau ),\\ \frac{1}{2h}(c_{i,m+1}^{\mu +1}-c_{i,m-1}^{\mu +1})\frac{1}{2h}(F_{m+1}^{\mu +1}-F_{m-1}^{\mu +1})-{\partial }_x{c}_{i,m}^{\mu }{\partial }_x{F}_m^{\mu }=O(\tau +h^2),\\ c_{i,m}^{\mu }\frac{1}{h^2}(F_{m+1}^{\mu +1}-2F_m^{\mu +1}+F_{m-1}^{\mu +1})-c_{i,m}^{\mu }{\partial }_{\mathrm{xx}}{F}_m^{\mu }=O(\tau +h),\\ \frac{1}{h}(D_{i,m+\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m+1}^{\mu +1}-c_{i,m}^{\mu +1})-D_{i,m-\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m}^{\mu +1}-c_{i,m-1}^{\mu +1}))-{\partial }_x(D_{i,m}^{\mu }{\partial }_x{c}_{i,m}^{\mu })=O(\tau +h).\end{array}$$

(37)

Analogously, in the boundary points:

$$\begin{array}{c}\frac{1}{h}(F_1^{\mu +1}-F_0^{\mu +1})-{\partial }_x{F}_0^{\mu }=O(\tau +h),\\ (D_{k,\frac{1}{2}}^{\mu }-D_{s,\frac{1}{2}}^{\mu })\frac{1}{h}(c_{k,1}^{\mu }-c_{k,0}^{\mu })-(D_{k,0}^{\mu }-D_{s,0}^{\mu }){\partial }_x{c}_{k,0}^{\mu }=O(h),\\ \frac{1}{h}(F_{M+1}^{\mu +1}-F_M^{\mu +1})-{\partial }_x{F}_{M+1}^{\mu }=O(\tau +h),\\ (D_{k,M+\frac{1}{2}}^{\mu }-D_{s,M+\frac{1}{2}}^{\mu })\frac{1}{h}(c_{k,M+1}^{\mu }-c_{k,M}^{\mu })-(D_{k,M+1}^{\mu }-D_{s,M+1}^{\mu }){\partial }_x{c}_{k,M+1}^{\mu }=O(h),\\ -D_{i,0}^{\mu }\frac{1}{h}(c_{i,1}^{\mu +1}-c_{i,0}^{\mu +1})+c_{i,1}^{\mu +1}\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,0}^{\mu }-D_{s,0}^{\mu })\frac{1}{h}(c_{k,1}^{\mu }-c_{k,0}^{\mu })+D_{i,0}^{\mu }{\partial }_x{c}_{i,0}^{\mu }-c_{i,0}^{\mu }\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,0}^{\mu }-D_{s,0}^{\mu }){\partial }_x{c}_{k,0}^{\mu }=O(\tau +h),\\ -D_{i,M+1}^{\mu }\frac{1}{h}(c_{i,M+1}^{\mu +1}-c_{i,M}^{\mu +1})+c_{i,M}^{\mu +1}\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+1}^{\mu }-D_{s,M+1}^{\mu })\frac{1}{h}(c_{k,M+1}^{\mu }-c_{k,M}^{\mu })+D_{i,M+1}^{\mu }{\partial }_x{c}_{i,M+1}^{\mu }\\ -c_{i,M+1}^{\mu }\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+1}^{\mu }-D_{s,M+1}^{\mu }){\partial }_x{c}_{k,M+1}^{\mu }=O(\tau +h).\end{array}$$

(38)

We will justify the most difficult last formula in equation (37) only. We will use the Taylor formula for the concentrations $c_i$ , the mole fractions $N_i$ and the diffusion coefficients $D_i$ . Note that

$$\begin{array}{c}{\partial }_x(D_{i,m}^{\mu }{\partial }_x{c}_{i,m}^{\mu })=D_{i,m}^{\mu }{\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu }+{\partial }_x{D}_{i,m}^{\mu }{\partial }_x{c}_{i,m}^{\mu }\hfill \\ =D_{i,m}^{\mu }{\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu }+({D\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }{\partial }_x{c}_{i,m}^{\mu },\hfill \end{array}$$

where

$${\partial }_x{N}_{i,m}^{\mu }=\frac{{\partial }_x{c}_{i,m}^{\mu }\sum _{k=1}^s{c}_{k,m}^{\mu }-c_{i,m}^{\mu }\sum _{k=1}^s{\partial }_x{c}_{k,m}^{\mu }}{{(\sum _{k=1}^s{c}_{k,m}^{\mu })}^2}.$$

Moreover,

$$\begin{array}{c}{D}_{i,m+\frac{1}{2}}^{\mu }=D_i\left(\frac{N_{i,m}^{\mu }+N_{i,m+1}^{\mu }}{2}\right)=D_i\left(\frac{N_{i,m}^{\mu }+N_{i,m}^{\mu }+{\partial }_x{N}_i(P_{i,m}^{\mu })h}{2}\right)\hfill \\ =D_i(N_{i,m}^{\mu }+\frac{1}{2}{\partial }_x{N}_i(P_{i,m}^{\mu })h)=D_i\left(N_{i,m}^{\mu }\right)+\frac{1}{2}{D\prime }_i\left(Q_{i,m}^{\mu }\right){\partial }_x{N}_i\left(P_{i,m}^{\mu }\right)h\hfill \\ =D_{i,m}^{\mu }+O(h)\hfill \end{array}$$

for some intermediate points $P_{i,m}^{\mu }$ , $Q_{i,m}^{\mu }$ . In the same way, we get

$$D_{i,m-\frac{1}{2}}^{\mu }=D_{i,m}^{\mu }+O(h).$$

In consequence,

$$\frac{1}{2}(D_{i,m+\frac{1}{2}}^{\mu }+D_{i,m-\frac{1}{2}}^{\mu })=D_{i,m}^{\mu }+O(h).$$

(39)

Making similar calculations as before, we obtain

$$\begin{array}{c}{D}_{i,m+\frac{1}{2}}^{\mu }=D_i(N_{i,m}^{\mu }+\frac{1}{2}{\partial }_x{N}_{i,m}^{\mu }h+\frac{1}{4}{\partial }_{\mathrm{xx}}{N}_i\left(P_{1,i,m}^{\mu }\right)h^2)\hfill \\ =D_i\left(N_{i,m}^{\mu }\right)+{D\prime }_i\left(N_{i,m}^{\mu }\right)(\frac{1}{2}{\partial }_x{N}_{i,m}^{\mu }h+\frac{1}{4}{\partial }_{\mathrm{xx}}{N}_i\left(P_{1,i,m}^{\mu }\right)h^2)+{D\prime }_i\left(Q_{1,i,m}^{\mu }\right){(\frac{1}{2}{\partial }_x{N}_{i,m}^{\mu }h+\frac{1}{4}{\partial }_{\mathrm{xx}}{N}_i\left(P_{1,i,m}^{\mu }\right)h^2)}^2\hfill \\ =D_{i,m}^{\mu }+\frac{1}{2}({D\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }h+O(h^2)\hfill \end{array}$$

for some intermediate points $P_{1,i,m}^{\mu }$ , $Q_{1,i,m}^{\mu }$ . Similarly,

$$D_{i,m-\frac{1}{2}}^{\mu }=D_{i,m}^{\mu }-\frac{1}{2}(D{\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }h+O(h^2).$$

Hence,

$$\frac{1}{h}(D_{i,m+\frac{1}{2}}^{\mu }-D_{i,m-\frac{1}{2}}^{\mu })=(D{\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }+O(h).$$

(40)

Equations (39) and (40) imply

$$\begin{array}{c}\frac{1}{h}(D_{i,m+\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m+1}^{\mu +1}-c_{i,m}^{\mu +1})-D_{i,m-\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,m}^{\mu +1}-c_{i,m-1}^{\mu +1}))\hfill \\ =\frac{1}{h^2}(D_{i,m+\frac{1}{2}}^{\mu }{c}_{i,m+1}^{\mu +1}-(D_{i,m+\frac{1}{2}}^{\mu }+D_{i,m-\frac{1}{2}}^{\mu })c_{i,m}^{\mu +1}+D_{i,m-\frac{1}{2}}^{\mu }{c}_{i,m-1}^{\mu +1})\hfill \\ =\frac{1}{h^2}(D_{i,m+\frac{1}{2}}^{\mu }(c_{i,m}^{\mu +1}+{\partial }_x{c}_{i,m}^{\mu +1}h+\frac{1}{2}{\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu +1}{h}^2+O(h^3))-(D_{i,m+\frac{1}{2}}^{\mu }+D_{i,m-\frac{1}{2}}^{\mu })c_{i,m}^{\mu +1}\hfill \\ +D_{i,m-\frac{1}{2}}^{\mu }(c_{i,m}^{\mu +1}-{\partial }_x{c}_{i,m}^{\mu +1}h+\frac{1}{2}{\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu +1}{h}^2+O(h^3)))\hfill \\ =\frac{1}{2}(D_{i,m+\frac{1}{2}}^{\mu }+D_{i,m-\frac{1}{2}}^{\mu }){\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu +1}+\frac{1}{h}(D_{i,m+\frac{1}{2}}^{\mu }-D_{i,m-\frac{1}{2}}^{\mu }){\partial }_x{c}_{i,m}^{\mu +1}+O(h)\hfill \\ =(D_{i,m}^{\mu }+O(h))({\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu }+O(\tau ))+(({D\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }+O(h))({\partial }_x{c}_{i,m}^{\mu }+O(\tau ))+O(h)\hfill \\ =D_{i,m}^{\mu }{\partial }_{\mathrm{xx}}{c}_{i,m}^{\mu }+({D\prime }_i{)}_m^{\mu }{\partial }_x{N}_{i,m}^{\mu }{\partial }_x{c}_{i,m}^{\mu }+O(\tau +h).\hfill \end{array}$$

It follows immediately from the form of equations (25) to (27), (29), (30), (37) and (38) that our difference method is consistent and the truncation errors $r,R=O(\tau +h)$ . □

Remark 3. The difference method preserves the Vegard rule (equation (13)) by equation (24). Unfortunately, we have not yet been able to prove the convergence and stability formally. In our previous paper [24], where the diffusion coefficients were constant, we have proved that the convergence and stability with respect to the concentrations were equivalent to the convergence and stability with respect to the potential, in the suitable difference maximum norms. We do not know how to use that technique in the case of the nonlinear diffusion coefficients considered in this paper.

3.2. Case n = 2

Let $\mathrm{\Omega }=(-\mathrm{\Lambda },\mathrm{\Lambda })\times (-\mathrm{\Lambda },\mathrm{\Lambda })$ . We assume that each diffusion coefficient $D_i$ depends on the mole fraction $N_i=c_i/(\sum _{k=1}^s{c}_k)$ of the $i$ th component of the mixture and the physical system is closed, i.e. $j_i(t,x_1,x_2)\equiv 0$ for $(t,x_1,x_2)\in [0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega }$ , $i=1,\dots ,s$ . We will write, for simplicity, $D_i(N_i)$ , where $D_i:[0,1]\to {\mathbb{R}}_{+}$ . Reasoning similarly as in the case $n=1$ , the initial–boundary differential problem (equations (14) to (16)) takes a form

$$(\begin{array}{c}{\partial }_t{c}_i={\partial }_{x_1}(D_i(N_i){\partial }_{x_1}{c}_i)+{\partial }_{x_2}(D_i(N_i){\partial }_{x_2}{c}_i)+{\partial }_{x_1}{c}_i{\partial }_{x_1}F+{\partial }_{x_2}{c}_i{\partial }_{x_2}F+c_i({\partial }_{x_1{x}_1}F+{\partial }_{x_2{x}_2}F)\hfill \\ \multicolumn{1}{c}{\mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },}\\ \multicolumn{1}{c}{-{\partial }_{x_1{x}_1}F-{\partial }_{x_2{x}_2}F=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k({\partial }_{x_1}\left((D_k(N_k)-D_s(N_s)){\partial }_{x_1}{c}_k\right)+{\partial }_{x_2}\left((D_k(N_k)-D_s(N_s)){\partial }_{x_2}{c}_k\right))}\\ \multicolumn{1}{c}{\mathrm{on}[0,t_{\mathrm{end}}]\times \mathrm{\Omega },}\\ \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}F\mathrm{d}x=0\mathrm{on}[0,t_{\mathrm{end}}],}\end{array}$$

(41)

$$c_i(0,x)=c_{0i}(x)\mathrm{on}\mathrm{\Omega },$$

(42)

$$(\begin{array}{cc}-D_i{\partial }_{x_j}{c}_i+c_i\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_k(D_k)-D_s(N_s)){\partial }_{x_j}{c}_k=0\hfill & \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \\ \multicolumn{1}{c}{{\partial }_{x_j}F=-\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_k(N_k)-D_s(N_s)){\partial }_{x_j}{c}_k}& \mathrm{on}[0,t_{\mathrm{end}}]\times \partial \mathrm{\Omega },\hfill \end{array}$$

(43)

with $j=1,2$ , $i=1,\dots ,s-1$ . We calculate the concentration $c_s$ from the formula

$$c_s=\frac{1}{{\mathrm{\Omega }}_s}(1-\sum _{i=1}^{s-1}{\mathrm{\Omega }}_i{c}_i).$$

We assume, additionally, that

$$\sum _{i=1}^{s-1}{\mathrm{\Omega }}_i{c}_{0i}(x)\le 1\mathrm{on}\mathrm{\Omega }.$$

(44)

We define a mesh on $\mathrm{\Omega }$ in a similar way as in the case $n=1$ . Let $h=2\mathrm{\Lambda }/(M+1)$ and $\tau =t_{\mathrm{end}}/K$ stand for the space and time steps of the mesh, respectively, where $M,K\in \mathbb{N}$ are given. Define nodal points $t^{\mu }$ , $x_{1l}$ , $x_{2m}$ , $x_{1,l+1/2}$ , $x_{2,m+1/2}$ as follows:

$$\begin{array}{c}{t}^{\mu }=\mu \tau ,\mu =0,1,\dots ,K,\hfill \\ x_{1l}=-\mathrm{\Lambda }+\mathrm{lh},x_{2m}=-\mathrm{\Lambda }+\mathrm{mh},l,m=0,1,\dots ,M+1,\hfill \\ x_{1,l+\frac{1}{2}}=-\mathrm{\Lambda }+(l+\frac{1}{2})h,x_{2,m+\frac{1}{2}}=-\mathrm{\Lambda }+(m+\frac{1}{2})h,l,m=0,1,\dots ,M.\hfill \end{array}$$

The terms ${\partial }_{x_1}(D_i(N_i){\partial }_{x_1}{c}_i)$ , ${\partial }_{x_2}(D_i(N_i){\partial }_{x_2}{c}_i)$ in equation (41) will be approximated at a point $(t^{\mu },x_{1l},x_{2m})$ with the use of a modification of the method given in [43], namely by the difference quotients, respectively,

$$\begin{array}{c}\frac{1}{h}(D_{i,l+\frac{1}{2},m}^{\mu }\frac{1}{h}(c_{i,l+1,m}^{\mu +1}-c_{i,l,m}^{\mu +1})-D_{i,l-\frac{1}{2},m}^{\mu }\frac{1}{h}(c_{i,l,m}^{\mu +1}-c_{i,l-1,m}^{\mu +1})),\hfill \\ \frac{1}{h}(D_{i,l,m+\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,l,m+1}^{\mu +1}-c_{i,l,m}^{\mu +1})-D_{i,l,m-\frac{1}{2}}^{\mu }\frac{1}{h}(c_{i,l,m}^{\mu +1}-c_{i,l,m-1}^{\mu +1})),\hfill \end{array}$$

where

$$\begin{array}{c}{D}_{i,l-\frac{1}{2},m}^{\mu }=D_i\left(\frac{1}{2}(N_{i,l-1,m}^{\mu }+N_{i,l,m}^{\mu })\right),\hfill \\ D_{i,l+\frac{1}{2},m}^{\mu }=D_i\left(\frac{1}{2}(N_{i,l,m}^{\mu }+N_{i,l+1,m}^{\mu })\right),\hfill \\ D_{i,l,m-\frac{1}{2}}^{\mu }=D_i\left(\frac{1}{2}(N_{i,l,m-1}^{\mu }+N_{i,l,m}^{\mu })\right),\hfill \\ D_{i,l,m+\frac{1}{2}}^{\mu }=D_i\left(\frac{1}{2}(N_{i,l,m}^{\mu }+N_{i,l,m+1}^{\mu })\right),\hfill \end{array}$$

for $l,m=1,\dots ,M$ . The terms ${\partial }_{x_1}((D_k(N_k)-D_s(N_s)){\partial }_{x_1}{c}_k)$ , ${\partial }_{x_2}((D_k(N_k)-D_s(N_s)){\partial }_{x_2}{c}_k)$ will be approximated analogously, but the concentrations will be taken at points $t^{\mu }$ instead of $t^{\mu +1}$ . To approximate the remaining derivatives, the central, forward and back difference quotients will be used. Moreover, the terms $D_i(N_i)$ in the first $(s-1)$ equations in equation (43) will be approximated by

$$D_{i,l,m}^{\mu }=D_i(N_{i,l,m}^{\mu })$$

for $l=0,M+1$ , $m=1,\dots ,M$ and $m=0,M+1$ , $l=1,\dots ,M$ , and, similarly, $D_k(N_k)$ , $D_s(N_s)$ . But the terms $D_k(N_k)$ , $D_s(N_s)$ in the last equation in equation (43) will be approximated by $D_{k,l,1/2}^{\mu }$ , $D_{s,l,\frac{1}{2}}^{\mu }$ , $D_{k,1/2,m}^{\mu }$ , $D_{s,\frac{1}{2},m}^{\mu }$ , $D_{k,l,M+1/2}^{\mu }$ , $D_{s,l,M+1/2}^{\mu }$ , $D_{k,M+1/2,m}^{\mu }$ , $D_{s,M+1/2,m}^{\mu }$ for $l,m=1,\dots ,M$ . We will find an approximate solution of equations (41) to (43) using a splitting linearized implicit difference method.

We define an implicit difference scheme for the elliptic subsystem on the potential $F$ in equations (41) to (43). It is a system of linear algebraic equations of the form

$$(\begin{array}{c}{F}_{0,0}^{\mu +1}-\frac{1}{2}{F}_{0,1}^{\mu +1}-\frac{1}{2}{F}_{1,0}^{\mu +1}=P_{0,0}^{\mu }:=0,\hfill \\ \multicolumn{1}{c}{-F_{l,0}^{\mu +1}+F_{l,1}^{\mu +1}=P_{l,0}^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,l,\frac{1}{2}}^{\mu }-D_{s,l,\frac{1}{2}}^{\mu })(c_{k,l,0}^{\mu }-c_{k,l,1}^{\mu }),}\\ \multicolumn{1}{c}{F_{M+1,0}^{\mu +1}-\frac{1}{2}{F}_{M+1,1}^{\mu +1}-\frac{1}{2}{F}_{M,0}^{\mu +1}=P_{M+1,0}^{\mu }:=0,}\\ \multicolumn{1}{c}{-F_{0,m}^{\mu +1}+F_{1,m}^{\mu +1}=P_{0,m}^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,\frac{1}{2},m}^{\mu }-D_{s,\frac{1}{2},m}^{\mu })(c_{k,0,m}^{\mu }-c_{k,1,m}^{\mu }),}\\ \multicolumn{1}{c}{-4F_{l,m}^{\mu +1}+F_{l+1,m}^{\mu +1}+F_{l-1,m}^{\mu +1}+F_{l,m+1}^{\mu +1}+F_{l,m-1}^{\mu +1}=P_{l,m}^{\mu }}\\ \multicolumn{1}{c}{:=-\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k((D_{k,l+\frac{1}{2},m}^{\mu }-D_{s,l+\frac{1}{2},m}^{\mu })(c_{k,l+1,m}^{\mu }-c_{k,l,m}^{\mu })-(D_{k,l-\frac{1}{2},m}^{\mu }-D_{s,l-\frac{1}{2},m}^{\mu })(c_{k,l,m}^{\mu }-c_{k,l-1,m}^{\mu }))}\\ \multicolumn{1}{c}{-\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k((D_{k,l,m+\frac{1}{2}}^{\mu }-D_{s,l,m+\frac{1}{2}}^{\mu })(c_{k,l,m+1}^{\mu }-c_{k,l,m}^{\mu })-(D_{k,l,m-\frac{1}{2}}^{\mu }-D_{s,l,m-\frac{1}{2}}^{\mu })(c_{k,l,m}^{\mu }-c_{k,l,m-1}^{\mu })),}\\ \multicolumn{1}{c}{-F_{M+1,m}^{\mu +1}+F_{M,m}^{\mu +1}=P_{M+1,m}^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+\frac{1}{2},m}^{\mu }-D_{s,M+\frac{1}{2},m}^{\mu })(c_{k,M+1,m}^{\mu }-c_{k,M,m}^{\mu }),}\\ \multicolumn{1}{c}{F_{0,M+1}^{\mu +1}-\frac{1}{2}{F}_{0,M}^{\mu +1}-\frac{1}{2}{F}_{1,M+1}^{\mu +1}=P_{0,M+1}^{\mu }:=0,}\\ \multicolumn{1}{c}{-F_{l,M+1}^{\mu +1}+F_{l,M}^{\mu +1}=P_{l,M+1}^{\mu }:=\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,l,M+\frac{1}{2}}^{\mu }-D_{s,l,M+\frac{1}{2}}^{\mu })(c_{k,l,M+1}^{\mu }-c_{k,l,M}^{\mu }),}\\ \multicolumn{1}{c}{F_{M+1,M+1}^{\mu +1}-\frac{1}{2}{F}_{M+1,M}^{\mu +1}-\frac{1}{2}{F}_{M,M+1}^{\mu +1}=P_{M+1,M+1}^{\mu }:=0,}\\ \multicolumn{1}{c}{\frac{1}{4}(F_{0,0}^{\mu +1}+F_{0,M+1}^{\mu +1}+F_{M+1,0}^{\mu +1}+F_{M+1,M+1}^{\mu +1})}\\ \multicolumn{1}{c}{+\frac{1}{2}\sum _{k=1}^M(F_{0,k}^{\mu +1}+F_{M+1,k}^{\mu +1}+F_{k,0}^{\mu +1}+F_{k,M+1}^{\mu +1})+\sum _{j=1}^M\sum _{k=1}^M{F}_{j,k}^{\mu +1}=P_{M+2,0}^{\mu }:=0,}\end{array}$$

(45)

for $l,m=1,\dots ,M$ , $\mu =0,\dots ,K-1$ . For each $\mu$ , this system has $M^2+4M+5$ equations and $M^2+4M+4$ unknowns.

Then we define an implicit difference scheme for the parabolic subsystem on the concentrations $c_i$ , $i=1,\dots ,s-1$ in equations (41) to (43). It is a system of linear algebraic equations of the form

$$(\begin{array}{c}{c}_{i,0,0}^{\mu +1}-\frac{1}{2}{c}_{i,0,1}^{\mu +1}-\frac{1}{2}{c}_{i,1,0}^{\mu +1}=Q_{i,0,0}^{\mu }:=0,\hfill \\ \multicolumn{1}{c}{s_{i,l,0}^{\mu }{c}_{i,l,0}^{\mu +1}+s_{i,l,1}^{\mu }{c}_{i,l,1}^{\mu +1}=Q_{i,l,0}^{\mu }:=0,}\\ \multicolumn{1}{c}{c_{i,M+1,0}^{\mu +1}-\frac{1}{2}{c}_{i,M+1,1}^{\mu +1}-\frac{1}{2}{c}_{i,M,0}^{\mu +1}=Q_{M+1,0}^{\mu }:=0,}\\ \multicolumn{1}{c}{q_{i,0,m}^{\mu }{c}_{i,0,m}^{\mu +1}+q_{i,1,m}^{\mu }{c}_{i,1,m}^{\mu +1}=Q_{i,0,m}^{\mu }:=0,}\\ \multicolumn{1}{c}{v_{i,l,m}^{\mu }{c}_{i,l,m}^{\mu +1}+r_{i,l+1,m}^{\mu }{c}_{i,l+1,m}^{\mu +1}+p_{i,l-1,m}^{\mu }{c}_{i,l-1,m}^{\mu +1}+u_{i,l,m+1}^{\mu }{c}_{i,l,m+1}^{\mu +1}+d_{i,l,m-1}^{\mu }{c}_{i,l,m-1}^{\mu +1}=Q_{i,l,m}^{\mu }:=e_{i,l,m}^{\mu }{c}_{i,l,m}^{\mu },}\\ \multicolumn{1}{c}{q_{i,M+1,m}^{\mu }{c}_{i,M+1,m}^{\mu +1}+q_{i,M,m}^{\mu }{c}_{i,M,m}^{\mu +1}=Q_{i,M+1,m}^{\mu }:=0,}\\ \multicolumn{1}{c}{c_{i,0,M+1}^{\mu +1}-\frac{1}{2}{c}_{i,0,M}^{\mu +1}-\frac{1}{2}{c}_{i,1,M+1}^{\mu +1}=Q_{i,0,M+1}^{\mu }:=0,}\\ \multicolumn{1}{c}{\begin{array}{c}{s}_{i,l,M+1}^{\mu }{c}_{i,l,M+1}^{\mu +1}+s_{i,l,M}^{\mu }{c}_{i,l,M}^{\mu +1}=Q_{i,l,M+1}^{\mu }:=0,\hfill \\ \multicolumn{1}{c}{c_{i,M+1,M+1}^{\mu +1}-\frac{1}{2}{c}_{i,M+1,M}^{\mu +1}-\frac{1}{2}{c}_{i,M,M+1}^{\mu +1}=Q_{i,M+1,M+1}^{\mu }:=0,}\end{array}}\end{array}$$

(46)

where

$$\begin{array}{c}\kappa =\frac{\tau }{h^2},\hfill \\ \multicolumn{1}{c}{v_{i,l,m}^{\mu }=1+\kappa (D_{i,l+\frac{1}{2},m}^{\mu }+D_{i,l-\frac{1}{2},m}^{\mu }+D_{i,l,m+\frac{1}{2}}^{\mu }+D_{i,l,m-\frac{1}{2}}^{\mu }),}\\ \multicolumn{1}{c}{p_{i,l-1,m}^{\mu }=\kappa (-D_{i,l-\frac{1}{2},m}^{\mu }+\frac{1}{4}(F_{l+1,m}^{\mu +1}-F_{l-1,m}^{\mu +1})),}\\ \multicolumn{1}{c}{r_{i,l+1,m}^{\mu }=\kappa (-D_{i,l+\frac{1}{2},m}^{\mu }-\frac{1}{4}(F_{l+1,m}^{\mu +1}-F_{l-1,m}^{\mu +1})),}\\ \multicolumn{1}{c}{d_{i,l,m-1}^{\mu }=\kappa (-D_{i,l,m-\frac{1}{2}}^{\mu }+\frac{1}{4}(F_{l,m+1}^{\mu +1}-F_{l,m-1}^{\mu +1})),}\\ \multicolumn{1}{c}{u_{i,l,m+1}^{\mu }=\kappa (-D_{i,l,m+\frac{1}{2}}^{\mu }-\frac{1}{4}\kappa (F_{l,m+1}^{\mu +1}-F_{l,m-1}^{\mu +1})),}\\ \multicolumn{1}{c}{e_{i,l,m}^{\mu }=1+\kappa (F_{l+1,m}^{\mu +1}+F_{l-1,m}^{\mu +1}-4F_{l,m}^{\mu +1}+F_{l,m+1}^{\mu +1}+F_{l,m-1}^{\mu +1}),}\\ \multicolumn{1}{c}{s_{i,l,0}^{\mu }=D_{i,l,0}^{\mu },}\\ \multicolumn{1}{c}{s_{i,l,M+1}^{\mu }=D_{i,l,M+1}^{\mu },}\\ \multicolumn{1}{c}{q_{i,0,m}^{\mu }=D_{i,0,m}^{\mu },}\\ \multicolumn{1}{c}{q_{i,M+1,m}^{\mu }=D_{i,M+1,m}^{\mu },}\\ \multicolumn{1}{c}{s_{i,l,1}^{\mu }=-D_{i,l,0}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,l,0}^{\mu }-D_{s,l,0}^{\mu })(c_{k,l,1}^{\mu }-c_{k,l,0}^{\mu }),}\\ \multicolumn{1}{c}{s_{i,l,M}^{\mu }=-D_{i,l,M+1}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,l,M+1}^{\mu }-D_{s,l,M+1}^{\mu })(c_{k,l,M}^{\mu }-c_{k,l,M+1}^{\mu }),}\\ \multicolumn{1}{c}{q_{i,1,m}^{\mu }=-D_{i,0,m}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,0,m}^{\mu }-D_{s,0,m}^{\mu })(c_{k,1,m}^{\mu }-c_{k,0,m}^{\mu }),}\\ \multicolumn{1}{c}{q_{i,M,m}^{\mu }=-D_{i,M+1,m}^{\mu }+\sum _{k=1}^{s-1}{\mathrm{\Omega }}_k(D_{k,M+1,m}^{\mu }-D_{s,M+1}^{\mu })(c_{k,M,m}^{\mu }-c_{k,M+1,m}^{\mu }),}\end{array}$$

for $l,m=1,\dots ,M$ , $i=1,\dots ,s-1$ and $\mu =0,\dots ,K-1$ . For each $\mu$ and $i$ , this system has $M^2+4M+4$ equations and $M^2+4M+4$ unknowns.

Remark 4. Similar theorems and remarks to Theorems 1 and 2 and Remark 3 are true. The proofs are analogous.

4. Experiment and simulations

In this section, we explicitly show the results, experiment and simulations. We start with the diffusion multiple experiment and perform the computations of the concentration maps and profiles and the drift velocity.

4.1. Diffusion multiple

To validate the model and numerical simulations, a diffusion multiple was prepared from pure Fe, Co and Ni, each with a purity of at least 99.95%, provided by Alfa Aesar. In subsequent stages, the materials were polished down to $1$ $\mathrm{\mu }\mathrm{m}$ according to the standard Struers procedure, cleaned and dried. Every annealing was carried out at a temperature of $1373$ K under a pressure below ${10}^{-4}$ Pa. The samples were prepared following the procedure shown in Figure 1. First the Fe|Ni diffusion couple was prepared by annealing for $30$ min and quenching in a cold Al tube, still in a vacuum (Figure 1(a)). During annealing, the metals were placed in a molybdenum holder, thanks to which a good join between the metals was provided. The obtained couple was cut and connected with a Co plate, annealed for $30$ min and quenched (Figure 1(b)). The obtained multiple was cut into $2$ mm slices (Figure 1(c)). An arbitrarily chosen slice was subjected to $199.5$ h annealing. Finally, the diffusion multiple was thinned, polished and subjected to an examination using a FEI Versa 3D scanning electron microscope equipped with an EDAX Apollo XP Silicon Drift detector. The resulting multiple is shown in Figure 1(d). The inset presents a fragment chosen for the measurements and the simulations.

Figure 1.

Steps to prepare the Fe-Co-Ni multiple sample: (a) Fe|Ni couple, (b) multiple, (c) cutting, (d) sample; the inset shows an examined area.

The Fe, Co and Ni concentration maps were measured with a resolution of $1024\times 926$ points.

4.2. Calculations

We consider a ternary Fe-Co-Ni mixture and assume the following data:

$$\begin{array}{c}\mathrm{components}:s=3,1\equiv \mathrm{Fe},2\equiv \mathrm{Co},3\equiv \mathrm{Ni},\mathrm{\Lambda }=0.042\mathrm{cm},\hfill \\ ({\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2,{\mathrm{\Omega }}_3)=(7.09,6.67,6.59)\mathrm{c}{\mathrm{m}}^3/\mathrm{mol},\hfill \\ (c_{01},c_{02},c_{03})=(0.0,0.149925,0.0)\mathrm{mol}/\mathrm{c}{\mathrm{m}}^3\mathrm{in}[-\mathrm{\Lambda },\mathrm{\Lambda }]\times (\mathrm{\Lambda }/3,\mathrm{\Lambda }],\hfill \\ (c_{01},c_{02},c_{03})=(0.141044,0.0,0.0)\mathrm{mol}/\mathrm{c}{\mathrm{m}}^3\mathrm{in}[-\mathrm{\Lambda },0]\times [-\mathrm{\Lambda },\mathrm{\Lambda }/3],\hfill \\ (c_{01},c_{02},c_{03})=(0.0,0.0,0.151745)\mathrm{mol}/\mathrm{c}{\mathrm{m}}^3\mathrm{in}(0,\mathrm{\Lambda }]\times [-\mathrm{\Lambda },\mathrm{\Lambda }/3],\hfill \\ t_{\mathrm{end}}=200\mathrm{h}.\hfill \end{array}$$

For the mixture, the diffusion coefficients $D_i$ in our differential problems in Section 3 depend on the mole fractions $N_i=\frac{c_i}{\sum _{k=1}^s{c}_k}$ . To find them, we use the method by Wierzba et al. [44]. For the temperature of the experiment ( $1373$ K), they can be defined as the following functions:

$$\begin{array}{c}{D}_1(N)={10}^{6.361N^2-6.4106N-10.3428}\mathrm{c}{\mathrm{m}}^2/\mathrm{s},\hfill \\ D_2(N)={10}^{-2.1363N^2+2.2129N-10.8572}\mathrm{c}{\mathrm{m}}^2/\mathrm{s},\hfill \\ D_3(N)={10}^{0.7582N^2+0.7734N-10.9198}\mathrm{c}{\mathrm{m}}^2/\mathrm{s},\hfill \end{array}$$

where $N\in [0,1]$ . Figure 2 shows the adequate graphs.

Figure 2.

Diffusion coefficients, cm²/s, for the components in the ternary Fe-Co-Ni mixture.

The one-dimensional concentration profiles were calculated using the difference method presented in case $n=1$ in Section 3. The method presented in case $n=2$ in Section 3 was used to compute two-dimensional maps and to generate the profiles.

4.2.1. Maps and profiles: comparison with experiment

Experimental concentration maps, processed with the Savitzky–Golay filter to reduce noise, are shown in Figure 3(a). The scanning electron micrographs reveal the Frenkel porosity in iron. For binary couples chosen from the multiple, the positions of the Matano’s planes were calculated using the procedure of Leszczyński et al. [45]. They were used to simulate the interdiffusion in the multiple. The resulting concentrations maps are shown in Figure 3(b). Together with the maps, we show the simulated fluxes in Figure 3(c).

Figure 3.

Experimental (a) and simulated (b) concentration maps, and the fluxes (c) in the Fe-Co-Ni multiple after annealing at 1373 K for 200 h.

For better presentation of the effectiveness of the present model, the concentration profiles for various one-dimensional cross-sections were simulated (Figures 4 to 8). In the figures, we compare the experimental with simulated profiles. In simulations, the one- and two-dimensional difference methods were used (Section 3). Generally, a good qualitative agreement of all sets of data is seen.

Figure 4.

Scanning electron micrographs and concentration profiles at the cross-sections: Fe|Ni, Co|Ni, Co|Fe. The Matano plane position is indicated.

Figure 5.

Concentration profiles at the one-dimensional Co|Fe|Ni cross-section (shown on the left) of the multiple.

Figure 6.

Concentration profiles at the one-dimension Fe|Ni|Co cross-section (shown on the left) of the multiple.

Figure 7.

Concentration profiles at the one-dimension Co|Ni cross-section (shown on the left) of the multiple.

Figure 8.

Concentration profiles at the one-dimension Fe|Co cross-section (shown on the left) of the multiple.

In Figure 4, the results for the three pseudo one-dimension couples are shown. We present two sets of simulated data, computed using the one- and two-dimensional methods.

The computed and experimental profiles agree perfectly with experimental data only for the Co|Ni and Co|Fe couples. In the case of the Fe|Ni couple, the data simulated with the one- and two-dimensional models overlap but the experimental results clearly differ. The discrepancy can be attributed to a big fraction of the voids in iron next to the contact plane.

In Figures 5 and 6, we compare the results for the Co|Fe|Ni and Fe|Ni|Co cross-sections passing the three-metal contact point. Here we compare the results from the one- and two-dimensional methods with the experimental ones. For all components, a clearly better agreement with experimental results is provided by the two-dimensional model, which is primarily seen in the case of the central components, Fe and Ni, respectively.

An advantage of the two-dimensional method over the one-dimensional method is particularly seen in Figures 7 and 8, where the concentration profiles at the one-dimensional cross-sections going through the three-metal joint point are shown. Only the two-dimensional method has been used in this case. The one-dimensional model neglects one of the multiple components.

4.2.2. Time evolution

An advantage of any numerical model of the diffusion is that it can be applied to simulate the processing. In Figures 9 to 12, we present the simulations of the time evolution and we show the results for the times $0$ h, $100$ h and $200$ h.

Figure 9.

Time evolution of the concentrations of Fe ( $c_{551}$ ), Co ( $c_{552}$ ) and Ni ( $c_{553}$ ), for 0, 100 and 200 h.

Figure 10.

Time evolution of the norm of the drift velocity $\parallel v^{55\mathrm{d}}\parallel =\parallel -\nabla F\parallel$ , for 0, 100 and 200 h.

Figure 11.

Time evolution of the divergence of the drift velocity div $(v^{56\mathrm{d}})=-\mathrm{\Delta }F$ , for 0, 100 and 200 h.

Figure 12.

Modules of difference between the experimental and simulated (the last column of Figure 9) concentrations, for 200 h.

5. Conclusions

The results of the paper are summarized as follows:

1. We presented a theoretical study of the model introduced in [23, 24] of the one- and multidimensional transport in an $s$ -component solid solution. The model was expressed as a nonlinear parabolic–elliptic system of strongly coupled differential equations with initial and nonlinear coupled boundary conditions.

2. The implicit finite difference methods in one and two dimensions were constructed, for the case of nonlinear diffusion coefficients depending on mole fractions. Let us stress that in [24] only the case of constant diffusion coefficients was studied.

3. The existence and uniqueness of solutions of the implicit difference schemes and consistency of the difference methods were proved. The Vegard rule was preserved.

4. The numerical results were compared with experimental data for a ternary Fe-Co-Ni system. A good agreement of both sets was revealed, which confirmed a force of the model and numerical methods. The comparison particularly demonstrates that our methods are convergent and stable.