Difference between revisions of "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 $. ---------------------------------------------------- \ref{steady_eq} $\rightarrow n(+\infty) = \beta = n_\infty \rightarrow n(r) = -\frac{\alpha}{r} + n_\infty$ Problem : determination of $\alpha$. Quelles conditions aux limites choisir pour une solution physique ? Unsteady diffusion : $\partial_t n= \frac{D}{r^2} \, \partial_r(r^2 \partial_r n)$ \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} it is necessarily negative,necessary condition for the decrease of $n$ in time On pose $n(r,t)=T(t)*R(r)$ \ref{eq:continuity_eq} $\rightarrow$ $\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]$ $\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]$ Temporal part \ref{eq:unsteady_lambda} $\rightarrow$ $\frac{1}{DT} \frac{\mathrm{d}T}{\mathrm{d}t} = - \lambda^2$ $ \rightarrow$ $\frac{DT}{T} = - \lambda^2 D \mathrm{d}t \rightarrow T(t) = \mathrm{e}^{-\lambda^2 D t}$ Radial part \ref{eq:unsteady_lambda} $\rightarrow$ $\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$ $\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 \,$ Eq diff de Bessel -> résolution avec méthode de Newmann General formula of the differential equation: \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 $(5) +(6) \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]$ $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]$ $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]$