Asymptotic Expansion Method with Respect to Small Parameter for Ternary Diffusion Models.

Ternary diffusion models lead to strongly coupled systems of PDEs. We choose the smallest diffusion coefficient as a small parameter in a power series expansion whose components fulfill relatively simple equations. Although this series is divergent, one can use its finite sums to derive feasible numerical approximations, e.g. finite difference methods (FDMs).

Introduction

The subject of divergent series dates back to the distant past and is linked with the names of outstanding mathematicians such as Euler, Poincaré, Borel, Padé, and Birkhoff. For more details and additional references we refer to articles by Tucciarone [1] and Ferraro [2]. The idea of the asymptotic series comes from Poincaré, who introduced the definition of an asymptotic expansion. Nowadays, the importance of methods which are based on asymptotic expansions of solutions in small or large parameters series has grown considerably in many branches of physical, chemical, biological and engineering sciences. We mention several examples among a wide range of applications, and the list is not exhaustive. The fundamental problem of reconstruction of physical observables, based on divergent power series expansions, has been considered in many areas of quantum physics. Dyson already showed in his pioneering work [3] that expansions in quantum electro-dynamics diverge factorially. After that, it was found that divergent perturbation expansions of quantum physics are very common [4, 5]. This observation led to a substantial amount of research works in various areas of quantum physics, for example in quantum field theory [6-8] or in quantum mechanics [9, 10]. The asymptotic analysis of singularly perturbed problems is performed in the monograph [11], based on expansions of analytical solutions in the power series with respect to a small parameter. In the light of this approach there are discussed models of populations, epidemiological problems, classical models of fluid dynamics and many others. The paper [12] deals with the method of asymptotic expansions of solutions to a singularly perturbed system of integro-differential equations in epidemiology with two small parameters. Asymptotic solutions are constructed by the Tikhonov-Vasil’eva method of boundary functions.

This work concerns a system of diffusion equations. When the number of components exceeds two, the term ‘interdiffusion’ or ‘cross-diffusion’ is often used, see [13] in material sciences and [14] in life sciences. The drift in each of the PDEs depends on the gradient of all other variables, so these models exhibit the phenomenon of cross-diffusion. This makes the PDE system strongly coupled. In addition, the coefficients depend on the unknown functions and the parabolicity of this system is conditional, namely it can be proved only if solutions remain in an admissible set. The invariance of this set is connected to mass conservation laws. Any violation of this property leads to serious analytical and computational problems. We overcome these problems by constructing conservative numerical schemes, i.e. preserving a discrete version of mass. In spite of mathematical and numerical difficulties, three dimensional interdiffusion problems are still of interest. Interdiffusion becomes an important topic in electrochemistry and biochemistry, e.g. molecular channels (nano-channels) [15].

Preliminaries

Let be a bounded domain in with a smooth boundary . The closure of is denoted by . In this manuscript, we investigate the system of diffusion equationson , where and is the drift velocity given byThe system (1) is considered with the initial conditionsand the Neumann boundary conditionswhere denotes the outward normal vector to the boundary .

Our analysis of the strongly coupled system (1) starts from the simplest case . Consider two chemical substances and their normalized densities , described by the system of differential equationsIt is clear that the assumptions on initial functions and imply , . Indeed, consider new dependent variables p, q defined as linear combinations of the functions u, v:Then we obtain the following systemSince it is seen that solves the first equation of system (6). Based on the relation we obtain from (4) partial differential equationTherefore, the implication follows easily from the maximum principle.

There are various numerical techniques for approximations of solutions of the ternary diffusion models (1)–(3). In particular, we refer to the papers [16, 17]. The authors of [16] analyze the method of lines for diffusion equations. The paper [17] deals with iterative methods for this problem. The authors of both papers study stability and convergence of the constructed method in and in the Sobolev space . In this paper, as a tool for approximation of solutions to the initial boundary value problems for systems of equations obtained by an asymptotic expansion, we use a finite-difference methodology. In the construction of difference schemes we apply general ideas presented in [18, 19].

Small Parameter Method

Suppose that initial functions in (2) satisfy the relation for . In the sequel we assume that the diffusion coefficients in the system (1) satisfy the relation . If u, v, w satisfy (1)–(3), then and the system (1) reduces to the system of two equationsIt is more convenient to use the same change of variables u, v to p, q as in (5): , For these functions we obtain the system of differential equationswhere . Observe that if then (7) is equivalent to (6). We expand the functions p and q in the series with respect to the small parameter :Once we substitute these expansions to the system (7) and compare the coefficients we get the recurrence relationfor . We represent the initial conditions for the unknown functions and as follows:where , . It is easy to see that . If we compare the right-hand sides of (10) and (8), we obtain the initial conditions for the unknown coefficients for :for .

Using a small parameter expansion of functions p and q we replace the original system of equations (7) by a separate systems of equations (9) for each pair of functions , . The complexity of these systems is rising with increasing values of k. Therefore, the expansions (8) are practically limited to several terms. According to this we are mainly interested in the asymptotic nature of this truncated expansion for , , then the possible divergence of the whole series is in general not important.

