Team:Paris Saclay/Modeling

Modeling

The diffusion of particules is based on Fick's first law which is given by: \begin{equation} \textbf{j} = -D \, \bf{\nabla} n \label{eq:fick} \end{equation} In this equation, $\textbf{j}$ is the diffusion flux, $D$ the diffusion coefficient and $n$ the concentration of particles. This equation can be coupled with the continuity equation $\partial_t n = \mathbf{\nabla} \cdot \mathrm{j} \quad (+\sigma)$ expressing the conservation of the total number of diffusing particles. $\sigma$ is the net particle production rate. The beads being spherical, it is more interesting to work with spherical coordinates. The Laplace operator is then defined by : \begin{equation} \triangle a = \frac{1}{r^2} \frac{\partial}{\partial r} \Big(r^2 \frac{\partial a}{\partial r} \Big) + \frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \Big(\sin \theta \frac{\partial a}{\partial \theta} \Big) + \frac{1}{r^2 \sin \theta} \frac{\partial^2 \varphi}{\partial \varphi^2} \end{equation} Let's suppose that there is no dependence on angles in the beads, i.e. there is a spherical symmetry. We can write $n(r,\theta,\varphi,t) = n(r,t)$.

Steady diffusion

In steady diffusion, the equation $(2)$ is simpler as there is no dependence in time. \begin{equation} \frac{1}{r^2} \frac{\partial}{\partial r} \Big[ r^2 \frac{\partial n}{\partial r} \Big] =0 \end{equation} The equation above leads to the following differential equation : \begin{equation} 2r\frac{\mathrm{d}n}{\mathrm{d}r} + r^2 \frac{\mathrm{d}^2n}{\mathrm{d}r^2} = 0 \end{equation} Let's consider a variable m defined by $m=\frac{\mathrm{d}n}{\mathrm{d}r}$. The equation $(4)$ can be rewritten $ 2rm + r^2 \frac{\mathrm{d}m}{\mathrm{d}t}= 2 r +r^2 \frac{\mathrm{d}(\log m)}{\mathrm{d}r} =0$. We can show that the differential equation is then $m\frac{\alpha}{r^2}=\frac{\mathrm{d}n}{\mathrm{d}r}$. The concentration of particles is then : \begin{equation} n(r) = -\frac{\alpha}{r} + \beta \end{equation} $\alpha$ and $\beta$ being constants to be determined. One of the boundary conditions is determined by the density of the external medium, when the distance from the center of the bead is far enough compared to the size of a bead. When $r\rightarrow + \infty \, n=n_\infty $. The concentration equation then is : \begin{equation} n(r) = -\frac{\alpha}{r} + n_\infty$ \end{equation} The difficult part here is the determination of the the constant \alpha. There is a singularity when $r=0$. Another boundary condition can be determined by knowing the flux of particules crossing a given surface.

Unsteady diffusion

In the case of unsteady diffusion, there is a dependence in time of the diffusion equation $\partial_t n= \frac{D}{r^2} \, \partial_r(r^2 \partial_r n)$. Due to the conditions of decrease of $n$ in time, can be written as: \begin{equation} \frac{1}{D} \frac{\partial n}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial t} \Big(r^2 \frac{\partial n}{\partial r}\Big) = -\lambda^2 \label{eq:unsteady_lambda} \end{equation} If we consider $n$ as being the product of a spatial function and a temporal function, $n(r,t)=T(t)*R(r)$, the diffusion equation can be rewritten as : \begin{equation} \frac{1}{D} R(r) \frac{\mathrm{d}T(t)}{\mathrm{d}t} = \frac{1}{r^2} \frac{\mathrm{d}}{\mathrm{d}r} \Big[ r^2 T(t) \frac{\mathrm{d}R(r)}{\mathrm{d}r}\Big] \end{equation} or \begin{equation} \frac{1}{DT} \frac{\mathrm{d}T}{\mathrm{d}t} = \frac{1}{r^2 R} \frac{\mathrm{d}}{\mathrm{d}r} \Big[ r^2 \frac{\mathrm{d}R(r)}{\mathrm{d}r}\Big] \end{equation}

Temporal part

If we consider only the temporal part, we have: \begin{equation} \frac{1}{DT} \frac{\mathrm{d}T}{\mathrm{d}t} = - \lambda^2 \end{equation} This leads us to exponential dependence of the temporal part of the evolution of concentration in time: \begin{equation} T(t) = \mathrm{e}^{-\lambda^2 D t} \end{equation}

Radial part

Let's consider now the radial part of the diffusion equation: \begin{equation} \frac{1}{r^2 R} \frac{\mathrm{d}}{\mathrm{d}r} \Big[ r^2 \frac{\mathrm{d}R(r)}{\mathrm{d}r}\Big] = -\lambda^2 \end{equation} \begin{equation} \rightarrow r^2 \frac{\mathrm{d}^2R}{\mathrm{d}r^2} + 2 r \frac{\mathrm{d}R}{\mathrm{d}r} + \lambda^2 r^2 R =0 \, \end{equation} This is a Bessel differential equation. To solve it, we can use Newmann's method. The general formula of the differential equation is: \begin{equation} x^2 y''+ (2p+1)xy' + (\alpha^2x^{2r} + \beta^2)y = 0 \, \, \, \mathrm{with} \, y=y(x) \end{equation} Solution of the differential equation is : \begin{equation} y=x^{-p}[C_1 J_{^q/_r}(\frac{\alpha}{r}x^r)+ C_2 Y_{^q/_r}(\frac{\alpha}{r}x^r)] \end{equation} with $q=\sqrt{p^2-\beta^2}$ Here we have $2p+1 =2 \rightarrow p=1/2$ $\alpha^2x^{2r} + \beta^2 = \lambda^2 x^2 \rightarrow \beta=0, r=1, \alpha = \lambda$ The relations we had before then becomes : \begin{equation} R(r) = \frac{1}{\sqrt{r}} [C_1 J_{^1/_2}(\lambda r)+ C_2 Y_{^1/_2}(\lambda r)] \end{equation} $J_{^1/_2}$ and $Y_{^1/_2}$ are respectively Bessel functions of the first and second kind \begin{equation} \rightarrow \qquad n(r,t)=\frac{\mathrm{e}^{-\lambda D^2 t}}{\sqrt{r}} \Big[C_1 J_{^1/_2}(\lambda r)+ C_2 Y_{^1/_2}(\lambda r) \Big] \end{equation} $J_{^1/_2}(x) = \sqrt{\frac{2}{\pi x}} \sin x \qquad$ and $Y_{^1/_2}(x)=-J_{^{-1}/_2}(x) = \sqrt{\frac{2}{\pi x}} \cos x]$ Finally, the concentration is given by : \begin{equation} n(r,t) = \frac{\mathrm{e}^{-\lambda D^2 t}}{\sqrt{r}} \, \sqrt{\frac{2}{\lambda \pi}} \Big[C_1 \sin(\lambda r) - C_2 \cos(\lambda r) \Big] \end{equation}