A new analytical formulation of retention effects on particle diffusion processes

A new analytical formulation of retention effects on particle diffusion processes

Luiz Bevilacqua^I; Augusto C.N.R. Galeão^II; Flavio P. Costa^III

^IUniversidade Federal do Rio de Janeiro, Instituto Alberto Coimbra, COPPE, Centro de Tecnologia, Bloco B, Sala 101, Ilha do Fundão, Caixa Postal 68506, 21945-970 Rio de Janeiro, RJ, Brasil

^IILaboratório Nacional de Computação Científica, Coordenação de Mecânica Computacional, Av. Getúlio Vargas, 333, Quitandinha, 25651-075 Petrópolis, RJ, Brasil

^IIIUniversidade Estadual de Santa Cruz, Departamento de Ciências Exatas e Tecnológicas, BR 415, Rodovia Ilhéus/Itabuna, km 16, 45662-900 Ilhéus, BA, Brasil

^{Correspondence to}

ABSTRACT

The ultimate purpose of this paper is to present a new analytical formulation to simulate diffusion with retention in a reactive medium under stable thermodynamic conditions. The analysis of diffusion with retention in a continuum medium is developed after the solution of an equivalent problem using a discrete approach. The new law may be interpreted as the reduction of all diffusion processes with retention to a unifying phenomenon that can adequately simulate the retention effect namely a circulatory motion. It is remarkable that the governing equation requires a fourth order differential term as suggested by the discrete approach. The relative fraction of diffusion particles β is introduced as a control parameter in the diffusion-retention law as suggested by the discrete approach. This control parameter is essential to avoid retention isolated from the diffusion process. Two matrices referring to material properties are introduced and related to the real phenomenon through the circulation hypothesis. The governing equation may be highly non-linear even if the material properties are constant, but the retention effect is a function of the concentration level, that is, β is a function of the concentration.

Key words: diffusion, retention, constitutive law, discrete approach, population dynamics, fourth order differential equations.

RESUMO

O objetivo último desse trabalho é apresentar uma nova formulação analítica para simular difusão com retenção em um meio reativo sob condições termodinamicamente estáveis. A análise da difusão com retenção em um meio contínuo é desenvolvida a partir da solução de um problema equivalente usando uma abordagem discreta. A nova lei pode ser interpretada como a redução de todos os processos de difusão com retenção a um fenômeno unificador que pode simular adequadamente a retenção. O propósito principal desse trabalho é apresentar uma nova formulação analítica para simular difusão com retenção em meio reativo termodinamicamente estável. A análise da difusão com retenção em um meio contínuo é desenvolvido a partir da solução de um problema equivalente usando uma abordagem discreta. A nova lei pode ser interpretada como a redução de todos os processos com retenção a um fenômeno unificado que pode simular adequadamente o efeito de retenção como um movimento circulatório. É notável que a equação que governa o fenômeno exige a inclusão de um termo diferencial de quarta ordem conforme sugerido pela abordagem discreta. A fração relativa β das partículas difundindo-se no meio é introduzida como um parâmetro de controle no processo de difusão-retenção. Esse parâmetro de controle é essencial para evitar retenção independente do processo de difusão. Duas matrizes contendo parâmetros relativos às propriedades materiais são introduzidas e relacionadas com a hipótese de circulação. A equação que governa o fenômeno pode ser altamente não linear mesmo se as propriedades materiais forem constantes desde que o efeito de retenção seja função da concentração, isto é que β seja função da concentração.

Palavras-chave: difusão, retenção, leis constitutivas, formulação discreta, equações diferenciais de quarta ordem.

INTRODUCTION