The study of asymptotic expansions was first made by Poincaré. According to Poincaré’s definition, a series is said to be asymptotic to a function , i.e.iff for every and , , we haveThis definition is equivalent to the following conditionsTaking the first limit under consideration, we haveSince , the functions p(t, x; 0), become solutions of the same initial boundary value problem which leads to the conclusion that the limit is equal to 1 for every N. The same result we obtain for the expansion of . Therefore, the expansions (8) are asymptotic if only . Even if the series do not converge, the asymptotic expansions give a good approximation of solutions if we take under consideration only few terms of these expansions with a very small value of [11, 20]. This conclusion is confirmed by our numerical analysis in one dimensional case presented in Example 1. It turns out that our method is very effective for multidimensional couples, even in the case of complicated geometries.

Numerical Analysis

This section is devoted to a numerical analysis of solutions to the problem (9), (10), obtained by the small parameter expansion method. We construct a conservative difference scheme and present the numerical behavior of ternary solid solutions based on this method.

Finite Difference Method

We discretize the system of partial differential equations (9). Consider a uniform partition of , where . Let are positive integers and Let , , stand for steps of a regular mesh on , such that and . For we put andWe convey that will represent an approximation to the exact value of the functions at the point for and .

Let , , where 1 is the jth coordinate. We employ standard and linear difference operators , , :Denote by the discrete equivalent of the gradient operator . Let then we defineWith these conventions, the finite-difference scheme will be given by the following discrete systemwith the discrete initial conditionsThe discrete Neumann boundary conditions have the formfor .

Numerical Simulations

In this section we explore the capability of the constructed method (11)–(13) to provide good approximations to the exact solution of the PDE problem (1)–(3). We present some examples of tracer and up-hill diffusion in and tracer diffusion in . Examples in one-dimensional case are similar to those presented in [16]. Therefore, one can compare the results of the method of lines obtained in [16] with the results obtained by our numerical scheme. We also provide more complicated two-dimensional cases which are absent from the literature.

Example 1

(Tracer at interface) Let . In the numerical analysis we take and the diffusion parameters . The initial distribution is a small perturbation of the Heaviside function and is a symmetry of about the ordinate axis. The solution w is initiated by a smoothed small peakwhere denotes the characteristic function of a set A, andfor . We analyze the results for the PDE system (1) with the initial distributions given by (14) with respect to the small parameter . Illustrative simulations are provided with: and . Based on formula (5) we can deriveWe consider expansions (8) for :

, one term in each expansions: , . Numerical results for the respective , , calculated by (15) are presented in Fig. 1. Note that we obtain an expected behavior of densities u, v, w regardless of the value selection.

Fig. 1

(Tracer at interface) the case : functions for different values of at time

, two terms in each expansions: . In Fig. 2 we observe the distributions of computed by (15) with different values of . Figure 3 presents the behavior of partial sums with different values of . Although oscillations of increase with decreasing values of , we obtain reasonable approximations of u, v, w by means of partial sums , respectively. Notice that these partial sums take values in the interval [0, 1] for all .

Fig. 2

(Tracer at interface) the case : functions for different values of at time

Fig. 3

(Tracer at interface) the case : partial sums for different values of at time

, three terms in each expansions: , . Figure 4 presents the behavior of , computed by (15). Numerical results of partial sums are presented in Fig. 5. Similarly as in the case , it can be observed that the oscillations of , increase. Unlike the case , the partial sums for u, v, w take values in [0, 1] for sufficiently small . If , the partial sums for u, v are not monotone. In Fig. 6 one can observe additional results for partial sums at times and with fixed . We present these results to illustrate the dynamics of approximate distributions.

Fig. 4

(Tracer at interface) the case : functions for different values of at time

Fig. 5

(Tracer at interface) the case : partial sums for different values of at time

Fig. 6

(Tracer at interface) the case : partial sums for u, v, w at time and

We conclude that the number of terms in the expansions (8) is closely related to the size of the parameter Once we use more terms in the expansions, we need to decrease the value of to get reasonable approximation of u, v, w by means of respective partial sums, i.e. taking values in the interval [0, 1] and preserving monotonicity of partial sums for u, v.

Example 2

(Up-hill diffusion) We consider the one-dimensional case and we set similarly as in Example 1. We take . We put some monotone functions as initial distributions of substances u, v and we assume a constant initial distribution w:whereResults of numerical simulations for partial sums computed with are given in Fig. 7(b) at time . Figure 7a presents initial functions .

Fig. 7

(Up-hill diffusion) a initial functions , b partial sums at time

No numerical analysis for tracer diffusion in two-dimensional case have been done so far in the literature because of the computational complexity, not to mention any efficient comparison of numerical results with real-life experiments. Due to our approach we constructed a numerical method which allows to perform computations even for multi-dimensional cases. We present two examples of numerical simulations for tracer diffusion in .

Example 3

(Tracer diffusion ) Consider the case . Set and . Initial functions (Fig. 8) are given bywhereIn Fig. 9a we observe approximate distributions of partial sums for u, v, w with and the step sizes , at time . Figure 9b presents enlargement of w.

