Abstract

The method of lines (MOL) for diffusion equations with Neumann boundary conditions is considered. These equations are transformed by a discretization in space variables into systems of ordinary differential equations. The proposed ODEs satisfy the mass conservation law. The stability of solutions of these ODEs with respect to discrete L2 norms and discrete norms is investigated. Numerical examples confirm the parabolic behaviour of this model and very regular dynamics.

1. Introduction

Diffusion is one of several transport phenomena that occur in nature. Many physical processes in metallurgy are controlled by diffusion, for example, oxidation and sintering. There are two ways to introduce the notion of diffusion: an atomistic approach and continuum approach. The central role in the diffusion theory of multicomponent systems plays the Kirkendall effect [1], i.e., the discovery that drift in solids can be generated by diffusion. The explanation of the effect and an original method of its mathematical analysis was proposed by Darken [2]. The Kirekndall effect describes what happens when two solids diffuse into each other at different rates. This subject is treated in many textbooks on solid-state and publications [3–7]. In [8, 9] the authors deal with stability of Kirkendall planes. The peculiarities of the deformation in one-dimensional binary diffusion couple with a moving interface are described in [10].

Our ternary system arises from the Fick second law. If the overall molar density and component molar ratio are defined by
is the molar concentration of the th component in the mixture, [mole/dm3])
where this quatity is a ratio of th component atoms to averall number of atoms. Then we have
where denotes the intrinsic diffusivity of the th component and is the drift velocity. In the paper , will be denoted by , , and , respectively.

The method proposed by Darken [2], was extended by Danielewski et al. in [11, 12]. Their generalization of Darken’s method describes the interdiffusion process in bounded mixtures with constant concentrations and with variable diffusivities of the elements.

We are interested in establishing an approximation method of solutions to the ternary system mentioned above by solutions of associated systems of ordinary differential equations. These systems of ordinary differential equations are obtained by using a discretization in spatial variables.

From the abundant literature concerning the numerical method of lines (MOL) for classical PDEs we mention the monographs [13–15]. MOL for time-dependent one-dimensional systems of parabolic partial differential equations is proposed in [16]. Numerical examples for problems such as Burger’s equation or nonlinear diffusion equation with nonlinear boundary conditions were presented in [16]. A posteriori error estimates for evolution problems (general nonlinear parabolic problem with a strongly monotonic elliptic operator, a linear nonstationary convection-diffusion problem, and linear second-order hyperbolic problem) solved by the method of lines are derived in [17].

The aim of the paper is to construct a method of lines for diffusion equations in three-component system with Neumann boundary conditions. We prove stability of approximate solutions with respect to discrete norms and discrete norms.

Our research is concerned with the situation where the initial data and solutions are at least of class . This regularity is assumed throughout the paper. It is observed that even very irregular data in a short time give highly smooth solutions (concentrations of the components). Moreover, these concentrations are mixed, which means that they are inside the interval for . Therefore, we assume that the initial data are of class . This observation also justifies our assumptions that the Lipschitz constants and second derivatives of solutions are bounded by some constant independent of the discretization parameter. Under these assumptions we prove that solutions of ODEs satisfy the maximum principle and remain inside the interval . We also show convergence and stability of the approximate solutions. Our results confirm the parabolic nature of the equations in question. The existence and uniqueness of solutions of the interdiffusion problem in space was shown by the variational methods in [12]. In [18] authors constructed an example of an interdiffusion equation that does not fulfill a peculiar type of parabolicity, meant in the sense of a decreasing norm. In Section 4 we show the convergence of MOL in discrete norms, while in Section 3 we attempt to show convergence in discrete norms. It turns out that we have stability in norms for certain values of the diffusion coefficients , , , which highlights the parabolic nature of the equation.

In Section 5 we present numerical examples where we consider Darken’s trajectory besides solutions. Since the right-hand sides of ODEs are polynomial, we attempt in the appendix to estimate the radius of convergence of Taylor series solutions. It turns out that this radius is proportional to . This means that Cauchy-Kovalevskaya theorem has limited (restricted) applicability to our system. Because R-K methods are close to Taylor expansions, our observation points out that one should be careful with these methods.

2. Formulation of the Problem

We consider the system of diffusion equations:
on with the initial boundary conditions
where the diffusion coefficients , , are different and positive, say , , and for , and
is the drift velocity. These data imply .

We formulate the method of lines corresponding to (4) on a regular mesh with step such that . The difference operators , , are defined in the following way:

Denote by the discrete drift velocity . We will consider the following ODE system:
for with the initial conditions
The discrete Neumann boundary conditions
can be regarded as a convenient definition of auxiliary quantities: , , , , , .

Remark 1. Local existence and uniqueness of solutions of (8) with the above conditions follows from the general theory of ODE systems.

Suppose that we have two constants , , independent of , satisfying the inequality . This assumption is valid throughout the paper. will be the Lipschitz constant for the initial functions , , . will play the same role for , ,
We will refer to these Lipschitz-type inequalities in the definitions of and as -Lipschitz and -Lipschitz, respectively. It is assumed that approximate solutions are Lipschitz continuous with respect to spatial variable. The Lipschitz continuity of approximate solutions can be proved rigorously; however, we omit the proof for brevity. To justify our proceeding, consider solutions of ternary systems. They are bounded together with their first and second derivatives.