Spreading of particles or microorganisms immersed in a given medium or deployed on a given substrata is frequently modeled as a diffusion process. The mathematical formulation of the classical diffusion problem in a three-dimensional space x = [x₁ ,x₂ ,x₃ [^T can be written as follows:

defined in V for t > 0 together with the boundary and initial conditions:

where ℜ and ℑ represent respectively a differential operator and a function related to given initial conditions, p(x,t) stands for the particle concentration and K is the diffusion coefficient.

The system above represents quite satisfactorily the behavior of several physical phenomena related to dispersion processes. For some cases, however, the simplified approach above fails to represent the real physical behavior. Population spreading or dispersing particles may be partially and temporarily blocked when immersed in some particular media. An invading species, for instance, may hold a fraction of the total population stationary on the conquered territory in order to guarantee territorial domain. Chemical reactions may induce adsorption processes for the diffusion of solutes in liquid solvents in the presence of adsorbent material (Brandani et al. 2000). Multiphase flow, flows through porous media (Liu and Thompson 2002) and particulate systems are also among physicochemical phenomena that need improvement of the analytical formulation due to side effects not included in the classical diffusion modeling.

It is also important to remark that the determination of the diffusion parameters for more complex phenomena requires at least a reliable analytical representation. Certainly models are nothing but an approximation of the real world. However the more sophisticated the experimental setups, the more advanced the analytical formulation should be. There is evidence that particle retention may occur in diffusion processes for some dispersing substances immersed in particular supporting media (Atsumi 2002, D'Angelo et al. 2003, Deleersnijder et al. 2006, Green 1996, Liu and Thompson 2002, Muhammad 2004). The classical diffusion equation doesn't represent satisfactorily the retention phenomenon since there is a new phenomenon involved in the process that cannot be characterized simply by manipulating the diffusivity parameters. The main purpose of this paper is to present a new model that explicitly includes retention processes. The theory introduced here shows that a higher order partial differential equation is required that, consequently, involves the inclusion of new parameters. These new parameters, besides the diffusion coefficient, characterize the blocking process. Specific experimental techniques need to be settled to determine those complementary coefficients. Trying to solve the retention problem by imposing an artificial dependence of the diffusion coefficient on the particle concentration p(x,t) or introducing extra differential terms while keeping the second order rank of the governing equations, or both, disguises the real physical phenomenon occurring in the process. We will show that there is a more precise modeling that explicitly takes into account the retention effect in the dispersion process.

In the following sections the diffusion in a complex medium - admitting retention - is developed, first with a discrete approach and afterwards using a continuum approach. As will be explained later, in the context of the present theory retention is a dynamical rather than a stationary state.

POPULATION DYNAMICS: A DISCRETE APPROACH

This section is introduced here to make the paper self-contained. The formulation below is part of another paper, to appear shortly, dealing with the discrete approach of transport phenomena (Bevilacqua et al. in press).

Consider the process depicted in Figure 1. The rule governing the redistribution of contents of each cell indicates that a fraction of the contents αp_n is retained in the n^th cell and the exceeding volume is evenly transferred to the neighboring cells, that is, 0.5βp_n to the left, (n - 1)^th cell, and 0.5βp_n to the right, (n + 1)^th cell, at each time step, where β = (1 - α). This means that the dispersion runs slower than for the classical diffusion problem. Note that if β = 1, the problem is reduced to the classical Gaussian distribution.

Translating this rule into algebraic expressions we get:

It is now possible to manipulate Eq. (1a) and Eq. (1b) in order to obtain finite difference terms. It is important to notice that the rearrangement of terms must be done carefully otherwise we could obtain an equation that doesn't reproduce the required solutions for the limits β = 0 and β = 1. That is, the differential terms involved in the final expression may not reproduce the behavior required by the fundamental set of equations for the entire range of the parameters. Therefore, it is necessary to test the intermediate expressions at critical steps to verify if the underlying assumptions related to Eq. (1a) and Eq. (1b) are preserved. This means that, for β = 1, corresponding to no retention, the classical Gaussian distribution is obtained, that is, the contents of each cell is equally redistributed to the neighboring ones at each time step. If β = 0 on the other hand the solution is stationary that is, the total content of each cell remains the same all time.

Subtracting Eq. (1a) from Eq. (1b), after reducing all terms to the time t - 1, leads to:

After a relatively lengthy and careful algebraic manipulation - see appendix - the following expression is finally obtained:

Note that taking β = 1 Eq. (3) reduces to the classical diffusion problem, that is no retention, and with β = 0 Eq. (3) represents a stationary solution since the right hand side term vanishes. Consequently the time rate variation of the contents p_n equals zero for all t for all cells.

Calling K₂ = /2T₀ and K₄ = /4T₀, and taking the limit as Δx → 0 and Δt → 0, we obtain:

The fourth order term with negative sign introduces the retention effect. Clearly this term comes in naturally, without any artificial assumption, as an immediate consequence of the temporary retention imposed by the redistribution law. We would like to call the attention that similar equations can be obtained by introducing non-linear correlations in the Fick's law (Cohen and Murray 1981). These equations, however, are derived to take into account physical phenomena that hardly could be associated to the retention effect.

The result presented above justifies the following statement:

PROPOSITION. Suppose a one-dimensional diffusion process consisting of N particles moving in a homogeneous, isotropic supporting medium according to the Fick's law. If a fraction f(N) = (1 - β) of the diffusing particles is temporarily trapped along their trajectories due to some mechanical, physical, chemical or physicochemical interaction with the medium, the governing equation must include a fourth order differential term of the form:

where Γ is a time-independent material constant, (1 - β) is the fraction of the particles temporarily trapped, with 0 < β < 1 constant, and p(x,t) stands for the particle concentration in the medium.

It will be shown in the sequel that this proposition may be generalized in order to allow β and Γ to be time and space dependent parameters. In this case, however, a more complex differential equation is obtained including non-linear terms, and it would be necessary a more detailed discussion that we are going to postpone for a future paper.

PHENOMENOLOGICAL CONSIDERATIONS OF DIFFUSION PROCESSES

Before stating the laws that govern retention, it is worthwhile describing some basic ideas sustaining the formulation that will be introduced later on. Scientists and particularly engineers have frequently to deal with a class of problems that, behind an apparent simplicity, hides complex events. A problem with hidden complexity results from the simultaneous convergence of several phenomena involving a relatively large number of transformation processes. Most of the problems belonging to this class that appear in engineering practice require the determination of certain variables, that we will call here macro-state variables, representing complex physical, chemical or physicochemical interactions at micro scales. For the engineering point of view, within certain limits, it is not necessary to analyze the phenomenon at the smallest possible scale, because frequently the basic phenomena are not yet fully understood. It suffices knowing the behavior of the macro-state variables. That is, for a problem with hidden complexity, in several circumstances, the main task is to determine the relationship between the input variables and the output variables flowing in and out of a "black box".

The determination of a method or a theory adequate to establish the behavior of the macro-variables is by no means an easy task. The main difficulty arises from the fact that the theory must play a unification role. That is, the theory must as far as possible translate, with a relatively simple set of rules and laws, the response of the macro-state variables to the perturbations introduced in the system under consideration. More than that, a successful theory should give adequate response to several similar phenomena belonging to the same class. For instance, dispersion of gas, liquid or solid particles in different supporting media should be accurately modeled by the same theoretical framework because they belong to the same class of problems independently of the particle nature and of the medium. The macro-variable for this class of phenomenon is the particle concentration denoted by p(x,t). The theory proposed in this paper is derived from the following ideal process.

Consider a collection of particles immersed in some medium and moving in a preferred direction driven by repulsion forces such that they displace from high concentration regions towards low concentration regions free from constraints. The motion is not subjected to any critical barrier that could induce a significant perturbation in the process. Under this circumstance it is plausible to admit that the average speed of the particles is proportional to the concentration gradient, the particles being driven from regions with high concentration levels, strong repulsive forces, to regions with low concentration levels, weak repulsive forces. This is the fundamental idea supporting Fick's law. For the one-dimensional problem in a homogeneous medium, the absolute value of the diffusion speed may be defined as Dgrad(p(x,t)), where D is the diffusion coefficient.

The diffusion coefficient cannot be considered as a parameter characterizing all physicochemical phenomena associated to the diffusion processes. It is instead a coefficient associate to an ideal event that simulates satisfactorily a wide range of phenomena. The diffusion coefficient may have different physical interpretations depending on the particular phenomenon under consideration. Another way to see this methodology is to consider a "black box" that transforms input variables into output variables according to well-established laws, irrespectively from the physics that may be going on inside the box. In this sense, the Fick's law is a decomplexification process.

Now suppose that other events interfere in the diffusion mechanism. In this case the law must be modified accordingly. For instance, the effects of sources or sinks in a diffusion process cannot be modeled by some artificial modification of the flow velocity or distorting the interpretation of the diffusion coefficient. A new term is necessary to be incorporated in the governing equation. Diffusion with sinks or sources belongs to a different class of phenomena.

To avoid possible misunderstandings referring to more complex diffusion problems, we would like to remark that it is possible within a same class of phenomena to improve the respective laws incorporating extra terms. For the case of diffusion, for instance, Fick's law may be extended to incorporate higher order terms in the definition of the flux b(x,t). For the one-dimensional case Cohen and Murray, starting from a generalized law for the internal energy, arrived to the following expression for the flux:

This means that the driving forces inducing the particle motion are not related only to the gradient of the concentration grad(Ap(x,t)) as in Fick's law, but require the consideration of higher order terms. Since the corresponding governing equation includes for this case two new parameters B and k, it would be presumably possible to obtain a better representation of experimental results by adjusting properly thevalues of these two parameters. This is a refinement of the theoretical framework, but the basic event remains the same, belonging to the class of diffusion problems without taking into account any other phenomenon as retention for instance.

If however there is some essential modification of the phenomenology governing a given process, it is necessary to review the theoretical framework to find some universal model adequate to decomplexify the corresponding class of phenomenon. Unfortunately attempts to extend a theory fitted to represent a given class of events to a new phenomenological class are not seldom found in the literature. The results are certainly not universal and may fail to express a reliable response of the macro-state variables. This is the case of diffusion with temporary retention. The phenomenon of temporary retention cannot be adequately modeled by Fick's law or any refinement of the classical diffusion law. Diffusion with temporary retention belongs to another class of phenomena. A new suitable law could be stated axiomatically. It would be better to find a kind of universal event that may encompass a large number of phenomena and may also provide a consistent evaluation of the macro-state variables.

The main purpose of this paper is to propose a model that can be used to simulate with good accuracy the retention effect in a diffusion process. Based on that model, a new universal theoretical framework to decomplexify the temporary process is derived.

THE DIFFUSION CONSTITUTIVE LAW

The diffusion law adopted here is the classical Fick's law. Define Ψ_D as:

Ψ_D = d · n

where n is the unit vector normal to ∂V and d is the diffusivity vector. The diffusivity vector consistently with the Fick's law is defined as:

d = -Dgrad (p(x,t)).

Let us call D the diffusion matrix that is related to the physical properties of the medium, and p(x,t) is the mass concentration function. The components of D characterize the diffusion properties of the mass concentration in the embedding medium. For the general case, the matrix D is full (MacElroy 1996) and the diffusivity vector reads:

The minus sign comes from the well-known fact that under normal conditions the motion direction is oriented towards the regions of decreasing mass concentration.

THE RETENTION CONSTITUTIVE LAW

In this section we will introduce the key contribution of this paper, which is the derivation of the function Ψ_R(p(x,t)) representing the retention rate of the particle concentration trapped on the boundary layer in the vicinity of ∂V. The reasoning that will follow was guided by the result obtained with the discrete approach in the sense that the final equation governing diffusion with retention should include fourth order differential terms.

Instead of using an axiomatic approach through the definition of analytical expressions that wouldlead to the expected results, we will explore a plausible physical interpretation along with the development of the mathematical formulation. The physics of the problem postulated here is not unique, but is adequate to simulate a large range of cases. A consistent interpretation can also be given to the material parameters that will be introduced later on. The reader must have in mind that our purpose is to achieve a general law, that is, a law that would fit possibly all cases of diffusion including temporary retention. This means that the basic phenomenon introduced here to explain retention, has to be seen as a universal phenomenon or, more precisely, a universal virtual phenomenon.

The study of particular cases trying to simulate as close as possible the physical and chemical interactions governing the retention effect belongs to other approaches focusing very specific problems. Each problem requires a particular solution that cannot be extrapolated to other cases.

Two main ideas are enclosed in the definition of Ψ_R:

i) Recall that the terms representing the retention effect according to the previous derivation must be multiplied by β(1 - β) as shown in Eq. (4). To be consistent with the discrete approach, the effective mass temporarily trapped per unit of time Ψ_R must be proportional to the relative mass of the diffusing particles β as will become explicit in the next section from Eq. (5a). Clearly the factor β(1 - β) is a control parameter imposing variation limits on the blocking effect to avoid unacceptable mass accumulation on the boundary layer next to ∂V, where otherwise the density might grow beyond physically acceptable limits. Moreover it is also physically consistent to admit the retention effect directly proportional to the diffusing fraction. Indeed, the probability of particles being trapped increases with increasing intensity of particles diffusing in the medium.

ii) Now let us introduce a vector b for the retention state similar to the vector d for the diffusion state. The effective mass trapped in the boundary layer per unit of time, that is Ψ_R, is proportional to the component of the vector b orthogonal to the boundary surface ∂V. Considering the proportionality law proposed in i) above we may write:

