Format And Type Fonts CCHHEEMMIICCAALL EENNGGIINNEEEERRIINNGG TTRRAANNSSAACCTTIIOONNSS VOL. 39, 2014 A publication of The Italian Association of Chemical Engineering www.aidic.it/cet Guest Editors: Petar Sabev Varbanov, Jiří Jaromír Klemeš, Peng Yen Liew, Jun Yow Yong Copyright © 2014, AIDIC Servizi S.r.l., ISBN 978-88-95608-30-3; ISSN 2283-9216 DOI:10.3303/CET1439124 Please cite this article as: Reverberi A.P., Ferrando R., Fabiano B., Dovì V.G., 2014, Discrete modelling of a multiparticle diffusion process and front propagation: dynamics and source effects, Chemical Engineering Transactions, 39, 739-744 DOI:10.3303/CET1439124 739 Discrete Modelling of a Multiparticle Diffusion Process and Front Propagation: Dynamics and Source Effects Andrea P. Reverberi *a , Riccardo Ferrando b , Bruno Fabiano c , Vincenzo G. Dovì a a DCCI - Department of Chemistry and Industrial Chemistry, via Dodecaneso 31, 16146 Genova (ITALY) b DIFI - Department of Physics, via Dodecaneso 33, 16146 Genova (ITALY) c DICCA - Department of Civil, Chemical and Environmental Engineering – Chemical Engineering Section, via Opera Pia 15, 16145 Genova (ITALY) reverb@dichep.unige.it A multiparticle diffusion process is simulated in strip geometry. The pointlike particles are generated in a zone at constant concentration and allowed to diffuse in a semi-infinite medium with periodic boundaries. The average particle position, its concentration trend and the particle ensemble width are investigated at short and long times. Systematic deviations are observed with respect to the Fickian behaviour at the beginning of the diffusion process where a transient superdiffusive trend appears, despite the absence of particles interactions. Besides, the classical continuum approach shows a persistent lack of fit as far as the concentration at the diffusion front is concerned. The results are related with previous studies concerning collective on-lattice diffusion in a concentration gradient. 1. Introduction Modelling of diffusion processes is basically important for its implications in many fields of research, namely in heat diffusion (Solisio et al., 2012), in chemically reacting systems (Fabiano et al., 2012) and in combustion propagation (Palazzi and Fabiano, 2012). As a schematic classification, the aforementioned modelling can be divided in mesoscale and microscale approaches. Mesoscale simulation techniques can be adopted in the study of particle dynamics (Xu et al., 2011), in atmospheric gas-aerosol transport (Aloyan et al., 2010), in modelling diffusion through composite materials (Pascariu et al., 2013) and in surface-adsorbate catalytic processes. The basic idea of mesoscopic models relies upon applying continuum equations derived from the microscopics on larger cells with respect to the dimension of the single molecule. As a parallel approach, lattice models of diffusion processes on microscopic scale were object of intensive work in discrete or continuous time. In particles diffusion simulations, it was pointed out that the presence of repeated reflections may trigger an anomalous (superdiffusive) behaviour in the particle square displacement with time. Montecarlo (MC) simulations in modelling diffusion processes are traditionally related to the development of diffusion fronts (Sapoval et al., 2005) in close connection with topologically anomalous geometries of percolative structures and with random packings of solid agglomerates (Bertei et al., 2013). A common feature in Montecarlo approaches is the use of on-lattice techniques which prove to be useful from a computational point of view, but have the drawback of adding finite size effects that can be crucial at short times (Chappa and Albano, 2004). Semi-lattice and off-lattice simulations represent a good choice to overcome the previously mentioned disadvantages; however their use in modelling diffusion processes is limited to a small number of cases mainly related to random or correlated motion in trapping media. In this paper, which represents the corresponding off-lattice version of Kolb et al. (1987) scheme, we propose a simulation describing a collective motion of pointlike particles in strip geometry, where a generator kept at a fixed particle concentration behaves as a source of random walkers undergoing diffusion. The paper is divided as follows: in section 2, we outline the computational details of the present algorithm. In section 3, we discuss the results with particular attention to the average width <(t)> of the mailto:reverb@dichep.unige.it 740 0-  • • • • •• • • • •• • ••• • • • • p • • • • • • • • • • • • • • • • • • • •• m x y source 0-  • • • • •• • • • •• • ••• • • • • p • • • • • • • • • • • • • • • • • • • •• m x y source Figure 1: Scheme of the off-lattice diffusion process here considered particle ensemble, its average position  )t(x and to the asymptotic concentration  )]t(x[c . Finally, in section 4 we draw the conclusions and we trace the direction for future works. 2. The model Simulations are realized according to the scheme visualized in Figure 1. All particles are allowed to move in a semi-infinite rectangular region of width m; their x and y coordinates are stored in an array where they are ranked according to the distance from the origin. The algorithm consists of the following steps: a) At time t=0, N0 points are randomly located in the area spanning in the range - (dashed lines) and < x > (solid lines) versus the number of Monte Carlo steps. The upper two lines refer to p=510 -4 ; the lower ones to p=510 -5 . Dotted lines are tangent to the curves for long times and have slope =0.5. N0=256 ;  =0.05 ; m=1 2 0t CMCM 0 )0(R)t(R Ntd2 1 LimD;D S 1 D   (3) where DCM is the centre of mass diffusion coefficient,  i i r)t(R is the collective coordinate, <> is the average operation on the number of realizations and S0 is the thermodynamic factor defined as:    20 )N( N S  (4) S0 is related to the particle number fluctuation in a simulation box embedded in the d-dimensional space where diffusion occurs. As expected, we find a quadratic dependence D~p 2 consistent with the one observed in on-lattice simulations (Reed and Ehrlich, 1981). From this point on, we will consider the dynamics of the collective motion. In Figure 2, the average position    N 1i i x N 1 )t(x and the width          21 1 21 / N i i )]t(xx[ N )t( of N(t) particles are plotted versus the number of Monte Carlo steps, linearly depending on time, for different values of the jump length p. As expected, z/1 ht)t(  and z/1 kt)t(x  with 1/z=0.5 asymptotically, namely for t, according to a Brownian motion. However, the trend for short times is dominated by a superdiffusive behaviour which persists at longer times for lower values of the diffusion coefficient. Analogous tendency was observed in an anomalous spreading of a density front from an infinite continuous source in a concentration-dependent lattice gas. Küntz and Lavallée (2003) used a two dimensional lattice gas automaton and observed that the slope of the cumulative absorption curve was significantly larger than 0.5 for a wide time range during the simulation. Here, the situation is even more intriguing as we observe a transient superdiffusive trend in the absence of interactions, owing to the constancy of the jump length. The plot of Figure 3 is analogous to the previous one, with the difference that, in this case, the particle concentration in the source is varied while the diffusion coefficient remains constant. Again, we observe an anomalous superdiffusive trend which is more pronounced for lower concentration values within the generator. We know that an enhanced probability of long-scale motion (Lévy flights) and memory effects conditioning the following steps of a walker may trigger a superdiffusive behaviour (Lacasta et al., 2004). Heuristically, we believe that, at least in the first phases of the motion when the return probability is negligible, the random walker is subject to a source effect which tends to disappear when the layers close to the generator attain a saturation. 742 Figure 3: Log-log plot of < > (dashed lines) and < x > (solid lines) versus the number of Monte Carlo steps. All lines are plotted for p=10 -4 .  and m as in Figure 2 In the context of gradient multiparticle diffusion, the role of the classical Fick II law that reads: 2 2 x c D t c      (5) with the following initial and boundary conditions:         1 0 1 c)t,(c c)t,0(c c)0,0x(c 01 cc  (6) was deeply discussed in literature and some inadequacies were pointed out for short timescales (McGahay, 2004). In the following, we will show that a mismatch apparently limited to early times leads to further and more evident shortcomings for long times. Eq(5), with the boundary condition Eq(6), has an analytical solution that reads:        Dt2 x erf)cc(c)t,x(c 010 (7) Owing to its parabolic character, the previous equation does not show a sharply cut advancing front typical of a propagatory hyperbolic problem, where shock waves allows one to intuitively identify the front as the zone of the moving discontinuity of the diffusion profile. Nevertheless, we define the abscissa s(t) of the diffusion front as follows: 01 c c cc dcx )t(s 1 0    (8) In parabolic diffusion equations, a suitable change of variables is generally adopted according to: Dt2 x  (9) with such a change, Eq(7) is differentiated at constant time and it is replaced in Eq(8), obtaining the following expression for the abscissa s(t) of the moving front:     0 de/Dt4)t(s 2   = /Dt2 (10) Finally, the concentration c[s(t)] on the moving front is determined by replacing Eq(10) in Eq(7), which, for c1=0, assumes the form:   )](erf1[c)]t(s[c 1 0  0.4249 c0 = constant (11) 743 Figure 4: Plot of u(t) versus the number of Monte Carlo steps for p=10 -4 and N0=1024 (solid line) ; p=10 -4 and N0=512 (dashed line) ; p=10 -3 and N0=1024 (dotted line).  and m as in Figure 2 We remember that c[s(t)]/c0 in a continuum approach corresponds to the scaled concentration gen c/)]t(x[c)t(u  at the average position of the N1 particles ensemble in the discrete approach. However, we have found a significant lack of fit between Eq(11) and the results of our discrete simulations. To prove this statement, we have plotted, in Figure 4, the dimensionless concentration at the average position u(t) versus time for different values of particles N0 contained in the generator and for different values of p. We note that u(t) is constant in time at u0.5320.005, irrespective of p and N0 . Kolb et al. (1987) proposed a lattice-gas model of diffusing particles generated at the left of a source line. The particles could move only on empty sites according to the energies of the nearest-neighbour sites. The authors reported the particles concentration versus the scaled variable  x/x and investigated the effects of attractive, null and repulsive interactions on the concentration at discrete average distances from the source. However, they overlooked the fact that, in the absence of interactions, the concentration at the average position was significantly different from the one predicted by the classical approach related to Eq(6-7). We stress that the asymptotically constant value of u obtained in the present simulations is very close to the one reported in Kolb et al. (1987) study for null interactions. Hence, both their on-lattice model and the present one describing a discrete multiparticle diffusion process show a persistent discrepancy for long times with respect to the corresponding continuum approach. A possible origin of the mismatch reported in Figure 4 can be ascribed to an unphysical limitation related to the phenomenological Fick’s law that reads: cDJ  (12) A discontinuity between initial and boundary conditions, as the one appearing in Eq(6), gives an unlimited local flux J at the start and a consequent infinite propagation speed x/t for t0. Despite this intrinsic weakness, typical of Eq(5-6), we stress that the classical approach still keeps to be utilized in many technical problems, such as in heat transfer (Reverberi et al., 2013), in the diffusion of polluting agents in the environment (Tagliabue et al., 2014) and in modelling fire propagation in explosion disasters (Palazzi et al., 2014). To circumvent the aforementioned drawbacks, Eq(5) can be replaced by the hyperbolic non- Fickian (wave-type) heat/mass transfer equation, generally known as the Maxwell-Cattaneo equation, containing a relaxation time  accounting for a finite speed of signal propagation, namely:   2 2 t c cD t c       (13) Despite Eq(13) is a valid tool to fill the aforementioned gaps of the traditional scheme, the instabilities related to the numerical solution of this equation in two or three spatial dimensions limit its use in real applications. 744 4. Conclusions The most important results can be summarized as follows: - A transient superdiffusive behaviour, whose duration depends on the collective diffusion coefficient and on the particle concentration in the generator, is observed despite the absence of particles interactions. - For long times, the concentration at the average position gen c/)]t(x[c)t(u  is asymptotically constant but it differs from the value predicted by the classical Fick II law, giving a value less than 20.1% of the one obtained by the present simulations. We believe that this is a long-lasting effect of the initial superdiffusive behaviour. As a real-world application, this result may be useful to describe the moving front concentration for alloy elements diffusion in metals. Acknowledgements The financial support from the EC FP7 project EFENIS (ENER/FP7/296003) is gratefully acknowledged. References Ala-Nissila T., Ferrando R., Ying S.C., 2002, Collective and single particle diffusion on surfaces, Advances in Physics, 51(3), 949-1078. Aloyan A.E., Arutyunyan V.O., Yermakov A.N., 2010, Regional-scale numerical modeling of gas-aerosol dynamics, Chemical Engineering Transactions, 22, 173-178, DOI: 10.3303/CET1022028 Bertei A., Nucci B., Nicolella C., 2013, Effective transport properties in random packings of spheres and agglomerates, Chemical Engineering Transactions, 32, 1531-1536, DOI: 10.33032/CET1332256 Chappa V.C., Albano E.V., 2004, A finite-size dynamic-scaling approach for the diffusion front of particles, Journal of Chemical Physics 121, (1), 328-332. Fabiano B., Reverberi A.P., Del Borghi A., Dovi V.G., 2012, Biodiesel production via transesterification: Process safety insights from kinetic modeling, Theoretical Foundations of Chemical Engineering, 46, (6), 673-680, DOI: 10.1134/S0040579512060097 Kolb M., Gouyet J.F., Sapoval B., 1987, Diffusion of interacting particles in a concentration gradient: scaling, critical slowing down and phase separation, Europhysics Letters, 3(1), 33-38. Küntz M., Lavallée P., 2003, Anomalous spreading of a density front from an infinite continuous source in a concentration-dependent lattice gas automaton diffusion model, Journal of Physics D: Applied Physics, 36(9), 1135-1142. Lacasta A.M., Sancho J.M., Romero A.H., Sokolov I.M., Lindenberg K, 2004, From subdiffusion to superdiffusion of particles on solid surfaces, Physical Review E 70(5), 051104-1-051104-10. McGahay V., Inertial effects and diffusion, 2004, Journal of Non-Crystalline Solids, 349(1-3), 234-241. Palazzi E., Fabiano B., 2012, Analytical modelling of hydrocarbon pool fires: Conservative evaluation of flame temperature and thermal power, Process Safety and Environmental Protection, 90(2), 121-128, DOI: 10.1016/j.psep.2011.06.009 Palazzi E., Currò F., Reverberi A., Fabiano B., 2014, Development of a theoretical framework for the evaluation of risk connected to accidental oxygen releases, Process Safety and Environmental Protection, DOI: 10.1016/j.psep.2014.02.015 (in press) Pascariu V., Avadanei O., Gasner P., Stoica I., Reverberi A.P., Mitoseriu L., 2013, Preparation and characterization of PbTiO 3-epoxy resin compositionally graded thick films, Phase Transitions, 86(7), 715-725, DOI: 10.1080/01411594.2012.726727 Reed D.A., Ehrlich G., 1981, Surface diffusion, atomic jump rates and thermodynamics, Surface Science, 102(2-3), 588-609. Reverberi A.P., Fabiano B., Dovì V.G., 2013, Use of inverse modelling techniques for the estimation of heat transfer coefficients to fluids in cylindrical conduits, International Communications in Heat and Mass Transfer, 42, 25-31, DOI: 10.1016/j.icheatmasstransfer.2012.12.005 Sapoval B., Andrade J.S. Jr., Baldassarri A., Desolneux A., Devreux F., Filoche M., Grebenkov D., Russ S., 2005, New simple properties of a few irregular systems, Physica A, 357(1), 1-17. Solisio C., Reverberi A.P., Del Borghi A., Dovi' V.G., 2012, Inverse estimation of temperature profiles in landfills using heat recovery fluids measurements, Journal of Applied Mathematics, Volume 2012, Article number 747410, DOI: 10.1155/2012/747410 Tagliabue M., Reverberi A.P., Bagatin R., 2014, Boron removal from water: needs, challenges and perspectives, Journal of Cleaner Production, DOI: 10.1016/j.jclepro.2013.11.040 Xu J-B., Zhang S-F., Wu H., Zhao Y-H., Wen H., 2011, Mesoscopic simulation of aggregate structure and stability of heavy crude oil by GPU accelerated DPD, Chemical Engineering Transactions, 24, 1531- 1536, DOI: 10.3303/CET1124256