Fig. 8

(Tracer diffusion) initial distributions: a , b enlargement of

Fig. 9

(Tracer diffusion) partial sums with at time : a for u, v, w, b enlargement of w

In the above example initial distribution of u and v is symmetrical with respect to the plane . Now we present two-dimensional example where the symmetry of initial distribution is broken.

Example 4

(Tracer diffusion) Consider the case . Set and . Initial functions (Fig. 10) are given bywhere and are the same as in Example 3. We perform simulations with parameters: the step sizes and the length of series expansions . In Fig. 11(a) the approximate distributions of partial sums for u, v, w at time are presented. Figure 11b shows enlargement of w.

Fig. 10

(Tracer diffusion) initial distributions: a , b enlargement of

Fig. 11

(Tracer diffusion) partial sums with at time : a for u, v, w, b enlargement of w

We conduct a brief analysis of marker positions in time, see [21]. By ‘marker’ we understand a small amount of substance placed in the diffusion couple, this substance moves according to locally not-balanced diffusive fluxes. Recall that diffusive fluxes are proportional to gradients of u, v, w. In the above tracer diffusion examples the marker is represented by the substance w. Referring to Example 1 the trajectory of the marker position is described by the ODE initial value problemwhere is the drift velocity. Since the diffusion coefficient is small, it can be omitted,then the solution of the Cauchy problem (16) can be approximated by the solution of the following ODE problemTo solve the problems (16) and (17) numerically, we first use a polynomial fitting for discrete functions u, v, w with respect to spatial variable. Then approximate solutions to the above ODE problems are computed by a second-order Runge–Kutta method. The results of computations with step sizes: are illustrated in Fig. 12. A mesh refinement leads to better accuracy, i.e. the solution to problem (17) gives a very good approximation of the solution to problem (16). An analogous discussion can be conducted for Examples 3 and 4.

Fig. 12

The trajectory of marker position—approximate solutions of problems (16) (solid line) and (17) (dashed line) for

Conclusions and Remarks

Before we close this paper with concluding remarks we would like to mention other applications of asymptotic expansion methods, in addition to those given in the introduction, to point out the importance of these techniques in a different branches of knowledge. Some applications can be observed in the fluid and solid mechanics theory. In [22] the asymptotic expansion method was used as a tool to derive a new two-dimensional shallow water model from the time-averaged non-dimensional Navier–Stokes equations. Asymptotic simplification procedures for linear and nonlinear wave propagation problems that contain large parameters are examined in [20]. A divergent asymptotic series found applications also in the boundary layer theory [23] and in the nonlinear wave theory [24].

We use the asymptotic expansion techniques to the ternary diffusion problem. Modeling of interdiffusion phenomena is not yet unified, and various methods differ from each other in the arbitrary choice of the reference velocity, i.e. the drift or convection velocity. In computational solid mechanics and physics it is defined as being equal to the mass average velocity [13]. In gases and fluids the drift is defined based on the volume average velocity [25]. It is still common to neglect drift (convection) and assume Fickian diffusion [26]. Such simplified approach to the mass transport was used there to a three-dimensional model of interdiffused quantum dots. The above methods do not allow considering all interdiffusion effects [26]. The use of widely accepted Onsager method [27, 28] is narrow by very limited data on transport coefficients. In material sciences and chemistry the drift definition bases on the Darken method [29, 30]. The Darken’s basic postulate is that the reference velocity in multicomponent systems coincides with the drift velocity [30]. The drift in solids is often the vacancy flux generated during an interdiffusion process caused by the difference in the diffusion coefficients. The compatibility of Darken and Onsager formalisms has been already proved [31]. Onsager’s and Fickian laws are known in life sciences, see [14].

In this paper we consider a strongly coupled system of diffusion equations with Neumann boundary conditions. This system is strongly coupled, highly nonlinear and possesses a parabolic type only in a very restricted subregion of the phase space, which is related to the mass conservation. This causes analytical and numerical difficulties. To overcome these difficulties we use asymptotic expansions of solutions in divergent series with respect to a small parameter . We obtain another system of equations which is more convenient to analyze. Our examples described in the previous section are based on this method and illustrate the behavior of ternary solid solutions (alloys), in particular: up-hill diffusion in , tracer diffusion in and tracer diffusion in a complex geometry in . We observe that solutions of the new system obtained by the expansion method with several terms, e.g. or , reflect an expected behavior of solutions of the original diffusion problem. There are many open questions in the theory and practice of small parameter expansions:

Whether or not it is possible to extend the small parameter method onto the case of singularities [11]? In the case of ternary diffusion problems singularities may be caused by irregularity of initial distributions u, v, w, e.g. Heaviside or smoothed peak type functions.

Is it possible to apply the small parameter method to a derivation of interfaces, e.g. via Matano’s method [32]?

If indeed replacing Neumann boundary conditions by Stefan free boundary conditions is easy to conduct?

Is it possible to generalize our method onto the case of media with variable diffusion coefficients and multi-phase media?

Is it possible to formulate and perform computations by means of a small parameter method within Onsager’s formalism? [33-35]