Ψ_R = β(b · n).

The vector b will be called blocking vector. We introduce now the fundamental hypothesis to derive the formation law of the blocking vector.

First of all let us recall that we are dealing with temporary retention. Permanent retention corresponds to a sink that belongs to another class of phenomenon. Temporary retention may occur under different forms that are very difficult to represent analytically. As a matter of fact, if specific formulations for each kind of phenomenon are to be derived - temporary interstitial retention, temporary blocking by mechanical obstacles, chemical reactions, tortuosity of the substrata - it would be necessary to develop very complex theories, each one corresponding to a specific class of phenomenon.

What we intend to do here is to detect the key phenomenon that is common to all or at least to the majority of the diffusion processes under the effect of retention. That is, we are seeking for a unified model that represents the fundamental effects of the retention process for a large range of cases. The retention coefficients appearing in the evolution equation will in some cases represent the parameters corresponding to a virtual phenomenon simulating the real one.

The unifying phenomenon that will be assumed to simulate the retention effect for all cases is the circulation in a boundary layer in the neighborhood of the surface ∂V. That is, a fraction of the particles follows a helicoidal trajectory consisting of n spirals distributed on very thin boundary layer Δs in the neighborhood of the surface ∂V (Fig. 2a). While travelling along this helicoidal trajectory, the particle will be considered trapped in the boundary layer in the vicinity of the surface ∂V. The retention effect is therefore simulated by a temporary circulation in the neighborhood of the surface enclosing the volume V. The time delay for a particle to cross the boundary layer is proportional to the number of loops n and the angular velocity ω. In order to simplify the analysis, consider the loops compressed into a single closed trajectory contained in a plane Π tangent to the boundary. That is, the particle moves with velocity v on a circumferential path with radius a on the plane Π. To simulate the displacement of the particle along the direction parallel to n after the completion of n loops covering a distance Δs in a time Δt, the plane Π is assumed to move inwards the volume V with a velocity u. The retention effect is therefore simulated by a circulation with angular velocity ω = (|v|/a)n on the plane Π coupled with a translation inwards the volume V with velocity u. The modulus of u may be easily calculated |u| = Δs/Δt = ωλ/2π, and therefore u = (ωλ/2π)n. The modulus of u is very small since it is assumed that the pitch λ of the helicoidally motion is very small consistent with the hypothesis of retention within a very thin boundary layer. That is, the fraction of the particles moving on a temporarily closed trajectory (Fig. 2b) coupled with a very slow motion u of the circulation plane inwards the volume V will be considered trapped on the surface ∂V. Thus, the retention phenomenology defined here requires the particles to be temporarily confined in a boundary layer in the vicinity of the surface ∂V.

