Transport processes in porous media. 2. Numerical modeling

Călin Vamoş
(original), paper

Abstract

The  paper proposes a numerical model for transport in heterogeneous porous media, built on the background of the continuous modeling from the first part of this work. The macroscopic behaviour of a microscopic particles ensamble is obtained by numerical simulation of their microscopic motion, in the molecular dynamics manner [Kiplik and Banavar, 1995].
The particles motion is governed by a random walk on a grid, similarly to the cellular automaton presented by Nishidate and Baba [1996]. The macroscopic description is given by space-time coarse-grained averages which provide a continuous description of the system [Vamos et al., 1996, our first paper, in this issue]. A first test was achieved by an accurate numerical soluiton of the one-dimensional diffusion equation. The number of particles and the averaging space-time scale needed  for a macroscopical description of the diffusion process with a given precision and the behaviour of systems with small concentrations are discussed in [Vamos et al., 1997b]. The model for diffusion in random environments was obtained by embededing the particles system into a random advection field. Numerical results are in good agreement with analytical ones obtained by Matheron and de Marsily [1980], using their model for statified aquifers.

Authors

C. Vamos
Tiberiu Popoviciu Institute of Numerical Analysis

N. Suciu
Tiberiu Popoviciu Institute of Numerical Analysis

U. Jaeckel
Forschungszentrum Julich GmbH, Institut fur Chemie und Dynamik der Geosphare, Deutschland

H. Vereecken
Forschungszentrum Julich GmbH, Institut fur Chemie und Dynamik der Geosphare, Deutschland

Keywords

Paper coordinates

C. Vamoş, N. Suciu, U. Jaeckel, H. Vereecken (1998), Transport processes in porous media. 2. Numerical modeling, Rom J. of Hydr. & Water Resour., 5(1-2), 85-97

References

PDF

Scanned paper.

About this paper

Journal

Rom. J. Hydr. & Water Resour.

Publisher Name
DOI

Not available yet.

Print ISSN

Not available yet.

Online ISSN

Not available yet.

Google Scholar Profile

soon

[1] C. W Gardiner,  1983, Handbook of Stochastic Methods (for  Physics), Chemistry and Natural Science),  Springer-Verlag, New York.
[2] J. Koplik and J. R. Banavarm, 1995,  Continuum deductions form molecular hydrodynamics,  Annu. Rev. Fluid Mech. 27, 257-292.
[3] L D. Landau and E. M. Lifchitz, 1988,  Statistical Physics (in romanian), Ed. Tehnică, București.
[4] G. Matheron and G. de Marsily, 1980,  Is Transport in porous media always diffusive?  Water Resour. Res. 16, 901-917.
[5] K. M. Nishidate and R. J. Gaylord, 1996, Cellular automaton model for random walkers,  Phys. Rev. Lett., 77 (9), 1675-78.
[6] L. Z. Rumshinski, 1974,  Mathematical Processing of Experimental Data (in romanian), Ed. Tehnică, București.
[7] H.  Schwarze, U. Jaekel and H. Vereecken, 1998,  Estimation of macrodispersion by different approximation methods for flow and transport in randomly hetergeneous media (to be published).
[8] N. Suciu, C. Vamoș, A. Georgescu, U. Jaekel and H. Vereecken, 1998, Transport processes in porous media. 1. Continuous modeling (this issue).
[9] B. F. A Tompson and L. W. Gelhar, 1990, Numerical simulation in three-dimensional randomly hetergeneous media,  Water Resour. Res. 26 (10), 2541-2562.
[10] C. Vamoș, A. Georgescu, N. Suciu and I. Turcu, 1996a, Balance equations for physical with corpuscular structure, Physica A, 227, 81-92.
[11] C. Vamoș, A. Georgescu and N. Suciu, 1996b, Balance equations for a finite number of particles,  St. Cerc. Mat., 48, 115-127.
[12] C. Vamoș, N. Suciu and A. Georgescu, 1997a, Hydrodynamic equations for one-dimensional systems of inelastic particles, Phys. Rev. E. 55, 6277-6280.
[13] C. Vamoș, N. Suciu and M. Peculea, 1997b,  Numerical modelling of the one-dimensional diffusion by random walkers, Rev. Anal. Numer. Theorie Approximation, 26 (1-2), 237-247.
[14] N. G.. van Kampen, 1981,  Stochastic Processes in Physics and Chemistry,  North-Holland, Amsterdam.

1998b Suciu-Vamos-Jaekel-Vereecken - Transport processes in porous media-Numerical modeling

TRANSPORT PROCESSES IN POROUS MEDIA. 2. NUMERICAL MODELING

C. VAMOŞ a a ^(a){ }^{a}a, N. SUCIU a a ^(a){ }^{a}a, U. JAEKEL b b ^(b){ }^{b}b, H. VEREECKEN b b ^(b){ }^{b}ba) "Tiberiu Popoviciu" Institute of Numerical Analysis, Romanian Academy, P. O. Box 68, 3400 Cluj-Napoca 1, România, e-mail: nsuciu@ictp.math.ubbcluj.ro, cvamos@ictp.math.ubbcluj.rob) Forschungszentrum Jülich GmbH, Institut für Chemic und Dynamik der Geosphäre, ICG-4, D52425-Jülich, Deutschland,c-mail: H.Vereecken@fz-juelich.de, U.Jaekel@fz-juelich.de

Abstract

The paper proposes a numerical model for transport in heterogeneous porous media, built on the background of the continuous modeling from the first part of this work. The macroscopic behavior of a microscopic particles ensamble is obtained by numerical simulation of their microscopic motion, in the molecular dynamics manner [Koplik and Banavar, 1995]. The particles motion is governed by a random walk on a grid, similarly to the cellular automaton presented by Nishidate and Baba [1996]. The macroscopic description is given by space-time coarse-grained averages which provide a continuous description of the system [Vamoss et al., 1996, our first paper, in this issue]. A first test was achieved by an accurate numerical solution of the one-dimensional diffusion equation. The number of particles and the averaging space-time scale needed for a macroscopical description of the diffusion process with a given precision and the behavior of systems with small concentrations are discussed in [Vamoş et al., 1997b]. The model for diffusion in random environments was obtained by embedding the particles system into a random advection field. Numerical results are in good agreement with analytical ones obtained by Matheron and de Marsily [1980], using their model for statified aquifers.

1 THE MICROSCOPIC DEFINITION OF THE CONCENTRATION FIELD