If the initial data belong to the equations for can be eliminated from (8). The following lemma allows us to express (8) as the system of ordinary differential equations.

Lemma 2. Suppose that is the solution of (8) with the initial condition . Assume that , , are -Lipschitz. Then .

Proof. First, we prove that for . Let . We have
for with the initial condition and the discrete Neumann boundary conditions. It follows from the general ODE theory that there exists exactly one solution .We next show that on for . Suppose first that . We claim that on for . To obtain a contradiction, we suppose that there is such that for and and there is such that . We have for . Then
We exploit (8) at the point to get the following inequality:
The expressions in the brackets are positive for small . Hence
We see at once that for . It follows from the general ODE theory [19, page 108] that for , contrary to strict initial inequalities. This argument also settles the case of or . We now turn to the case . It follows from the theorem on the continuous dependence on initial data [19, page 145] that for .

Denote by the trapezoidal sum. More precisely,
This quantity corresponds to the total mass of . The MOL system is conservative. It preserves the mass of , .

Proof. Note that for . Hence
where and . It follows from these equations that

3. Stability with respect to Discrete Norms

The stability of the method of lines with respect to discrete norm is established by the following theorem.

Theorem 4. Suppose that (1) is the unique solution of system (8) with the initial condition ,(2) is an approximate solution of (8) with the initial condition and with perturbations of the right-hand sides denoted by , , for ,(3)consider
Then we have
where are constants independent of and

Proof. Denote by , , the differences between the approximate and exact solutions of (8).The proof is based on the following observation:
The terms and are calculated similarly as . Set
We apply the summation by parts and obtain
From Lemma 3 we have . Set
Hence the sign of depends on the quadratic form :
which is negative definite for , , satisfying condition (20) of Theorem 4. We have
It follows from the inequality
that
Hence
We conclude from this differential inequality that
where
This finishes the proof.

4. Stability with respect to Discrete Norms

We define the norm in as follows:
where (with the discrete Neumann convention).

Theorem 5. Suppose that (1) is the exact solution of (8) with the initial condition ,(2) is an approximate solution of (8) with the initial condition and with perturbations of the right-hand sides denoted by , , for .Then we have
where is some constant independent of .

Proof. It is easily seen that satisfies differential equation
satisfies a similar differential equation.Suppose that we have Green functions , corresponding to the differential-difference operator:
such that
where and (see [20]). Using the above Green functions we express the differences and as follows:
where , depend on , , , , , and satisfy the estimates
where depends on , , . Hence
where
The same inequalities are satisfied by , . We get
It follows from Lemma A.1 that
where depends on , .

5. Numerical Simulations

Numerical Example 1. Set , , , and the initial distributions , are almost step functions, regularized near zero, whereas starts from a droplet:
for where stands for a characteristic function. We use a second-order Runge-Kutta method with time step . Here , , satisfy condition (20) of Theorem 4. Note that with . In Figure 1(a) we give experimental values of the solution of (8) for and with the initial condition (45) (a). Since we divide into intervals. We put the values of , , at . Note that is increasing on the interval . Furthermore is decreasing on the interval . The even function (“marker”) is increasing on the interval and decreasing on the interval . We consider the differential equations describing the marker position , where is the Darken velocity. We use a second-order Runge-Kutta method. The approximate position of the marker denoted by (Figure 2(b)) forms a trajectory in -plane shown in Figure 2(a). This object traces the evolution of concentrations of diffunding substances , .

Figure 2: Evolution of concentrations , (a) and experimental position of marker for (b).

Numerical Example 2. Set , , , and the initial distribution of is constant:
for where is the characteristic function of . Note that with . In Figure 3(a) we give initial values of the solution of (8) with the initial condition (46) for and . Starting with a constant distribution (Figure 3(a)) we give experimental values of , , for and (Figure 3(b)). It is seen that diffusion causes a perturbation for . This observation is coherent with engineers’ practice.

Figure 3: Experimental values of , , and for , .

Our numerical experiments confirm the theory and stability analysis. It is seen that the proposed numerical method is stable, mass-conservative, and dissipative. Its solutions do not leave the interval . Since our stability analysis relies on the -Lipschitz condition for , , , we show (numerically) that the Lipschitz constants for , , decrease in time. We focus on numerical Example 1. Denote these time-dependent constants by , , , respectively, and show them in Figure 4: (a) and (b). Considering the same numerical example we demonstrate convergence in Figure 5 where numerical solutions for different steps , , and are very close to each other. We take the concentrations (a) and (b).

Figure 4: Experimental values of Lipschitz constants and for , , and .

Appendix

Lemma A.1. Suppose that ,
and . Then there is a positive constant such that
for , .

Since the right-hand sides of system (8) are polynomials, we can apply the Cauchy-Kovalevskaya theorem. Suppose that solutions of (8) are given by

If we differentiate these expressions with respect to we get the obvious recursive relations for and . One can prove by induction on that
where
This means that the radius of convergence of the above series is lower bounded by , where

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the National Science Center (Poland) decision no. DEC-2011/02/A/ST8/00280.