Particle Diffusion

Particle Diffusion Distribution

Particle Analog of Diffusion

Having already solved the diffusion equation

`(deln)/(delt)=alpha(del^2n)/(delx^2)\ \ \ \ \ (0)`

for analog temperature or density diffusion, this will be the discrete particle version of the same physical behavior. Here we will concentrate on demonstrating the diffusion of particles and matching their computed density to the density expected from the diffusion equation. We've already demonstrated particle diffusion where the starting condition is with particles of type 1 on the left side of the container and type 2 on the right side. Here, instead, we will have a narrow swath of "dopant" particles centered in the container and immersed in a sea of host particles. Assuming that the initial dopant particle very time dependence of the dopant particle density, `n`, should then follow the expression (see Appendix 0) Assuming that the initial width of the dopant particle swath is very narrow, time dependence of the dopant particle density, `n`, should then follow the expression

`n(x,t)=sqrt(1/(4pialphat))exp(-x^2/(4alphat))\ \ \ \ (1)`

where `t` is time and `alpha` is the diffusivity of the combined dopant and host particle seas. So the program objective is to fit the actual spread of the dopant particle distribution to equation 1 by computing an appropriate value for `alpha` using least squares fitting of the dopant density to the time varying density of the dopant.

Description of the Animation

The diffusing particles (dopant spheres) are initially located in a very narrow band (swath) at the horizontal (x) center of the container. The blue dopant spheres start with three dimensional (3D) random velocity directions. The dopant spheres collide with other spheres of their own family as well as a much larger family of red host spheres. The blue dopant spheres progress away from the center band is by random walk and therefore follows the laws of diffusion. The x density numerical distribution (histogram) of the blue dop[ant pheres should form a gaussian as expressed by equation 1. The program fits this numerical distribution to an analog gaussian by minimizing the least squared deviation of the analog gaussian relative to the numerical density distribution. From the `sigma` of the fitted gaussian, the numerically derived diffusivity is computed.

This diffusivity is to be compared with the expected diffusivity (see Appendix 3):

`alpha=lambda/3 sqrt((v^2)_(avg)) `

where `(v^2)_(avg)=1/Nsum_1^N(v_i^2)` and `lambda` is the mean free path and is to be computed from

`lambda=1/(n*4pi r^2) `

where N is the total number of host+dopant spheres per unit volume and r is the radius of the both sphere colors.

Finding Values of the Diffusivity `alpha`

It is not unusual to use least squares fitting to find the coefficient for an exponential like

`f(x)=a exp(-x/lambda)`

where `lambda` is the parameter that we need to find. For equation 1, however, we need to find the corresponding value for

`n(x)=n_(max) exp(-x^2/sigma^2)`

So we need to use the arrays of data

`NA_i=n_i`
and
`XA_i=x_i^2`

The least squares method then chooses both `n_(max)` and `sigma^2` that minimize

`sum_1^Nn_i^2-[n_(max)exp(-x_i^2/(sigma^2))]^2`

and from `sigma^2` we can compute `alpha` as

`alpha=sigma^2/(4t)`

where time is the time elapsed since the start button was pushed.

Results for Diffusivity `alpha`

Typically the values found for `alpha differ from the expected value by about 30%. Some of this error is due to the fact that equation 1 is for 1 dimensional (1D) particle motion and our particles have 3D motion. In addition equation 1 assumes that the initial dopant swath width is infinitesimal but we have had to use very finite swath width.

For the computation of hard sphere velocities after collision see Appendix 1 below the Canvas. Also see Appendix 2 for a discussion of the rationale for equation 0. For the derivation of equation 1 see Appendix 0.

Appendix 0: Deriving Equation 1 from Equation 0

Since equation 1 assumes an infinite span in the x direction, the usual separation of variables method will not work. Instead we can use Fourier transforms to handle the second derivative with respect to x:

`-alpha t k^2+ ik(x)`
so: `a=-alphat`
and: `b=i(x)`
and the term to subtract from the completed square is: `b^2/(4a)=-(x^2)/(4alphat))`
Then the exponent can be written:
`-alpha t k^2+ ik(x)= -alpha t (k-i(x)/(2alphat))^2-(x)^2/(4alphat)`

`1/(2pi)int_-oo^ooexp(-k^2alphat)exp(ikx)dk=sqrt(1/(4pialphat))e^(-x^2/(4alphat))\ \ \ \ \ (A0.5)`

Appendix 1: Deriving Expressions for the Final Velocities after Collision

For hard sphere collisions, we need to use both momentum and energy conservation to find the final velocities of particles. The change of momentum of one particle is the negative of the change of momentum of the other:

`((m_1bb"v"_1+Mdeltavbbhat"u")cdot (m_1bb"v"_1+Mdeltavbbhat"u"))/m_1+ ((m_2bb"v"_2-Mdeltavbbhat"u")cdot (m_2bb"v"_2-Mdeltavbbhat"u"))/m_2= m_1bb"v"_1cdotbb"v"_1+m_2bb"v"_2cdotbb"v"_2`

`2bb"v"_1cdot(Mdeltav bbhat"u")-2bb"v"_2cdot(Mdeltav bbhat"u")+ (Mdeltav)^2/m_1+ (Mdeltav)^2/m_2=0`

Appendix 2: Discussion of Rationale for Diffusion Equation 0

Isaac Newton expressed the equation for heat flow across a very thin layer as the difference between temperatures from one side of the layer to the other divided by the layer thickness:

If the gradient at a certain level of the thickness is reduced, then the temperature or concentration beyond that level can increase. This is where the diffusion equation with its second derivative with respect to x comes into play. Let's consider a bell shaped (gaussian) density distribution function: At the top the second derivative is negative so equation 0 says that the density will be declining there. At the left or right sides, the second derivative is positive so, according to equation 0, the density will be increasing there and, considering that all diffusion is the result of random walk, that is what we expect. So in conclusion we can just say that equation 0 is equation A2-1 converted to non-constant gradients.