To discuss some problems related to the microscopic definition of the concentration, we present a method often used in thermodynamics [Landau and Lifchitz, 1988]. Consider N N NNN molecules in a volume V V VVV. The concentration at the point of the
position vector r r r\mathbf{r}r at the moment t t ttt is defined as the number of particles per unit volume in a domain of volume V < V V < V V < V\mathcal{V}<VV<V centered at this point. The domain should be spherical, otherwise the concentration would depend not only on the position but also on the orientation of the domain. We denote by S ( r , a ) S ( r , a ) S(r,a)S(\mathbf{r}, a)S(r,a) the sphere of center r r r\mathbf{r}r and radius a a aaa. Then the concentration is defined as the function c : R 3 × R → R + c : R 3 × R → R + c:R^(3)xxRrarrR_(+)c: \mathbf{R}^{3} \times \mathbf{R} \rightarrow \mathbf{R}_{+}c:R3×R→R+ given by
(1.1) c ( r , t ) = 1 V ∑ i = 1 N H + ( a 2 − ( r i ( t ) − r ) 2 ) (1.1) c ( r , t ) = 1 V ∑ i = 1 N   H + a 2 − r i ( t ) − r 2 {:(1.1)c(r","t)=(1)/(V)sum_(i=1)^(N)H^(+)(a^(2)-(r_(i)(t)-r)^(2)):}\begin{equation*} c(\mathbf{r}, t)=\frac{1}{\mathcal{V}} \sum_{i=1}^{N} H^{+}\left(a^{2}-\left(\mathbf{r}_{i}(t)-\mathbf{r}\right)^{2}\right) \tag{1.1} \end{equation*}(1.1)c(r,t)=1V∑i=1NH+(a2−(ri(t)−r)2)
where r i ( t ) , t ∈ R r i ( t ) , t ∈ R r_(i)(t),t inR\mathbf{r}_{i}(t), t \in \mathbf{R}ri(t),t∈R, is the position vector of the molecule i N i N i quad Ni \quad NiN at the moment t t ttt. The left continuous Heaviside function H + ( a 2 − ( r i ( t ) − r ) 2 ) H + a 2 − r i ( t ) − r 2 H^(+)(a^(2)-(r_(i)(t)-r)^(2))H^{+}\left(a^{2}-\left(\mathbf{r}_{i}(t)-\mathbf{r}\right)^{2}\right)H+(a2−(ri(t)−r)2) is equal to 1 if the molecule i i iii is inside the sphere S ( r , a ) S ( r , a ) S(r,a)S(\mathbf{r}, a)S(r,a) and vanishes otherwise. The function defined by (1.1) for given a a aaa and r i r i r_(i)\mathbf{r}_{i}ri, is a finite linear combination of Heaviside functions having null derivatives except, when their argument vanishes and the derivative does not exist. So (1.1) defines a step function which can not satisfy a partial derivative equation of the diffusion equation type. In thermodynamics one considers that for a large enough N N NNN the function c ( r , t ) c ( r , t ) c(r,t)c(\mathbf{r}, t)c(r,t) given by (1.1) is well approximated by a continuous field. To obtain the condition that N N NNN should satisfy we analyze the simple case of the thermodynamical equilibrium state in absence of the exterior fields. Then the molecules are uniformly distributed in the V V VVV volume and, taken as a continuous field, the concentration has a constant value at any point in the volume and at any moment, c o ( r , t ) = N / V c o ( r , t ) = N / V c_(o)(r,t)=N//Vc_{o}(\mathbf{r}, t)=N / Vco(r,t)=N/V. If we want to verify this equilibrium distribution by definition (1.1), we count the number of molecules n n nnn at the moment t t ttt in the sphere S ( r , a ) S ( r , a ) S(r,a)S(\mathbf{r}, a)S(r,a). It is obvious that the result is affected by fluctuations and the measured concentration e = n / V e = n / V e=n//Ve=n / \mathcal{V}e=n/V differs from c o c o c_(o)c_{o}co. In [Landau and Lifchitz, 1988, section 114], one shows that n n nnn satisfies a Poisson repartition with dispersion σ = n ¯ σ = n ¯ sigma=sqrt bar(n)\sigma=\sqrt{\bar{n}}σ=n¯, where n ¯ = N V / V n ¯ = N V / V bar(n)=NV//V\bar{n}=N \mathcal{V} / Vn¯=NV/V is the mean number of molecules in V V V\mathcal{V}V. So, the dispersion of the concentration fluctuation is equal to σ c = N / V V σ c = N / V V sigma_(c)=sqrt(N//VV)\sigma_{c}=\sqrt{N / \mathcal{V} V}σc=N/VV. If N N NNN is large enough, Δ c = c − c o Δ c = c − c o Delta c=c-c_(o)\Delta c=c-c_{o}Δc=c−co has a normal repartition of zero average and dispersion σ c σ c sigma_(c)\sigma_{c}σc. According to the "three sigma" rule of excluding rough errors [Rumshinski, 1974], we impose that the relative error should be smaller than a value ε ε epsi\varepsilonε fixed with a confidence level of 0.997 and we obtain
(1.2) 3 σ c ◻ ε c o ⟹ V V ≥ 9 ε 2 N (1.2) 3 σ c ◻ ε c o ⟹ V V ≥ 9 ε 2 N {:(1.2)3sigma_(c)◻epsic_(o)Longrightarrow(V)/(V) >= (9)/(epsi^(2)N):}\begin{equation*} 3 \sigma_{c} \square \varepsilon c_{o} \Longrightarrow \frac{\mathcal{V}}{V} \geq \frac{9}{\varepsilon^{2} N} \tag{1.2} \end{equation*}(1.2)3σc◻εco⟹VV≥9ε2N
For a large enough N N NNN, this formula gives the minimum volume V V V\mathcal{V}V (or the minimum radius a a aaa ) necessary to measure the concentration with the precision ε ε epsi\varepsilonε, i.e. the space scale for which the measured concentration behaves like a continuous field with approximation ε ε epsi\varepsilonε. If (1.2) is satisfied, we can write c ( r , t ) ∼ c o ( r , t ) + O ( ε ) c ( r , t ) ∼ c o ( r , t ) + O ( ε ) c(r,t)∼c_(o)(r,t)+O(epsi)c(\mathbf{r}, t) \sim c_{o}(\mathbf{r}, t)+\mathcal{O}(\varepsilon)c(r,t)∼co(r,t)+O(ε), for ε ⟶ 0 ε ⟶ 0 epsi longrightarrow0\varepsilon \longrightarrow 0ε⟶0.
The classical definition of the concentration is not applicable when N N NNN or V V V\mathcal{V}V is too small. To exemplify, we consider the case when there is a single molecule in the volume V V VVV. Then c ( r , t ) = V − 1 c ( r , t ) = V − 1 c(r,t)=V^(-1)c(\mathbf{r}, t)=\mathcal{V}^{-1}c(r,t)=V−1 for r ∈ S ( r 1 ( t ) , a ) r ∈ S r 1 ( t ) , a rin S(r_(1)(t),a)\mathbf{r} \in S\left(\mathbf{r}_{1}(t), a\right)r∈S(r1(t),a), and in the rest c c ccc vanishes.
Therefore the concentration is completely different from the equilibrium concentration which is as well c o ( r , t ) = 1 / V c o ( r , t ) = 1 / V c_(o)(r,t)=1//Vc_{o}(\mathbf{r}, t)=1 / Vco(r,t)=1/V in the entire volume V V VVV. These difficulties occur because in definition (1.1) one implicitly supposes an instantaneous measurement of the molecules number in the volume V V V\mathcal{V}V. The actual measurement has a duration which defines the temporal scale as the volume V V V\mathcal{V}V (or the radius a a aaa ) defines space scale. If we denote by ( t − τ , t + τ ) ( t − τ , t + τ ) (t-tau,t+tau)(t-\tau, t+\tau)(t−τ,t+τ) the averaging interval, we define the concentration by
(1.3) c ( r , t ) = 1 2 τ V ∑ i = 1 N ∫ t − τ t + τ H + ( a 2 − ( r i ( t ′ ) − r ) 2 ) d t ′ (1.3) c ( r , t ) = 1 2 τ V ∑ i = 1 N   ∫ t − τ t + τ   H + a 2 − r i t ′ − r 2 d t ′ {:(1.3)c(r","t)=(1)/(2tauV)sum_(i=1)^(N)int_(t-tau)^(t+tau)H^(+)(a^(2)-(r_(i)(t^('))-r)^(2))dt^('):}\begin{equation*} c(\mathbf{r}, t)=\frac{1}{2 \tau \mathcal{V}} \sum_{i=1}^{N} \int_{t-\tau}^{t+\tau} H^{+}\left(a^{2}-\left(\mathbf{r}_{i}\left(t^{\prime}\right)-\mathbf{r}\right)^{2}\right) d t^{\prime} \tag{1.3} \end{equation*}(1.3)c(r,t)=12τV∑i=1N∫t−τt+τH+(a2−(ri(t′)−r)2)dt′
In [Vamos et al. 1996, a , b a , b a,b\mathrm{a}, \mathrm{b}a,b ] one proves that, if the functions r i ( t ) r i ( t ) r_(i)(t)\mathbf{r}_{i}(t)ri(t) are at least piecewise analytic, then the function defined by (1.3) has a. e. continuous first order partial derivatives. That is, the temporal averaging transforms the step function (1.1) into a continuous field even if discontinuities remain in the first order partial derivatives.
The new definition (1.3) is meaningful even when V V V\mathcal{V}V or N N NNN is very small and the classical definition (1.1) can not be applied. Consider again the thermodynamical equilibrium state with an uniform distribution of the molecules in the volume V V VVV. The number of molecules in V V V\mathcal{V}V is measured at each moment over the interval ( t − τ , t + τ t − τ , t + τ t-tau,t+taut-\tau, t+\taut−τ,t+τ ) and then it is averaged. We attach to this continuous temporal average a discrete one, to which we can apply the same Poisson distribution. We denote by Δ t Δ t Delta t\Delta tΔt the mean time interval over which the molecule remains inside the volume V V V\mathcal{V}V. Consider the averaging interval ( t − τ , t + τ t − τ , t + τ t-tau,t+taut-\tau, t+\taut−τ,t+τ ) divided into 2 τ / Δ t 2 τ / Δ t 2tau//Delta t2 \tau / \Delta t2τ/Δt subintervals of Δ t Δ t Delta t\Delta tΔt length and suppose that the existence of the molecule in volume V V V\mathcal{V}V within an subinterval Δ t Δ t Delta t\Delta tΔt is independent from its existence in the same volume within another subinterval Δ t Δ t Delta t\Delta tΔt. Then the concentration fluctuations (1.3) for N N NNN molecules over a time interval of 2 τ 2 τ 2tau2 \tau2τ length, are equivalent to the fluctuations of 2 N τ / Δ t 2 N τ / Δ t 2N tau//Delta t2 N \tau / \Delta t2Nτ/Δt molecules in a Δ t Δ t Delta t\Delta tΔt interval. Therefore instead of formula (1.2) we have
(1.4) V τ ≥ 9 △ t 2 ε 2 c o (1.4) V τ ≥ 9 △ t 2 ε 2 c o {:(1.4)V_(tau) >= (9/_\t)/(2epsi^(2)c_(o)):}\begin{equation*} \mathcal{V}_{\tau} \geq \frac{9 \triangle t}{2 \varepsilon^{2} c_{o}} \tag{1.4} \end{equation*}(1.4)Vτ≥9△t2ε2co
This formula expresses the relation between the space scale ( V ) ( V ) (V)(\mathcal{V})(V) and the temporal one ( τ ) ( τ ) (tau)(\tau)(τ), necessary to obtain a continuous description of the concentration with a precision ε ε epsi\varepsilonε. The increase of the temporal scale can compensate the decrease of the space one. Unlike (1.2), the relation (1.4) is valid for any N N NNN and V V V\mathcal{V}V.

2 THE CELLULAR AUTOMATON NUMERICAL ALGORITHM

We consider a one-dimensional space lattice { x k = k δ x ∣ − m k m } x k = k δ x ∣ − m k m {x_(k)=k delta x∣-m quad k quad m}\left\{x_{k}=k \delta x \mid-m \quad k \quad m\right\}{xk=kδx∣−mkm}, where δ x δ x delta x\delta xδx is the space step of the lattice. On the lattice there are N N NNN noninteracting particles
which move according to the random walk law. If the i i iii-th particle is at the site k k kkk at the moment t t ttt, then at the end of a time step δ t δ t delta t\delta tδt the particle will be either at the site k − 1 k − 1 k-1k-1k−1, or at k + 1 k + 1 k+1k+1k+1, with probabilities equal to 1 / 2 1 / 2 1//21 / 21/2. So in the time interval ( t , t + δ t t , t + δ t t,t+delta tt, t+\delta tt,t+δt ) the k k kkk particle moves with − δ x / δ t − δ x / δ t -delta x//delta t-\delta x / \delta t−δx/δt velocity, respectively δ x / δ t δ x / δ t delta x//delta t\delta x / \delta tδx/δt, to neighboring site. These rules define a "random walkers cellular automaton", similar to that proposed by Nishidate and Baba [1996], is done.
We denote by n k n k n_(k)n_{k}nk the number of particles at the site k k kkk at the moment t t ttt and by n k ′ n k ′ n_(k)^(')n_{k}^{\prime}nk′ (respectively n k ′ ′ n k ′ ′ n_(k)^('')n_{k}^{\prime \prime}nk′′ ) the number of particles which move at the moment t + δ t t + δ t t+delta tt+\delta tt+δt to the site k − 1 k − 1 k-1k-1k−1 (respectively k + 1 k + 1 k+1k+1k+1 ). Then the number of particles at the site k k kkk at the moment t + δ t t + δ t t+delta tt+\delta tt+δt can be written
(2.1) n k ( t + δ t ) = n k + 1 ′ + n k − 1 ′ ′ (2.1) n k ( t + δ t ) = n k + 1 ′ + n k − 1 ′ ′ {:(2.1)n_(k)(t+delta t)=n_(k+1)^(')+n_(k-1)^(''):}\begin{equation*} n_{k}(t+\delta t)=n_{k+1}^{\prime}+n_{k-1}^{\prime \prime} \tag{2.1} \end{equation*}(2.1)nk(t+δt)=nk+1′+nk−1′′
If we note
δ n k = 1 2 n k − n k ′ = n ′ ′ − 1 2 n k δ n k = 1 2 n k − n k ′ = n ′ ′ − 1 2 n k deltan_(k)=(1)/(2)n_(k)-n_(k)^(')=n^('')-(1)/(2)n_(k)\delta n_{k}=\frac{1}{2} n_{k}-n_{k}^{\prime}=n^{\prime \prime}-\frac{1}{2} n_{k}δnk=12nk−nk′=n′′−12nk
then, the number of particles at the site k k kkk and the moment t + δ t t + δ t t+delta tt+\delta tt+δt, (2.1), writes
(2.2) n k ( t + δ t ) = n k + 1 ′ + n k − 1 ′ ′ = 1 2 ( n k − 1 + n k + 1 ) − δ n k + 1 + δ n k − 1 (2.2) n k ( t + δ t ) = n k + 1 ′ + n k − 1 ′ ′ = 1 2 n k − 1 + n k + 1 − δ n k + 1 + δ n k − 1 {:(2.2)n_(k)(t+delta t)=n_(k+1)^(')+n_(k-1)^('')=(1)/(2)(n_(k-1)+n_(k+1))-deltan_(k+1)+deltan_(k-1):}\begin{equation*} n_{k}(t+\delta t)=n_{k+1}^{\prime}+n_{k-1}^{\prime \prime}=\frac{1}{2}\left(n_{k-1}+n_{k+1}\right)-\delta n_{k+1}+\delta n_{k-1} \tag{2.2} \end{equation*}(2.2)nk(t+δt)=nk+1′+nk−1′′=12(nk−1+nk+1)−δnk+1+δnk−1
For a large enough number of particles in each site, the last two terms in (2.2) may be neglected. Indeed, for a large enough n k n k n_(k)n_{k}nk, according to the law of large numbers, the frequencies n k ′ / n k = n k ′ ′ / n k n k ′ / n k = n k ′ ′ / n k n_(k)^(')//n_(k)=n_(k)^('')//n_(k)n_{k}^{\prime} / n_{k}=n_{k}^{\prime \prime} / n_{k}nk′/nk=nk′′/nk tend to the jump probability 1 / 2 1 / 2 1//21 / 21/2 and δ n k δ n k deltan_(k)\delta n_{k}δnk tends to zero. If it is assumed that the concentration c c ccc exists as a smooth enough function of x ∈ R x ∈ R x inRx \in \mathbb{R}x∈R and t ∈ R t ∈ R t inRt \in \mathbb{R}t∈R, then n k n k n_(k)n_{k}nk can be approximated by c ( x , t ) δ x c ( x , t ) δ x c(x,t)delta xc(x, t) \delta xc(x,t)δx, and (2.2) can be approximated by
(2.3) c ( x , t + δ t ) = 1 2 c ( x − δ x , t ) + 1 2 c ( x + δ x , t ) (2.3) c ( x , t + δ t ) = 1 2 c ( x − δ x , t ) + 1 2 c ( x + δ x , t ) {:(2.3)c(x","t+delta t)=(1)/(2)c(x-delta x","t)+(1)/(2)c(x+delta x","t):}\begin{equation*} c(x, t+\delta t)=\frac{1}{2} c(x-\delta x, t)+\frac{1}{2} c(x+\delta x, t) \tag{2.3} \end{equation*}(2.3)c(x,t+δt)=12c(x−δx,t)+12c(x+δx,t)
where x = k δ x x = k δ x x=k delta xx=k \delta xx=kδx. We have to emphasize that (2.3) is only an approximation of the exact relation (2.2). The difference between the two relations is a fluctuant quantity which, for a large n k n k n_(k)n_{k}nk, can be neglected in comparison with the dominant terms. Relation (2.3) can be written as
c ( x , t + δ t ) − c ( x , t ) δ t = δ x 2 2 δ t c ( x + δ x , t ) − 2 c ( x , t ) + c ( x − δ x , t ) δ x 2 c ( x , t + δ t ) − c ( x , t ) δ t = δ x 2 2 δ t c ( x + δ x , t ) − 2 c ( x , t ) + c ( x − δ x , t ) δ x 2 (c(x,t+delta t)-c(x,t))/(delta t)=(deltax^(2))/(2delta t)(c(x+delta x,t)-2c(x,t)+c(x-delta x,t))/(deltax^(2))\frac{c(x, t+\delta t)-c(x, t)}{\delta t}=\frac{\delta x^{2}}{2 \delta t} \frac{c(x+\delta x, t)-2 c(x, t)+c(x-\delta x, t)}{\delta x^{2}}c(x,t+δt)−c(x,t)δt=δx22δtc(x+δx,t)−2c(x,t)+c(x−δx,t)δx2
Which is the expression in finite differences of the one-dimensional diffusion equation, if the diffusion coefficient has the value given by
(2.4) D = δ x 2 2 δ t (2.4) D = δ x 2 2 δ t {:(2.4)D=(deltax^(2))/(2delta t):}\begin{equation*} D=\frac{\delta x^{2}}{2 \delta t} \tag{2.4} \end{equation*}(2.4)D=δx22δt
Since the concentration coincides with the one-particle probability density, c = p c = p c=pc=pc=p,
using (2.4), the relation (2.3) can also be expressed as
∂ t p ( x , t ) ≈ D δ x 2 [ p ( x + δ x , t ) − 2 p ( x , t ) + p ( x − δ x , t ] ∂ t p ( x , t ) ≈ D δ x 2 [ p ( x + δ x , t ) − 2 p ( x , t ) + p ( x − δ x , t ] del_(t)p(x,t)~~(D)/(deltax^(2))[p(x+delta x,t)-2p(x,t)+p(x-delta x,t]\partial_{t} p(x, t) \approx \frac{D}{\delta x^{2}}[p(x+\delta x, t)-2 p(x, t)+p(x-\delta x, t]∂tp(x,t)≈Dδx2[p(x+δx,t)−2p(x,t)+p(x−δx,t]
i. e. the master equation of the one-dimensional random walk, with the transition probability per time unity D / δ x 2 D / δ x 2 D//deltax^(2)D / \delta x^{2}D/δx2. One proves that for δ x → 0 δ x → 0 delta x rarr0\delta x \rightarrow 0δx→0 the master equation solution approximates the diffusion equation solution ([Gardiner, 1983, cap. 3, 7, van Kampen, 1981, cap 9, 10).
To verify that the algorithm described by (2.2) approximates the diffusion equation one considers the concentration of N N NNN non-interacting particles in Brownian motion,
(2.5) c ( x , t ) = N 4 π D t e − x 2 4 D t , x ∈ ( − ∞ , ∞ ) , t > 0 (2.5) c ( x , t ) = N 4 π D t e − x 2 4 D t , x ∈ ( − ∞ , ∞ ) , t > 0 {:(2.5)c(x","t)=(N)/(sqrt(4pi Dt))e^(-(x^(2))/(4Dt))","x in(-oo","oo)","t > 0:}\begin{equation*} c(x, t)=\frac{N}{\sqrt{4 \pi D t}} e^{-\frac{x^{2}}{4 D t}}, x \in(-\infty, \infty), t>0 \tag{2.5} \end{equation*}(2.5)c(x,t)=N4πDte−x24Dt,x∈(−∞,∞),t>0
To avoid the initial condition used in deriving (2.5), given by the Dirac function, c ( x , 0 ) = δ ( x ) c ( x , 0 ) = δ ( x ) c(x,0)=delta(x)c(x, 0)=\delta(x)c(x,0)=δ(x), we take as the initial time moment t 0 = δ t t 0 = δ t t_(0)=delta tt_{0}=\delta tt0=δt and the initial condition given by (2.5) for t = δ t t = δ t t=delta tt=\delta tt=δt. If the total number of particles is N N NNN, then one estimates
(2.6) n k ( t 0 ) = [ c ( x k , δ t ) δ x ] (2.6) n k t 0 = c x k , δ t δ x {:(2.6)n_(k)(t_(0))=[c(x_(k),delta t)delta x]:}\begin{equation*} n_{k}\left(t_{0}\right)=\left[c\left(x_{k}, \delta t\right) \delta x\right] \tag{2.6} \end{equation*}(2.6)nk(t0)=[c(xk,δt)δx]
where by [ r ] [ r ] [r][r][r] we note the integer part of the real number r r rrr and c c ccc is considered as given by (2.5). Because of the truncations, the actual number of particles N ∗ = ∑ k = − m m n k N ∗ = ∑ k = − m m   n k N^(**)=sum_(k=-m)^(m)n_(k)N^{*}=\sum_{k=-m}^{m} n_{k}N∗=∑k=−mmnk is smaller than N N NNN, but the difference is very small, N − N ∗ ≪ N N − N ∗ ≪ N N-N^(**)≪NN-N^{*} \ll NN−N∗≪N. To do the time step according to (2.2) one acts as it follows. For each particle i i iii laying in x k x k x_(k)x_{k}xk a random number q i ∈ { 0 , 1 } q i ∈ { 0 , 1 } q_(i)in{0,1}q_{i} \in\{0,1\}qi∈{0,1} is generated with the probability 1 / 2 1 / 2 1//21 / 21/2 for both values. One consider that, if q i = 0 q i = 0 q_(i)=0q_{i}=0qi=0 the particle i i iii moves to x k − 1 x k − 1 x_(k-1)x_{k-1}xk−1 and if q i = 1 q i = 1 q_(i)=1q_{i}=1qi=1 then it moves to x k + 1 x k + 1 x_(k+1)x_{k+1}xk+1. Repeating this operation for all N ∗ N ∗ N^(**)N^{*}N∗ particles, the new values n k ( t + δ t ) n k ( t + δ t ) n_(k)(t+delta t)n_{k}(t+\delta t)nk(t+δt) are obtained according to (2.2).
The probability density is the solution of the diffusion equation in an infinite domain. Since the numerical spatial mesh is bounded, it is possible to obtain a realistic simulation of the Brownian motion only if the concentration at the extremities of the mesh is negligible, i.e.
(2.7) c ( ± m δ x , t ) δ x ≪ N . (2.7) c ( ± m δ x , t ) δ x ≪ N . {:(2.7)c(+-m delta x","t)delta x≪N.:}\begin{equation*} c( \pm m \delta x, t) \delta x \ll N . \tag{2.7} \end{equation*}(2.7)c(±mδx,t)δx≪N.
This inequality must be satisfied both at the initial moment t = δ t t = δ t t=delta tt=\delta tt=δt and at the subsequent moments of the simulation, therefore the time t t ttt is limited by the condition (2.7). Thus the influence of the boundary condition will be negligible. To keep the total number of particles constant, one considers that the particles getting the extremities of the mesh remain in those points, i.e. (2.2) becomes
(2.8) n m ( t + δ t ) − n m = n m − 1 ′ ′ şi n − m ( t + δ t ) − n − m = n − m + 1 ′ . (2.8) n m ( t + δ t ) − n m = n m − 1 ′ ′  şi  n − m ( t + δ t ) − n − m = n − m + 1 ′ . {:(2.8)n_(m)(t+delta t)-n_(m)=n_(m-1)^('')quad" şi "quadn_(-m)(t+delta t)-n_(-m)=n_(-m+1)^(').:}\begin{equation*} n_{m}(t+\delta t)-n_{m}=n_{m-1}^{\prime \prime} \quad \text { şi } \quad n_{-m}(t+\delta t)-n_{-m}=n_{-m+1}^{\prime} . \tag{2.8} \end{equation*}(2.8)nm(t+δt)−nm=nm−1′′ şi n−m(t+δt)−n−m=n−m+1′.
For the proper simulations the spatial range [ − 1 , 1 ] [ − 1 , 1 ] [-1,1][-1,1][−1,1] has been chosen. Then, for a m m mmm given, δ x = 1 / m δ x = 1 / m delta x=1//m\delta x=1 / mδx=1/m. The diffusion coefficient has been fixed at D = 0.5 D = 0.5 D=0.5D=0.5D=0.5. From (2.4) it follows that the time steep has to be δ t = 1 / m 2 δ t = 1 / m 2 delta t=1//m^(2)\delta t=1 / m^{2}δt=1/m2. Simulations have been done for different values of m m mmm and N N NNN.
In this way, the cellular automaton governed by (2.2) can be use to perform a numerical computation of the one-dimensional diffusion equation, whose theoretical solution is (2.5). We note here that it is possible to solve numerically other partial differential equations in a similar manner. To achieve that one assign to continuous fields numbers of fictitious particles, similarly to (2.6), and one look for the corresponding cellular automata rules (for the diffusion equation the rule is expressed by (2.2)).
To obtain a quantitative hint on the numerical approximation for the solution of the diffusion equation, we have compared the theoretical probability to find a particle in the site k , c ( x k , t ) δ x / N k , c x k , t δ x / N k,c(x_(k),t)delta x//Nk, c\left(x_{k}, t\right) \delta x / Nk,c(xk,t)δx/N, where c c ccc is the theoretical solution (2.5), with the numerical one n k / N ∗ n k / N ∗ n_(k)//N^(**)n_{k} / N^{*}nk/N∗. We have defined the global indicator I I III, given by the sum on the whole mesh of the absolute values of the differences of the two quantities, for a given moment,
(2.9) I ( t ) = ∑ k = − m m | n k N ∗ − c ( x k , t ) δ x N | (2.9) I ( t ) = ∑ k = − m m   n k N ∗ − c x k , t δ x N {:(2.9)I(t)=sum_(k=-m)^(m)|(n_(k))/(N^(**))-(c(x_(k),t)delta x)/(N)|:}\begin{equation*} I(t)=\sum_{k=-m}^{m}\left|\frac{n_{k}}{N^{*}}-\frac{c\left(x_{k}, t\right) \delta x}{N}\right| \tag{2.9} \end{equation*}(2.9)I(t)=∑k=−mm|nkN∗−c(xk,t)δxN|
The nearer I I III is about zero, the better the numerical solution approximates the theoretical one. The maximum value of I I III is 2 . In the Table 1 the values of I I III for different values of m m mmm and N N NNN are given. One notices that, in agreement with the theoretical results, the larger m m mmm and N N NNN, the better the approximation. We have to emphasize that, for D D DDD fixed, relation (2.4) imposes the existence of a relationship between δ x δ x delta x\delta xδx and δ t δ t delta t\delta tδt, so to get the mesh denser implies the decrease of the time step.
Table 1. The maximum and minimum values in the first 10 time steps of the indicator I I III with respect to the spatial mesh dimension
N = 10 3 N = 10 3 N=10^(3)N=10^{3}N=103 N = 10 4 N = 10 4 N=10^(4)N=10^{4}N=104 N = 10 5 N = 10 5 N=10^(5)N=10^{5}N=105 N = 10 6 N = 10 6 N=10^(6)N=10^{6}N=106
m = 10 m = 10 m=10m=10m=10 0.062 0.031 0.022 0.024
0.102 0.112 0.106 0.104
m = 20 m = 20 m=20m=20m=20 0.027 0.024 0.010 0.013
0.133 0.052 0.013 0.017
m = 30 m = 30 m=30m=30m=30 0.053 0.019 0.009 0.003
0.129 0.038 0.015 0.007
m = 40 m = 40 m=40m=40m=40 0.062 0.028 0.008 0.003
0.137 0.039 0.015 0.005
N=10^(3) N=10^(4) N=10^(5) N=10^(6) m=10 0.062 0.031 0.022 0.024 0.102 0.112 0.106 0.104 m=20 0.027 0.024 0.010 0.013 0.133 0.052 0.013 0.017 m=30 0.053 0.019 0.009 0.003 0.129 0.038 0.015 0.007 m=40 0.062 0.028 0.008 0.003 0.137 0.039 0.015 0.005| | $N=10^{3}$ | $N=10^{4}$ | $N=10^{5}$ | $N=10^{6}$ | | :---: | :---: | :---: | :---: | :---: | | $m=10$ | 0.062 | 0.031 | 0.022 | 0.024 | | | 0.102 | 0.112 | 0.106 | 0.104 | | $m=20$ | 0.027 | 0.024 | 0.010 | 0.013 | | | 0.133 | 0.052 | 0.013 | 0.017 | | $m=30$ | 0.053 | 0.019 | 0.009 | 0.003 | | | 0.129 | 0.038 | 0.015 | 0.007 | | $m=40$ | 0.062 | 0.028 | 0.008 | 0.003 | | | 0.137 | 0.039 | 0.015 | 0.005 |
Also, using the discretized form of the Fick's law, J = − D d c / d x J = − D d c / d x J=-Ddc//dxJ=-D d c / d xJ=−Ddc/dx, it is calculated the numerical diffusion coefficient approximating the chosen coefficient D = 0.5 D = 0.5 D=0.5D=0.5D=0.5 of the theoretical solution. The particle flux between the spatial sites k k kkk and k + 1 k + 1 k+1k+1k+1 is expressed by J k = ( n k ′ ′ − n k + 1 ′ ) / δ t J k = n k ′ ′ − n k + 1 ′ / δ t J_(k)=(n_(k)^('')-n_(k+1)^('))//delta tJ_{k}=\left(n_{k}^{\prime \prime}-n_{k+1}^{\prime}\right) / \delta tJk=(nk′′−nk+1′)/δt. Since the number of particles varies very much from a site of the mesh to another, an averaged diffusion coefficient is calculated, obtained by summing discretized Fick's law over the whole mesh. The results are presented in Table 2.
Table 2. The minimum and maximum values in the first 10 steps of the numerical coefficient with respect to the spatial mesh dimension
N = 10 3 N = 10 3 N=10^(3)N=10^{3}N=103 N = 10 4 N = 10 4 N=10^(4)N=10^{4}N=104 N = 10 5 N = 10 5 N=10^(5)N=10^{5}N=105 N = 10 6 N = 10 6 N=10^(6)N=10^{6}N=106
m = 10 m = 10 m=10m=10m=10 0.4873 0.4596 0.4861 0.4953
0.6870 0.5548 0.5141 0.5044
m = 20 m = 20 m=20m=20m=20 0.3347 0.4693 0.4795 0.4987
0.7077 0.5247 0.5142 0.5018
m = 30 m = 30 m=30m=30m=30 0.5659 0.4087 0.4944 0.5012
0.6212 0.5203 0.5041 0.5035
m = 40 m = 40 m=40m=40m=40 0.5659 0.4087 0.4944 0.5012
0.6212 0.5203 0.5041 0.5035
N=10^(3) N=10^(4) N=10^(5) N=10^(6) m=10 0.4873 0.4596 0.4861 0.4953 0.6870 0.5548 0.5141 0.5044 m=20 0.3347 0.4693 0.4795 0.4987 0.7077 0.5247 0.5142 0.5018 m=30 0.5659 0.4087 0.4944 0.5012 0.6212 0.5203 0.5041 0.5035 m=40 0.5659 0.4087 0.4944 0.5012 0.6212 0.5203 0.5041 0.5035| | $N=10^{3}$ | $N=10^{4}$ | $N=10^{5}$ | $N=10^{6}$ | | :---: | :---: | :---: | :---: | :---: | | $m=10$ | 0.4873 | 0.4596 | 0.4861 | 0.4953 | | | 0.6870 | 0.5548 | 0.5141 | 0.5044 | | $m=20$ | 0.3347 | 0.4693 | 0.4795 | 0.4987 | | | 0.7077 | 0.5247 | 0.5142 | 0.5018 | | $m=30$ | 0.5659 | 0.4087 | 0.4944 | 0.5012 | | | 0.6212 | 0.5203 | 0.5041 | 0.5035 | | $m=40$ | 0.5659 | 0.4087 | 0.4944 | 0.5012 | | | 0.6212 | 0.5203 | 0.5041 | 0.5035 |

3 THE COARSE-GRAINED SPACE-TIME AVERAGING

The random walk hypothesis is a simple kind of microscopic evolution law of the type of those used in the molecular dynamic. There, the macroscopic evolution law of the particles ensamble, taken as a continuum in which the particles loose their individuality, is obtained by space-time averaging [Koplik and Banavar, 1995]. In comparison with other methods of molecular dynamics type, this approach is different by the use of the coarse-grained averages with particular properties introduced in [Vamoş et al., 1996a,b, and Suciu et al. 1998].
The macroscopic interpretation of the results is made by the space-time averages (1.3) as it follows. Consider N ∗ N ∗ N^(**)N^{*}N∗ particles which move on the straight line, over a time interval [ 0 , T ] [ 0 , T ] [0,T][0, T][0,T]. The motion of each particle i i iii consists of a sequence of jumps between the spatial mesh sites { x k = k δ x ∣ − m ◻ k ◻ m } x k = k δ x ∣ − m ◻ k ◻ m {x_(k)=k delta x∣-m◻k◻m}\left\{x_{k}=k \delta x \mid-m \square k \square m\right\}{xk=kδx∣−m◻k◻m} at the moments equal to the integer multiples of the time step δ t δ t delta t\delta tδt. Over the time interval δ t δ t delta t\delta tδt the particles move between two neighbour sites of the spatial mesh, at constant speed δ x / δ t δ x / δ t delta x//delta t\delta x / \delta tδx/δt. Each particle can be introduced or extracted from the spatial mesh, therefore we denote by t i + t i + t_(i)^(+)t_{i}^{+}ti+the moment of its appearance and by t i − t i − t_(i)^(-)t_{i}^{-}ti−the moment of its disappearance. The position of the particle i i iii is given by the function x i : [ t i + , t i − ] → [ x − m , x m ] ⊂ R x i : t i + , t i − → x − m , x m ⊂ R x_(i):[t_(i)^(+),t_(i)^(-)]rarr[x_(-m),x_(m)]subRx_{i}:\left[t_{i}^{+}, t_{i}^{-}\right] \rightarrow\left[x_{-m}, x_{m}\right] \subset \mathbb{R}xi:[ti+,ti−]→[x−m,xm]⊂R, which presents discontinuities when the particle changes the motion direction at one of the sites of the spatial mesh. The particle velocity ξ i = x ˙ i : [ t i + , t i − ] → { − δ x / δ t , δ x / δ t } ξ i = x ˙ i : t i + , t i − → { − δ x / δ t , δ x / δ t } xi_(i)=x^(˙)_(i):[t_(i)^(+),t_(i)^(-)]rarr{-delta x//delta t,delta x//delta t}\xi_{i}=\dot{x}_{i}:\left[t_{i}^{+}, t_{i}^{-}\right] \rightarrow\{-\delta x / \delta t, \delta x / \delta t\}ξi=x˙i:[ti+,ti−]→{−δx/δt,δx/δt} is a step function undergoing jumps with the amplitude 2 δ x / δ t 2 δ x / δ t 2delta x//delta t2 \delta x / \delta t2δx/δt when the motion direction changes. Thus, both position and velocities are piece-wise analytic functions, as required by the the continuous modeling method from [Vamos e al., 1996a,b].
Particularizing (1.3) for the one-dimensional case, for the real parameters a a aaa and τ < T / 2 τ < T / 2 tau < T//2\tau<T / 2τ<T/2, we define the coarse-grained average of the particles number at a point ( x , t ) ∈ R × ( τ , T − τ ) ( x , t ) ∈ R × ( τ , T − τ ) (x,t)inRxx(tau,T-tau)(x, t) \in \mathbf{R} \times(\tau, T-\tau)(x,t)∈R×(τ,T−τ) becomes
(3.1) ⟨ 1 ⟩ ( x , t ) = 1 4 τ a ∑ i = 1 N ∗ ∫ t − τ t + τ H + ( a − | x i ( t ′ ) − x | ) d t ′ (3.1) ⟨ 1 ⟩ ( x , t ) = 1 4 τ a ∑ i = 1 N ∗   ∫ t − τ t + τ   H + a − x i t ′ − x d t ′ {:(3.1)(:1:)(x","t)=(1)/(4tau a)sum_(i=1)^(N^(**))int_(t-tau)^(t+tau)H^(+)(a-|x_(i)(t^('))-x|)dt^('):}\begin{equation*} \langle 1\rangle(x, t)=\frac{1}{4 \tau a} \sum_{i=1}^{N^{*}} \int_{t-\tau}^{t+\tau} H^{+}\left(a-\left|x_{i}\left(t^{\prime}\right)-x\right|\right) d t^{\prime} \tag{3.1} \end{equation*}(3.1)⟨1⟩(x,t)=14τa∑i=1N∗∫t−τt+τH+(a−|xi(t′)−x|)dt′
In [Vamoş et al. 1997a] it is shown that (3.1) verifies the balance equation
(3.2) ∂ t ⟨ 1 ⟩ + ∂ x ⟨ ξ ⟩ = g (3.2) ∂ t ⟨ 1 ⟩ + ∂ x ⟨ ξ ⟩ = g {:(3.2)del_(t)(:1:)+del_(x)(:xi:)=g:}\begin{equation*} \partial_{t}\langle 1\rangle+\partial_{x}\langle\xi\rangle=g \tag{3.2} \end{equation*}(3.2)∂t⟨1⟩+∂x⟨ξ⟩=g
where
(3.3) ⟨ ξ ⟩ ( x , t ) = 1 4 τ a ∑ i = 1 N ∗ ∫ t − τ t + τ ξ i ( t ′ ) H + ( a − | x i ( t ′ ) − x | ) d t ′ (3.3) ⟨ ξ ⟩ ( x , t ) = 1 4 τ a ∑ i = 1 N ∗   ∫ t − τ t + τ   ξ i t ′ H + a − x i t ′ − x d t ′ {:(3.3)(:xi:)(x","t)=(1)/(4tau a)sum_(i=1)^(N^(**))int_(t-tau)^(t+tau)xi_(i)(t^('))H^(+)(a-|x_(i)(t^('))-x|)dt^('):}\begin{equation*} \langle\xi\rangle(x, t)=\frac{1}{4 \tau a} \sum_{i=1}^{N^{*}} \int_{t-\tau}^{t+\tau} \xi_{i}\left(t^{\prime}\right) H^{+}\left(a-\left|x_{i}\left(t^{\prime}\right)-x\right|\right) d t^{\prime} \tag{3.3} \end{equation*}(3.3)⟨ξ⟩(x,t)=14τa∑i=1N∗∫t−τt+τξi(t′)H+(a−|xi(t′)−x|)dt′
is the particles flux density and
(3.4) g ( x , t ) = 1 4 τ a ∑ i = 1 N ∗ ∫ t − τ t + τ [ H + ( a − | x i ( t i + ) − x | ) H + ( τ − | t i + − t | ) − − H + ( a − | x i ( t i − ) − x | ) H + ( τ − | t i − − t | ) ] (3.4) g ( x , t ) = 1 4 τ a ∑ i = 1 N ∗   ∫ t − τ t + τ   H + a − x i t i + − x H + τ − t i + − t − − H + a − x i t i − − x H + τ − t i − − t {:(3.4){:[g(x","t)=(1)/(4tau a)sum_(i=1)^(N^(**))int_(t-tau)^(t+tau)[H^(+)(a-|x_(i)(t_(i)^(+))-x|)H^(+)(tau-|t_(i)^(+)-t|)-:}],[{:-H^(+)(a-|x_(i)(t_(i)^(-))-x|)H^(+)(tau-|t_(i)^(-)-t|)]]:}:}\begin{array}{r} g(x, t)=\frac{1}{4 \tau a} \sum_{i=1}^{N^{*}} \int_{t-\tau}^{t+\tau}\left[H^{+}\left(a-\left|x_{i}\left(t_{i}^{+}\right)-x\right|\right) H^{+}\left(\tau-\left|t_{i}^{+}-t\right|\right)-\right. \tag{3.4}\\ \left.-H^{+}\left(a-\left|x_{i}\left(t_{i}^{-}\right)-x\right|\right) H^{+}\left(\tau-\left|t_{i}^{-}-t\right|\right)\right] \end{array}(3.4)g(x,t)=14τa∑i=1N∗∫t−τt+τ[H+(a−|xi(ti+)−x|)H+(τ−|ti+−t|)−−H+(a−|xi(ti−)−x|)H+(τ−|ti−−t|)]
describes the variation of the particles number in the spatial mesh. The functions ⟨ 1 ⟩ , ⟨ ξ ⟩ ⟨ 1 ⟩ , ⟨ ξ ⟩ (:1:),(:xi:)\langle 1\rangle,\langle\xi\rangle⟨1⟩,⟨ξ⟩ and g g ggg have a.e. continuous partial derivatives of the first order.
For 2 a = δ x 2 a = δ x 2a=delta x2 a=\delta x2a=δx and τ ≪ δ t / 2 T τ ≪ δ t / 2 T tau≪delta t//2T\tau \ll \delta t / 2 Tτ≪δt/2T, from (3.1) and (3.3), estimations are obtained for the concentration field and flux density, given by
(3.5) ⟨ 1 ⟩ ( x k , t l ) = c ( x k , t l ) (3.5) ⟨ 1 ⟩ x k , t l = c x k , t l {:(3.5)(:1:)(x_(k),t_(l))=c(x_(k),t_(l)):}\begin{equation*} \langle 1\rangle\left(x_{k}, t_{l}\right)=c\left(x_{k}, t_{l}\right) \tag{3.5} \end{equation*}(3.5)⟨1⟩(xk,tl)=c(xk,tl)
respectively,
(3.6) ⟨ ξ ⟩ ( x k + δ x / 2 , t l + δ t / 2 ) = J ( x k , t l ) (3.6) ⟨ ξ ⟩ x k + δ x / 2 , t l + δ t / 2 = J x k , t l {:(3.6)(:xi:)(x_(k)+delta x//2,t_(l)+delta t//2)=J(x_(k),t_(l)):}\begin{equation*} \langle\xi\rangle\left(x_{k}+\delta x / 2, t_{l}+\delta t / 2\right)=J\left(x_{k}, t_{l}\right) \tag{3.6} \end{equation*}(3.6)⟨ξ⟩(xk+δx/2,tl+δt/2)=J(xk,tl)
The improvement of the random walkers numerical algorithm, described in § 2 , consists in the possibility to calculate the terms (3.1), (3.3) and (3.4) of the balance equation (3.2) at the points x ≠ x k x ≠ x k x!=x_(k)x \neq x_{k}x≠xk and t ≠ t l t ≠ t l t!=t_(l)t \neq t_{l}t≠tl. Also, as discussed in section 1 , the time averaging reduces the number of particles needed to approximate the concentration field with a given precision and, implicitly, the computing time.

4 NUMERICAL MODELING OF DIFFUSION AT SMALL CONCENTRATIONS

The concentration field and the particles flux have been obtained by space-time coarse-grained averages, from the trajectories of the particles of the one-dimensional cellular automaton. There was obtained a good correspondence between the numerical results and the theoretical estimations. This proves that, under the hypothesis of the independent motion of particles, even at very small concentrations, Fick's law is still verified, if the space-time scale is big enough. In the following we present a few results which show the effectiveness of this numerical approach.
The coarse-grained average ⟨ 1 ⟩ ⟨ 1 ⟩ (:1:)\langle 1\rangle⟨1⟩, given by (3.1), is the sum of the contributions ⟨ 1 ⟩ i ⟨ 1 ⟩ i (:1:)_(i)\langle 1\rangle_{i}⟨1⟩i of the N N NNN particles from the mesh, during the averaging time interval: ⟨ 1 ⟩ = ∑ i = 1 N ⟨ 1 ⟩ i ⟨ 1 ⟩ = ∑ i = 1 N   ⟨ 1 ⟩ i (:1:)=sum_(i=1)^(N)(:1:)_(i)\langle 1\rangle=\sum_{i=1}^{N}\langle 1\rangle_{i}⟨1⟩=∑i=1N⟨1⟩i. The smooth concentration field is the average of ⟨ 1 ⟩ ⟨ 1 ⟩ (:1:)\langle 1\rangle⟨1⟩ over all particles evolutions, from the center of the mesh to the extremities, ⟨ 1 ⟩ ― = N ⟨ 1 ⟩ i ― ⟨ 1 ⟩ ¯ = N ⟨ 1 ⟩ i ¯ bar((:1:))=N bar((:1:)_(i))\overline{\langle 1\rangle}=N \overline{\langle 1\rangle_{i}}⟨1⟩―=N⟨1⟩i―. (This corresponds to a numerical simulation of the average with respect to a statistical ensemble.) The
corresponding dispersion is σ = N σ i σ = N σ i sigma=sqrtNsigma_(i)\sigma=\sqrt{N} \sigma_{i}σ=Nσi. The condition that the relative error of ⟨ 1 ⟩ ⟨ 1 ⟩ (:1:)\langle 1\rangle⟨1⟩ with respect to ⟨ 1 ⟩ ― ⟨ 1 ⟩ ¯ bar((:1:))\overline{\langle 1\rangle}⟨1⟩― be less than ε , 3 σ ε ⟨ 1 ⟩ ― ε , 3 σ ε ⟨ 1 ⟩ ¯ epsi,3sigmaquad epsi bar((:1:))\varepsilon, 3 \sigma \quad \varepsilon \overline{\langle 1\rangle}ε,3σε⟨1⟩―, with a reliance level of 0.997 , allows the estimation of the necessary number of particles, by
(4.1) 3 N σ i ◻ E N ⟨ 1 ⟩ i ― ⟹ N ≥ ( 3 σ i ε ⟨ 1 ⟩ i ― ) 2 (4.1) 3 N σ i ◻ E N ⟨ 1 ⟩ i ¯ ⟹ N ≥ 3 σ i ε ⟨ 1 ⟩ i ¯ 2 {:(4.1)3sqrtNsigma_(i)◻EN bar((:1:)_(i))Longrightarrow N >= ((3sigma_(i))/(epsi bar((:1:)_(i))))^(2):}\begin{equation*} 3 \sqrt{N} \sigma_{i} \square \mathcal{E} N \overline{\langle 1\rangle_{i}} \Longrightarrow N \geq\left(\frac{3 \sigma_{i}}{\varepsilon \overline{\langle 1\rangle_{i}}}\right)^{2} \tag{4.1} \end{equation*}(4.1)3Nσi◻EN⟨1⟩i―⟹N≥(3σiε⟨1⟩i―)2
Using (4.1), relationships between the particles number N N NNN, time scale 2 τ 2 τ 2tau2 \tau2τ and space scale a a aaa, needed for a given precision of the simulated concentration field, can be obtained [Vamoş et al., 1997b].
The Table 3 presents the time scale, measured in 10 5 δ t 10 5 δ t 10^(5)delta t10^{5} \delta t105δt units, needed to get the concentration at the mesh sites with a relative error ε = 0.01 ― ε = 0.01 _ epsi=0.01 _\varepsilon=\underline{0.01}ε=0.01―, as function of the total number of particles N N N\mathcal{N}N, used for the computation of ⟨ 1 ⟩ i ⟨ 1 ⟩ i (:1:)_(i)\langle 1\rangle_{i}⟨1⟩i and σ i σ i sigma_(i)\sigma_{i}σi, for a fixed space scale a = 0.125 δ x a = 0.125 δ x a=0.125 delta xa=0.125 \delta xa=0.125δx and N = 100000 N = 100000 N=100000N=100000N=100000.
Table 3.
N = 100 N = 100 N=100\mathcal{N}=100N=100 N = 200 N = 200 N=200\mathcal{N}=200N=200 N = 400 N = 400 N=400\mathcal{N}=400N=400 N = 800 N = 800 N=800\mathcal{N}=800N=800 N = 1600 N = 1600 N=1600\mathcal{N}=1600N=1600 N = 3200 N = 3200 N=3200\mathcal{N}=3200N=3200
x = 0.0 x = 0.0 x=0.0x=0.0x=0.0 4.31 4.98 4.74 4.96 4.83 4.71
x = 0.1 x = 0.1 x=0.1x=0.1x=0.1 4.84 5.39 5.17 5.36 5.11 5.02
x = 0.2 x = 0.2 x=0.2x=0.2x=0.2 5.47 6.11 5.90 6.09 5.76 5.65
x = 0.3 x = 0.3 x=0.3x=0.3x=0.3 5.98 6.99 6.75 6.90 6.47 6.43
x = 0.4 x = 0.4 x=0.4x=0.4x=0.4 6.85 8.13 7.87 7.95 7.50 7.53
x = 0.5 x = 0.5 x=0.5x=0.5x=0.5 8.22 9.37 9.12 9.28 9.05 9.11
x = 0.6 x = 0.6 x=0.6x=0.6x=0.6 10.9 11.2 11.1 11.4 11.4 11.5
x = 0.7 x = 0.7 x=0.7x=0.7x=0.7 15.3 15.0 14.9 15.3 15.2 15.2
x = 0.8 x = 0.8 x=0.8x=0.8x=0.8 22.0 22.3 22.5 22.9 22.2 22.1
x = 0.9 x = 0.9 x=0.9x=0.9x=0.9 46.0 44.3 44.3 45.1 43.6 44.0
N=100 N=200 N=400 N=800 N=1600 N=3200 x=0.0 4.31 4.98 4.74 4.96 4.83 4.71 x=0.1 4.84 5.39 5.17 5.36 5.11 5.02 x=0.2 5.47 6.11 5.90 6.09 5.76 5.65 x=0.3 5.98 6.99 6.75 6.90 6.47 6.43 x=0.4 6.85 8.13 7.87 7.95 7.50 7.53 x=0.5 8.22 9.37 9.12 9.28 9.05 9.11 x=0.6 10.9 11.2 11.1 11.4 11.4 11.5 x=0.7 15.3 15.0 14.9 15.3 15.2 15.2 x=0.8 22.0 22.3 22.5 22.9 22.2 22.1 x=0.9 46.0 44.3 44.3 45.1 43.6 44.0| | $\mathcal{N}=100$ | $\mathcal{N}=200$ | $\mathcal{N}=400$ | $\mathcal{N}=800$ | $\mathcal{N}=1600$ | $\mathcal{N}=3200$ | | :--- | :--- | :--- | :--- | :--- | :--- | :--- | | $x=0.0$ | 4.31 | 4.98 | 4.74 | 4.96 | 4.83 | 4.71 | | $x=0.1$ | 4.84 | 5.39 | 5.17 | 5.36 | 5.11 | 5.02 | | $x=0.2$ | 5.47 | 6.11 | 5.90 | 6.09 | 5.76 | 5.65 | | $x=0.3$ | 5.98 | 6.99 | 6.75 | 6.90 | 6.47 | 6.43 | | $x=0.4$ | 6.85 | 8.13 | 7.87 | 7.95 | 7.50 | 7.53 | | $x=0.5$ | 8.22 | 9.37 | 9.12 | 9.28 | 9.05 | 9.11 | | $x=0.6$ | 10.9 | 11.2 | 11.1 | 11.4 | 11.4 | 11.5 | | $x=0.7$ | 15.3 | 15.0 | 14.9 | 15.3 | 15.2 | 15.2 | | $x=0.8$ | 22.0 | 22.3 | 22.5 | 22.9 | 22.2 | 22.1 | | $x=0.9$ | 46.0 | 44.3 | 44.3 | 45.1 | 43.6 | 44.0 |
One finds that 200 particles are enough for the estimation of the time scale.
The time scale can be obtained for any space scale a a aaa and any point in ( − 1 , 1 − 1 , 1 -1,1-1,1−1,1 ). For instance, Table 4 presents the time scale, measured in 10 5 δ t 10 5 δ t 10^(5)delta t10^{5} \delta t105δt units, as function of x x xxx and a a aaa, for fixed N = 200 N = 200 N=200\mathcal{N}=200N=200 and ε = 0.01 ε = 0.01 epsi=0.01\varepsilon=0.01ε=0.01.
Table 4.
a = δ x / 8 a = δ x / 8 a=delta x//8a=\delta x / 8a=δx/8 a = δ x / 4 a = δ x / 4 a=delta x//4a=\delta x / 4a=δx/4 a = δ x / 2 a = δ x / 2 a=delta x//2a=\delta x / 2a=δx/2 a = δ x a = δ x a=delta xa=\delta xa=δx a = 2 δ x a = 2 δ x a=2delta xa=2 \delta xa=2δx a = 4 δ x a = 4 δ x a=4delta xa=4 \delta xa=4δx a = 8 δ x a = 8 δ x a=8delta xa=8 \delta xa=8δx
x = 0.000 x = 0.000 x=0.000x=0.000x=0.000 4.98 4.98 4.98 4.98 2.59 1.45 0.95
x = 0.025 x = 0.025 x=0.025x=0.025x=0.025 10.1 10.1 6.24 4.16 2.40 1.41 0.95
x = 0.050 x = 0.050 x=0.050x=0.050x=0.050 10.1 10.1 10.1 3.91 2.36 1.40 0.95
x = 0.075 x = 0.075 x=0.075x=0.075x=0.075 10.1 10.1 6.64 4.29 2.46 1.43 0.95
x = 0.100 x = 0.100 x=0.100x=0.100x=0.100 5.39 5.39 5.39 5.39 2.74 1.49 0.96
a=delta x//8 a=delta x//4 a=delta x//2 a=delta x a=2delta x a=4delta x a=8delta x x=0.000 4.98 4.98 4.98 4.98 2.59 1.45 0.95 x=0.025 10.1 10.1 6.24 4.16 2.40 1.41 0.95 x=0.050 10.1 10.1 10.1 3.91 2.36 1.40 0.95 x=0.075 10.1 10.1 6.64 4.29 2.46 1.43 0.95 x=0.100 5.39 5.39 5.39 5.39 2.74 1.49 0.96| | $a=\delta x / 8$ | $a=\delta x / 4$ | $a=\delta x / 2$ | $a=\delta x$ | $a=2 \delta x$ | $a=4 \delta x$ | $a=8 \delta x$ | | :---: | :---: | :---: | :---: | :---: | :---: | :---: | :---: | | $x=0.000$ | 4.98 | 4.98 | 4.98 | 4.98 | 2.59 | 1.45 | 0.95 | | $x=0.025$ | 10.1 | 10.1 | 6.24 | 4.16 | 2.40 | 1.41 | 0.95 | | $x=0.050$ | 10.1 | 10.1 | 10.1 | 3.91 | 2.36 | 1.40 | 0.95 | | $x=0.075$ | 10.1 | 10.1 | 6.64 | 4.29 | 2.46 | 1.43 | 0.95 | | $x=0.100$ | 5.39 | 5.39 | 5.39 | 5.39 | 2.74 | 1.49 | 0.96 |
As we can see, as space scale rises, the time scale decreases significantly, but only for a ≥ δ x a ≥ δ x a >= delta xa \geq \delta xa≥δx. For a ◻ δ x a ◻ δ x a◻delta xa \square \delta xa◻δx, one finds that between the sites the needed time scale is twice larger, which indicate a larger variability of coarse-grained averages.
In [Suciu et al., 1998], a numerical estimation of the diffusion coefficient was proposed, by D ∗ = [ ⟨ x ⟩ ― ⋅ ⟨ ξ ⟩ ― − ⟨ x ⋅ ξ ⟩ ― ] / c D ∗ = [ ⟨ x ⟩ ¯ ⋅ ⟨ ξ ⟩ ¯ − ⟨ x ⋅ ξ ⟩ ¯ ] / c D^(**)=[ bar((:x:))* bar((:xi:))- bar((:x*xi:))]//cD^{*}=[\overline{\langle x\rangle} \cdot \overline{\langle\xi\rangle}-\overline{\langle x \cdot \xi\rangle}] / cD∗=[⟨x⟩―⋅⟨ξ⟩―−⟨x⋅ξ⟩―]/c. To check this relation, we simulated the one-dimensional diffusion for 21100 particles, moving inside the space range [ − 1 , 1 ] [ − 1 , 1 ] [-1,1][-1,1][−1,1]. A variable space step mesh was used so that, according (2.4), D ( { k < 0 } ) = 0.05 D ( { k < 0 } ) = 0.05 D({k < 0})=0.05D(\{k<0\})=0.05D({k<0})=0.05 and D ( { k > 0 } ) = 5.0 D ( { k > 0 } ) = 5.0 D({k > 0})=5.0D(\{k>0\})=5.0D({k>0})=5.0. We also imposed that the particles undergo elastic collisions with the "walls" of the space interval. The coarse-grained averages were computed for a = 0.5 δ x a = 0.5 δ x a=0.5 delta xa=0.5 \delta xa=0.5δx and τ = 10 δ t τ = 10 δ t tau=10 delta t\tau=10 \delta tτ=10δt. For example, the Table 5 gives the values of the concentration c c ccc and the computed diffusion coefficients D ∗ D ∗ D^(**)D^{*}D∗, at the time moment t = 19 δ t t = 19 δ t t=19 delta tt=19 \delta tt=19δt.
Table 5. Numerically derived diffusion coefficient
x x xxx c × 10 3 c × 10 3 c xx10^(3)c \times 10^{3}c×103 D ∗ D ∗ D^(**)D^{*}D∗
-0.9 1.478 0.04997
-0.7 1.835 0.04942
-0.5 2.926 0.04842
-0.3 4.893 0.04837
-0.1 7.853 0.04885
0.0 4.885 0.62864
0.1 1.068 4.96512
0.3 1.068 4.96512
0.5 1.061 4.96478
0.7 1.054 4.49644
0.9 1.054 4.49644
x c xx10^(3) D^(**) -0.9 1.478 0.04997 -0.7 1.835 0.04942 -0.5 2.926 0.04842 -0.3 4.893 0.04837 -0.1 7.853 0.04885 0.0 4.885 0.62864 0.1 1.068 4.96512 0.3 1.068 4.96512 0.5 1.061 4.96478 0.7 1.054 4.49644 0.9 1.054 4.49644| $x$ | $c \times 10^{3}$ | $D^{*}$ | | ---: | ---: | ---: | | -0.9 | 1.478 | 0.04997 | | -0.7 | 1.835 | 0.04942 | | -0.5 | 2.926 | 0.04842 | | -0.3 | 4.893 | 0.04837 | | -0.1 | 7.853 | 0.04885 | | 0.0 | 4.885 | 0.62864 | | 0.1 | 1.068 | 4.96512 | | 0.3 | 1.068 | 4.96512 | | 0.5 | 1.061 | 4.96478 | | 0.7 | 1.054 | 4.49644 | | 0.9 | 1.054 | 4.49644 |
One finds that D ∗ D ∗ D^(**)D^{*}D∗ is a good estimation for the diffusion coefficients. Also, for deterministic and uniform motions of the particles, the correlation-like quantity is ⟨ x ⋅ ξ ⟩ ― / ⟨ x ⟩ ― ⋅ ⟨ ξ ⟩ ― ≡ 1 ⟨ x ⋅ ξ ⟩ ¯ / ⟨ x ⟩ ¯ ⋅ ⟨ ξ ⟩ ¯ ≡ 1 bar((:x*xi:))// bar((:x:))* bar((:xi:))-=1\overline{\langle x \cdot \xi\rangle} / \overline{\langle x\rangle} \cdot \overline{\langle\xi\rangle} \equiv 1⟨x⋅ξ⟩―/⟨x⟩―⋅⟨ξ⟩―≡1, and the coefficients D ∗ D ∗ D^(**)D^{*}D∗ identically vanish, as it is expected.

5 DIFUSION IN RANDOM FIELDS. THE MODEL OF MATHERON AND DE MARSILY

In a similar manner, we built a 2 -dimensional random walkers cellular automaton inside the rectangular mesh { x k = k δ x , y l = l δ y ; − m ◻ k , l ◻ m } x k = k δ x , y l = l δ y ; − m ◻ k , l ◻ m {x_(k)=k delta x,y_(l)=l delta y;-m◻k,l◻m}\left\{x_{k}=k \delta x, y_{l}=l \delta y ;-m \square k, l \square m\right\}{xk=kδx,yl=lδy;−m◻k,l◻m}. On it, we superposed an advection given by the samples of a random field and we computed ensemble averages. For horizontal advection and constant velocities in each layer ( k ; l , l + 1 ) ( k ; l , l + 1 ) (k;l,l+1)(k ; l, l+1)(k;l,l+1), we get the numerical simulation of the model proposed by Matheron and de Marsily [1980],
(5.1) d x ( ω V ) ( t ) = V ( y ( t ) , ω V ) d t + D d w ( t ) , d y ( t ) = D d w ( t ) (5.1) d x ω V ( t ) = V y ( t ) , ω V d t + D d w ( t ) , d y ( t ) = D d w ( t ) {:(5.1)dx^((omega_(V)))(t)=V(y(t),omega_(V))dt+Ddw(t)","quad dy(t)=Ddw(t):}\begin{equation*} d x^{\left(\omega_{V}\right)}(t)=V\left(y(t), \omega_{V}\right) d t+D d w(t), \quad d y(t)=D d w(t) \tag{5.1} \end{equation*}(5.1)dx(ωV)(t)=V(y(t),ωV)dt+Ddw(t),dy(t)=Ddw(t)
where w w www is the Wiener process, D D DDD is the diffusion coefficient, ω V ω V omega_(V)\omega_{V}ωV labels the realizations of the random field V V VVV. This model describes the transport of polutants in stratified natural porous media. Due to the heterogeneity of hydraulic conductivities (which through the Darcy law give the filatation velocities in porous media) the advection velocities look like the sampels of a random field and the two-dimensional random walk describes the mixing between horizontal layers. The authors' aim was to answer to the question when transport of polutants in such heterogeneous media can be described by a diffusion-like equation. Diffusive behavior corresponds to linear time dependence, σ x 2 ( t ) ∼ t σ x 2 ( t ) ∼ t sigma_(x)^(2)(t)∼t\sigma_{x}^{2}(t) \sim tσx2(t)∼t, of the longitudinal dispersion σ x 2 σ x 2 sigma_(x)^(2)\sigma_{x}^{2}σx2 (defined by mean square displacements) [van Kampen, 1981]. In their paper, Matheron and de Marsily show that situations occur when the transport is not diffusive. For instance, if the correlation function of horizontal advection velocities has a Gaussian shape ( ∼ e − y 2 ∼ e − y 2 ∼e^(-y^(2))\sim e^{-y^{2}}∼e−y2 ) the dispersion behaves as σ x 2 ∼ t 3 / 2 σ x 2 ∼ t 3 / 2 sigma_(x)^(2)∼t^(3//2)\sigma_{x}^{2} \sim t^{3 / 2}σx2∼t3/2, i. e. the the transport is superdiffusive.
Using the method from [Schwarze et al., 1998] we produced several random advection fields with different correlations (first picture). The second picture contains the corresponding dispersions curves given by numerical computations. Unbiased random walk ("ubrw") obvious has diffusive behavior. For Gaussian correlated longitudinal field we get nondiffusive behavior reported by Matheron and de Marsily [1980], σ x 2 ∼ t 3 / 2 σ x 2 ∼ t 3 / 2 sigma_(x)^(2)∼t^(3//2)\sigma_{x}^{2} \sim t^{3 / 2}σx2∼t3/2 ("ubrw+vx(y)") and, when a constant transversal velocity is added, diffusive long-time behavior ("ubrwv"). For random fields with identical values in a given number of neighbor layers (the last curve correspond to a 2-layers correlation), σ x 2 σ x 2 sigma_(x)^(2)\sigma_{x}^{2}σx2 goes to ∼ t 2 ∼ t 2 ∼t^(2)\sim t^{2}∼t2 dependence for increasing correlation length.

Concentration fields for the simulated situations can be estimated by coarsegrained averages, similar to (3.1). When the values of the advection velocities and diffusion coefficients are associated with an advection-diffusion equation like in [Tompson and Gelhar, 1990] the coarse-grained averages give the numerical solution of the partial differential equation. Thus our numerical method is very similar

N. SUCIU, C. VAMOŞ, U. JAEKEL, H. VEREECKEN

to the "particle tracking" method of Tompson and Gelhar (the main difference is that the particles move in a grid not in continuous space). But this method is also useful in analyzing more complex diffusion problems in random fields (like those discussed in [Matheron and de Marsily, 1980] when no transpot equations are available). In such cases the existence of balance equations of form (3.2) could be a powerful numerical tool.

REFERENCES

1.Gardiner, C. W. (1983)
Handbook of Stochastic Methods (for Physics, Chemistry and Natural Science), Springer-
Verlag, New York
2.Koplik, J. and J.R. Banavar (1995)
Continuum deductions from molecular hydrodynamics, Annu. Rev. Fluid Mech., 27, 257-292.
3.Landau, L.D. and E.M. Lifchitz (1988)
Statistical Physics (in romanian), Ed. Tehnică, Bucureşti.
4.Matheron,G., and G. de Marsily (1980)
Is Transport in porous media always diffusive?, Water Resour. Res., 16, 901-917. 5.Nishidate, K., M. Baba and R. J. Gaylord (1996)
Cellular automaton model for random walkers, Phys. Rev. Lett., 77 (9), 1675-78. 6.Rumshinski, L.Z. (1974)
Mathematical Processing of Experimental Data (in romanian), Ed. Tehnică, Bucureşti.
7.Schwarze, H., U. Jaekel and H. Vereecken (1998)
Estimation of macrodispersion by different approximation methods for flow and transport in
randomly heterogeneous media (to be published).
8.Suciu, N., C. Vamoş, A. Georgescu, U. Jaekel and H. Vereecken (1998)
Trnasport processes in porous media. 1. Continuous modeling (this issue).
9.Tompson, A. F. B. and L. W. Gelhar (1990)
Numerical simulation in three-dimensional randomly heterogeneous media, Water Resour. Res.,
26 (10), 2541-2562.
10.Vamoş, C., A. Georgescu, N. Suciu and I. Turcu (1996a)
Balace equations for physical systems with corpuscular structure, Physica A, 227, 81-92.
11.Vamoş, C., A. Georgescu and N. Suciu (1996b)
Balance equations for a finite number of particles, St. Cerc. Mat., 48, 115-127.
12.Vamoş, C., N. Suciu and A. Georgescu(1997a)
Hydrodynamic equations for one-dimensional systems of inelastic particles, Phys. Rev. E, 55,
6277-6280.
13.Vamoş, C., N. Suciu and M. Peculea (1997b)
Numerical modelling of the one-dimensional diffusion by random walkers, Rev. Anal. Numér.
Théorie Approximation, 26(1-2), 237-247.
14.van Kampen, N. G. (1981)
Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam.
1998

Related Posts