In order to describe the equivalent circulation motion, two different matrices of material parameters will be introduced. The first one denoted by C characterizes the velocity, of the equivalent circulatory motion. The second material matrix called R expresses the effectiveness of the blocking effect through the projection of the enclosed area by one loop on the tangent plane (t₁xt₂) of ∂V multiplied by the number of loops, such that the time delay equivalent to the temporary retention is achieved. Cases where a fraction of the dispersing particles is physically trapped for a certain amount of time on some material obstacle, as is the case of porous media, can be interpreted as a singularity where the radius of the circulation tends to zero and the motion is reduced to a spin. Description of some specific physical phenomena referring to the blocking time can be found in (Deleersnijder et al. 2006, Green 1996, Mota et al. 2004).

Now given the above hypothesis let us introduce the retention law.

1. The retention phenomenon occurring simultaneously with diffusion in a reactive medium under a state of stable thermodynamic equilibrium may be simulated by a circulatory motion defined by the rotational of an appropriate vector t given by:

rot(t) = Cgrad(p(x,t))

where C is a matrix of equivalent physical constants. That is the components c_ij of C refer to a circulation effect that adequately simulates the real phenomenon with respect to the dynamical and geometric characteristics of the circulatory motion.

2.The blocking matrix B is given by:

B = grad(rot(t) = grad(Cgrad(p(x,t)))

and the blocking vector b representing the retention rate is given by:

b = Rdiv(grad(Cgrad(p(x,t)))).

The components of the matrix R are also related to material properties that define the retention effectiveness that is the time delay imposed by the blocking effect. The matrix R stands for the control mechanism regulating the circulation intensity that is the number of loops in the equivalent spiralmotion. The experimental determination of C and R defines the equivalent circulation associatedwith the effective actual phenomenon.

3. The net retention rate Ψ_R corresponding to the particles temporarily trapped on the surface ∂V is given by the projection of b on the normal to the surface ∂V multiplied by the particles diffusing fraction β. That is:

Ψ_R = β(b · n).

Now taking the physical interpretation explained before for the retention effect as a circulation plus a translation, it is obvious that b = u. Note that we are dealing with a process evolving under a thermodynamic stable equilibrium. For more complex cases involving thermodynamic coupling the energy conservation principle must be taken into account. The theory for those cases introducing the retention effect has still to be developed.

The matrices C and R containing the physical parameters c_ij and r_ij depend on the interaction properties between supporting medium and diffusing substance. The determination of the dimensions of the variables intervening in the retention law may help to find the physical meaning of the components of C and R. Let the physical meaning of t be defined as the angular momentum of the trapped particles per unit of volume. The corresponding dimension is written as:

The notations M,L,T mean the dimensions of mass, length and time, respectively. Combining the expression above for the dimension of t with the definition of rot(t) introduced in the retention law we have:

From which follows the dimension of the components of the matrix C;

Note that this is the same dimension as for the diffusion coefficients composing the diffusion matrix D. Now the dimension of B = grad(rot(t)) comes immediately:

Now combining the dimension of b required by the principle of mass conservation with the definition of the blocking vector b introduced in the retention law, b = R div( B ), we obtain:

Therefore, the dimension of the coefficients of the matrix R is [R[ = [L[².

This is a plausible definition of the dimensions of the variables related to the retention law. From the dimensions of rot(t) and b, the dimensions of the components of C and R can be readily obtained, namely [L²/T[ and [L[², respectively, which is a an expansion speed and an area as suggested by the phenomenological assumption introduced before. The dimension of the components of the blocking matrix B shows that B may be interpreted as the time rate of the specific mass variation. It would be possible to generate other slightly different combinations of the fundamental dimensions but, if the concentration and concentration gradient are preserved as key elements any other alternative would not introduce anyessential deviation of the interpretation presented here.

The physical components c_ij and r_ij must be determined experimentally. These coefficients depend exclusively on the physicochemical or mechanical properties of the reactive medium interacting with the dispersing particles. One has to keep in mind that the experimental determination of C and R based on the present theory fits to a circulation effect simulating the actual retention phenomenon. The coefficients c_ij and r_ij may represent the real physical reaction depending on how the circulation hypothesis fits to the real case under consideration. This simplification, or, as stated before, this decomplexification process, is by no means a disadvantage of the present theory as compared with other approaches available in the specialized literature. The other theories are also based on a set of hypotheses that rarely explains the connection with the actual phenomenon. So the experimental coefficients obtained using other theories equally fails to represent those associated to the real case. The present formulation has the advantage of making it explicit the kind of simulation adopted namely the circulatory motion. Also the same equivalent phenomenon is taken for all cases that allows for the comparison of the coefficients c_ij and r_ij associated to different cases since all of them are modeled through the same universal representation.

It was mentioned before that the dimensions of the matrices D and C are the same. Also both appear in the formulations of the diffusion law and the retention law as well, as factors of the same quantity composing similar expressions, Dgrad(p(x,t)) and Cgrad(p(x,t)). This suggests that it would be plausible, in the absence of other more significant reason, to assume that they are equal, that is C = D. If this is true, the retention coefficients determining the temporary blocking process are reduced to those composing the matrix R. This would simplify considerably the experimental investigation. For sake of generality we will keep both independent in the sequel since there is no significant simplification from the theoretical point of view to make them equal.

Maybe most important of all is the introduction of higher order differential terms resulting from the present formulation of the retention law. This law is clearly inspired by the results obtained with the discrete approach without which it would be hard to establish a consistent representation of temporary retention. The presence of fourth order terms, consistent with the discrete approach, naturally follows from the law proposed above.

Now for the sake of completeness let us write the expanded expressions for the rot(t), B and b according to the above functional representation. In expanded notation we have successively:

Up to now we have left the theory as general as possible within the framework established for the representation of diffusion with retention. In this sense, the matrices D, C and R have been assumed to be complete. It is quite reasonable, however, for several cases, to assume that D, C and R are diagonals matrices, that is d_ij = 0, c_ij = 0 and r_ij = 0 for i ≠ j. This hypothesis holds for the case of orthotropic processes. For the more general case of skewed anisotropy the matrices D, C and R cannot be reduced to the diagonal terms.

For the sake of simplicity we will restrict ourselves in the sequel to orthotropic diffusion processes with retention. This will not impose any critical limitations to the theory developed here.

The diffusivity vector can now be written:

Similarly for the variables associated to the retention process we get:

For this case, the matrix B coincides with the Jacobian of rot(t), that is:

The vector b reads:

THE CONTINUUM APPROACH FORMULATION OF DIFFUSION WITH RETENTION

Now once established the basic mechanisms for the diffusion with retention, let us derive the governing equation assuming, as said before, that the process takes place under stable thermodynamic conditions. Consider a sufficient large set of elementary particles scattered in some reactive medium, that is, a medium with the property of delaying the motion of a certain number of particles for a given period of time as reported by some authors in physicochemical experiments (Atsumi 2002, D'Angelo et al. 2003, Deleersnijder et al. 2006, Green 1996, Liu and Thompson 2002, Muhammad 2004). Consider the open set V bounded by the surface ∂V at a given time t (Fig. 3). The total number of particles in V ∪ ∂V is divided into two groups. The first group comprises the particles belonging to V, and the set of particles distributed in the neighborhood of ∂V constitutes the second group. The specific mass of the dispersing particles, considered as a continuum, at that time t referred to an arbitrary point x = [x₁x₂x₃[^T is defined by p(x,t). For the sake of simplicity we will use the Cartesian system for the analytical development.

Call ΔM₁ the net mass variation inside the volume V in the time interval Δt. That is:

Now, suppose that the medium interacts with the elementary particles in such a way that a fraction of the particles on ∂V scatters through the medium, and the complementary fraction is trapped in a state that we will call temporarily retention state. Particles belonging to this state are captured back into the volume V crossing the boundary at very low speed and remaining temporarily within a thin boundary layer next to ∂V as explained in the previous section. In this sense we say that, while in the retention state, the particles are trapped. The time rate of mass variation per area unit given by Ψ_D characterizes the diffusion process, that is the outward particle flux through ∂V, and the time rate of mass variation per area unit given by Ψ_R characterizes the retention effect, that is the inward particle flux through ∂V. Both processes are active simultaneously on ∂V. Let us assume that diffusion and retention are independent processes homogeneously distributed on the surface ∂V covering two disconnected partitions of ∂V, that is, β∂V for diffusion and (1 - β)∂V for retention. This means that the particles on ∂V are either spreading - diffusion - or trapped - retention. Therefore, calling ΔM₂ the mass variation in the volume V due to the net flux through the boundary ∂V in the time interval Δt, we have:

The principle of mass conservation in the absence of sources or sinks reads (ΔM₁ + ΔM₂)/Δt = 0 or recaling Eq. (5a) and Eq. (5b):

As stated above, the variables Ψ_D and Ψ_R are the effective diffusion rate and the effective retention rate, respectively. Using the constitutive relations for diffusion and retention introduced in the previous sections, Eq. (5c) may be rewritten:

Now assuming the function p(x,t) sufficiently smooth and using the divergence theorem we get:

Introducing the expressions of b and d previously derived and taking the limit as Δt → 0 we obtain:

Since V is arbitrary and expanding the representation of B:

Or using the component notation:

Equation (7b) is the governing equation for the case of diffusion with retention. It differs significantly from some others approaches found in the technical literature to simulate this type of phenomenon. The theory developed above is consistent with the straightforward discrete derivation of the population dispersion with temporarily retention to colonize the occupied territory. Equation (7a) is in fact the generalized version of Eq. (4) corresponding to the discrete derivation. The determination of the physical components R, C and D depends on experimental efforts. Anyway, starting from the Eq. (7a) it should be possible to devise meaningful experimental strategies since it is clear what should be detected and measured.

For one-dimensional diffusion with retention along the x axis, the Eq. (7a) reduces to a simpler expression, namely:

Despite that the above equation is considerably simpler, it keeps some very fundamental information that applies in general for more complex mass diffusion problems. Particularly interesting is the role of the parameter β. Note first that the diffusion term grows steadily with β, while the retention term varies according to the logistic equation (Fig. 4). As a consequence, low intensity diffusion processes do not mean high intensity retention. The control of the retention process introduced in the governing equation avoids the existence of pure retention as expected, but allows for diffusion only, as it is also consistent with the physics of the process.

Now it is also reasonable to expect that the retention mechanism saturates after some time. This idea may be introduced by defining the parameter β as a function of p(x,t). Let us assume β, a power function of the saturation quotient q(x,t) = p(x,t)/p^* , that is β = [q(x,t)[^m, where p^* is the saturation limit of the mass concentration. With this expression for β and with r, c and k constants, the Eq. (7c) gives:

where the coefficients ϕ_i, i = 1,2,3,4 are given by:

Figure 5 illustrates the variation of β with q(x,t) for several values of m. For m < 1 the retention effect reaches the maximum, that is β = 0.5 for very low saturation quotient. That is, for relatively low levels of mass concentration, the capacity of the supporting medium to block diffusion is close to the upper limit. For m > 1 the situation is reversed. The maximum blocking effect is reached for relativelyhigh saturation quotient (Fig. 5).

For the case of linear correlation between β and q, that is m = 1, the Eq. (7d) reduces to:

Close to the maximum retention effect that is when q(x,t) is in the neighborhood of 0.5, we may write q(x,t) ≈ 0.5 + ε(x,t). Then, close to this point, the Eq. (7e) can be written:

CONCLUDING REMARKS

The phenomenon of diffusion with retention has been already observed and reported by several authors. Analytical formulations have also been proposed by some authors based on certain assumptions focused on the particular problem under consideration. All of them assume consistently the Fick's law as the analytical representation of the diffusion process even in the presence of retention effects. Consequently, the governing equations for diffusion with retention appearing in the literature so far, to the best of our knowledge, belong to the class of the second order partial differential equation as in the classical diffusion approach. They are as a matter of fact modified versions of the classical diffusion equation. As a consequence it would be virtually impossible to derive a general theory since each problem requires a particularapproach. The classical solution could apply with the introduction of sources and sinks, such that theoverall effect on the number of particles present in the medium remains unaltered after sufficiently long time. But as explained before, this is a different phenomenon since it requires energy exchange with the surrounding medium.

We believe that the diffusion through the brain intracellular space (Nicholson and Sykova 1998, Mota et al. 2004) may be analyzed with the theory proposed here. The intracellular spaces, that is, the substratum where neuroactive substances and molecules circulate, introduce a three-dimensional diffusion processdue to the tortuosity of the free trajectories. If the tortuosity acts as an obstacle introducing only a kind of circulation orthogonal to the main diffusion direction, the process may be simulated by a temporary delay and the approach introduced here may represent this type of phenomenon. This is a typical example of a complex problem being simulated by a general decomplexfication model.

This paper shows that the second order differential equation is not adequate to represent the physics of diffusion with retention. The retention process requires a fourth order term multiplied by a control factor β(1 - β) that is essential to avoid the formation of unattainable mass concentration. It is remarkable that the solution of the discrete problem of population invasion with temporary territory occupation provided the fundamental clue to derive the solution of the respective continuous problem. It has been shown very clearly that a fourth order term appears in a straightforward manner to represent the retention phenomenon. We believe that this argument suffices to prove the need to introduce this term to achieve an accurate model for diffusion with retention. Problems with a source term are even more sensitive to the presence of the fourth order term. Indeed the stability conditions cannot be correctly formulated without the consideration of the retention coefficients r_ij and c_ij. Moreover, the redistribution parameter b plays a very important role on the stability conditions. Indeed, for β ≈ (1 - ε) where ε is small, solutions of the type ζ = ζ₀ exp(2πix/λ - ωt) may cause instability for a very large range of frequencies independently of the values of the diffusion and retention coefficients (Bevilacqua et al. in press).This means that a state of relative intense diffusion coupled with an external source is very sensitive to small retention perturbations for a very large frequency range.

Fourth order terms may also appear when the Fick's law is expanded to introduce nonlinear terms (Cohen and Murray 1981). This approach, however, has no connection with the retention phenomenon. Moreover, the nonlinear expansion of the Fick's law introduces always nonlinear terms in the governing equation and ignores the control parameter multiplying the fourth order term. The approach presented here doesn't require necessarily the introduction of nonlinear terms in the governing equation.

It is worthwhile mentioning the occurrence of fourth order terms in other engineering and physical problems as those related to the determination of thin films profiles. Examples are molecular beam epitaxy (Barabási and Stanley 1995) and surface coating with viscous fluids (Schwartz and Roy 2004, Myers 1998, Mullins 1957). Fourth order term also appears in the partial differential equation for the determination of the entropy evolution connected to the Shannon information density (Broadbridge 2008).

The theory presented here is satisfactory for moderate temperatures and small temperature variations. The extended theory including thermodynamic coupling has still to be developed. Classical approaches can be found in the current literature (Kirkland et al. 1990, Eon et al. 1997, Fitch et al. 1996, Schimpf and Giddings 1989), but there is no satisfactory theory including retention effects. The determination of the functional dependence of β on mass concentration, on temperature and on other intervening molecular forces is, for instance, a challenge for further investigation. We believe, however, that, before extending the theory to include energy equations, the experimental validation of the blocking effect as introduced here, for thermodynamically stable processes, should be done. To concentrate efforts to find experimental evidence of the functional dependence of β on the saturation quotient is of particular importance.The terms of the matrices C and R have also to be determined experimentally. As a first step it is advisable to focus the attention on homogeneous media and one-dimensional flow. For these cases, the terms of matrices R and C come together (rc) which makes the experimental analysis easier despite being limited regarding the phenomenology inherent in the process.

Analytical solution of the equations derived above could be found maybe only for very particularcases. We believe that efforts should be concentrated on qualitative behavior, asymptotic and numerical solutions. Similar equations have been studied by several authors (You and Kavet 2000, Lu et al. 2007, Rubinstein et al. 1989, Bernoff and Bertozzi 1995, Myers 1998).

We believe that the theory presented here is consistent, compatible with a very robust discrete approach and represents satisfactorily the retention phenomenon in reactive media. We have introduced a general retention matrix for the sake of completeness. Possibly the retention matrix is also a diagonal matrix,which simplifies considerable the formulation.

Finally other type of phenomena in reactive media as counter-diffusion may occur (Brandani et al. 2000). For this case a reactive medium is one that has the property of triggering a mechanism thatpartially reverses the flow direction. This means that particle dispersion may occur simultaneously as diffusion processes in opposite directions. This mechanism may be thought of as being induced by the presence of ionic channels distributed in the support medium that are able to partially reverse the flow direction. The governing equations for this case may be derived in a similar way. Fourth order terms, however, may not appear in the governing equations.

ACKNOWLEDGMENTS

This research was partially supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) Research Fellowship granted to the first two authors, by the CNPq Research Grant"Universal 14/2009" and by Fundação de Amparo à Pesquisa do Estado da Bahia (FAPESB).

Manuscript received on March 25, 2010; accepted for publication on September 23, 2010

APPENDIX

SYMMETRIC DISTRIBUTION WITH RETENTION

Consider the symmetric distribution with retention as displayed in Figure 1. The distribution scheme can be written as:

Reduction of the right hand side of Eq. (A1b) to time t - 1 leads successively to:

Subtracting given by (A1a) we get successively:

Note that the left hand side of Eq. (A2) represents the difference in time (t + 1 and t) of contents for the same cell n. In general the time associated to the term t + k means t + kΔt, where Δt is the time increment. All terms on the right hand side are referred to the same time t - 1, but involve sum and difference operations of cell's contents for cells ranging from n - 2 to n + 2.

The term on the left hand side is the first-order difference with respect to the variable t. We assign this difference with the notation:

There is one term on the right hand side of Eq. (4) of the form , and another one ofthe form . These expressions are the second- and fourth-order differences, respectively, with respect to the space variable, centered at the cell n. We assign these differences with the notation:

The Eq. (A2) can now be written in terms of the first and second order differentials:

Now assume = p(x_n,t) where p(x,t) is sufficiently differentiable with respect to the space variables and a smooth function of t. That is, we want to consider Eq. (4) as the discrete formulation of certain PDE.

For p(x) sufficiently differentiable, the first-order difference referred to x_n and x_n+1, respectively, differ by a quantity of order Δx² where Δx = |x_n+1 - x_n|. Certainly for p(x,t) the previous calculation is true provided that time is kept constant, that is:

In a similar way it can be shown that the deviation of the second order differences centered at points x_n-1 and x_n, respectively, is of the order of Δx³, that is O(Δx²). That is:

Introducing these expressions in (A3) we obtain after some simple operations:

or:

where Δ²p_n and Δ⁴p_n are the second and fourth order difference expressions. Let us take:

where L₀, L₁ and T₀ are scale factors, and Δx = L₀/m = L₁/ and Δt = T₀/m² are incrementsize and time interval, respectively. Substituting the above relations in (A5) we get:

The scale factors L₀, L₁ and T₀, together with the parameter m, provide very useful clues to define the sizes of space increment and time step for numerical integration of the finite difference equation.

• ATSUMI H. 2002. Hydrogen bulk retention in graphite and kinetics of diffusion. J Nuc Mat 307/311: 1466-1470.
• BARABÁSI AL AND STANLEY HE. 1995. Fractal Concepts in Surface Growth. Cambridge University Press.
• BERNOFF AJ AND BERTOZZI A. 1995. Singularities in a modified Kuramoto-Sivashinsky equation describing interface motion for phase transition. Physics D 85: 375-404.
• BEVILACQUA L, GALEÃO ACNR AND PIETROBON F (in press). On the Significance of Higher Order Differential Terms in Transport Processes. J Braz Soc Mech Sci Eng.
• BRANDANI S, JAMA M AND RUTHVEN D. 2000. Diffusion, self-diffusion and counter-diffusion of benzene and p-xylene in silicalite. Micropor Mesopor Mat 35/36: 283-300.
• BROADBRIDGE P. 2008. Entropy Diagnostics for Fourth Order Partial Differential Equations in Conservation Form. Entropy 10: 365-379.
• COHEN DS AND MURRAY JM. 1981. A Generalized Model for Growth and Dispersal in a Population. J Math Biol 12: 237-249.
• D'ANGELO MV, FONTANA E , CHERTCOFF R AND ROSEN M. 2003. Retention phenomena in non-Newtonian fluids flow. Physics A 327: 44-48.
• DELEERSNIJDER E, BECKERS J-M AND DELHEZ EM. 2006. The Residence Time of Setting in the SurfaceMixed Layer. Environ Fluid Mech 6: 25-42.
• EON SJ, SCHIMPF ME AND NYBORG A. 1997. Compositional effects in the retention of colloids by thermal field-flow fractionation. Anal Chem 69(17): 3442-3450.
• FITCH A, SONG J AND STEIN J. 1996. Molecular Structure Effects on Diffusion of Cations in Clays. Clays Clay Miner 44(3): 370-380.
• GREEN PF. 1996. Translational Dynamics of Macromolecules in Metals 251 p. In: NEOGI P (Ed), Diffusion in Polymers. Marcel Dekker Inc.
• KIRKLAND JJ, BOONE LS AND YAU WW. 1990. Retention effects in thermal field-flow fractionation. J Chromatogr A 517: 377-393.
• LIU H AND THOMPSON KE. 2002. Numerical modeling of reactive polymer flow in porous media. Comput Chem Eng 26: 1595-1610.
• LU L, WANG M AND LAI C-H. 2007. Image denoise by fourth order PDE based on the changes of Lapacian.J Algorithm Comput Tech 2(1): 99-110.
• MACELROY JMD.1996. Diffusion in Homogeneous Media 1 p. In: NEOGI P (Ed), Diffusion in Polymers. Marcel Dekker Inc.
• MOTA M, TEIXEIRA JA, KEATING JB AND YELSHIN A. 2004. Changes in diffusion through the brain extracellular space. Biotechnol Appl Biochem 39: 223-232.
• MUHAMMAD N. 2004. Hydraulic, Diffusion, and Retention Characteristics of Inorganic Chemicals in Bentonite, PhD dissertation, Department of Civil and Environmental Engineering, College of Engineering, University of South Florida, 251 p.
• MULLINS WW. 1957. Theory of Thermal Grooving. J Appl Phys 28(3): 333-339.
• MYERS TG. 1998. Thin films with high surface tension. SIAM Review 40(3): 441-462.
• NICHOLSON C AND SYKOVA E. 1998. Extracellular space structure revealed by diffusion analysis. TINS 21: 207-215.
• RUBINSTEIN J, STERNBERG P AND KELLER J. 1989. Fast reaction, Slow Diffusion, and Curve Shortening, SIAM. J Appl Math 49(1): 116-133.
• SCHIMPF ME AND GIDDINGS JC. 1989. Characterization of thermal diffusion in polymer solutions by thermal field-flow fractionation: Dependence on polymer and solvent parameters. J Polym Sci Part B: Polym Phys 27(6): 1317-1332.
• SCHWARTZ LW AND ROY V. 2004. Theoretical and numerical results for spin coating of viscous liquids. Phys Fluids 16(3): 569-584.
• YOU Y-L AND KAVET M. 2000. Fourth-Order Partial Differential Equations for Noise Removal. IEEE TransImage Process 9(10): 1723-1730.

Publication Dates

• Publication in this collection
  30 Sept 2011
• Date of issue
  Dec 2011

History

• Received
  25 Mar 2010
• Accepted
  23 Sept 2010