Tutorial 4

Simulate solute spreading in shear flows. Reference: https://doi.org/10.1063/1.3078518.

For the commented code start here

Advection-Diffusion Equation

\[ \dot{u} + div (u \mathbf{V}) = \dot{u} + \mathbf{V} \cdot \nabla u + u\ \mathrm{div} (\mathbf{V}) = 0 \]

\[ D_t u = - u\ \mathrm{div} (\mathbf{V})\]

\[ u \mathbf{V} = \mathbf{J} = \mathbf{J^{\mathrm{drift}} + J^{\mathrm{diff}}} =\boldsymbol{\beta}\ u - D \nabla u \]

Concentration on the grid

In MPM, information about \(u\) is concentrated on the particles.

\[ U = \sum_p U_p\ \delta (x - x_p) \]

to transfer this quantity on the grid we apply a Particle-to-Grid (P2G) transfer.

\begin{align} M_i u_i &= \sum_j u_j \int_\Omega \varphi_j(x) \ \varphi_i(x) \\ &= \int_\Omega u(x)\ \varphi_i(x) \\ &= \sum_p \int_\Omega U_p\ \delta (x - x_p) \ \varphi_i(x) \\ &= \sum_p U_p\ \varphi_i(x_p)\quad \forall i \end{align}

Convective Flux

\[ \mathbf{J}^{\mathrm{drift}} = \boldsymbol{\beta}\ u \]

To transfer this quantity on the grid we use a P2G, one for each component of the flux, for instance

\[ J^{\mathrm{drift}}_{x} = \mathbf{J}^{\mathrm{drift}} \cdot \hat{x} \]

mathematically P2G can be described as

\begin{align} M_i \mathbf{j}^{\mathrm{drift}}_i \cdot \hat{x} &= \sum_k \int_\Omega \mathbf{j}^{\mathrm{drift}}_k \cdot \hat{x}\ \varphi_k(x)\ \varphi_i(x) \\ &= \int_\Omega \mathbf{J}^{\mathrm{drift}} \cdot \hat{x}\ \varphi_i(x)\\ &= \sum_p \int_\Omega \boldsymbol{\beta}\cdot\hat{x}\ U_p\ \delta (x - x_p)\ \varphi_i(x)\\ &= \sum_p \boldsymbol{\beta}(x_p)\cdot\hat{x}\ U_p\ \varphi_i(x_p)\qquad \forall i \end{align}

Diffusive Flux

\[\mathbf{J}^{\mathrm{diff}} = -D\ \nabla u\]

To transfer this quantity on the grid we use a Particle-to-Grid-Derivative (P2GD) transfer, one for each component of the flux, for instanceo

\[ {J}^{\mathrm{diff}}_{x} = \mathbf{J}^{\mathrm{diff}} \cdot \hat{x} \]

mathematically the P2GD can be described as

\begin{align} M_i \mathbf{j}^{\mathrm{diff}}_i \cdot \hat{x} &= \sum_k \int_\Omega \mathbf{j}^{\mathrm{diff}}_k \cdot \hat{x}\ \varphi_k(x)\ \varphi_i(x) \\ &= \int_\Omega \mathbf{J}^{\mathrm{diff}} \cdot \hat{x}\ \varphi_i(x)\\ &= - \int_\Omega D\ \partial_x u\ \varphi_i(x)\\ &= \int_\Omega D\ u\ \partial_x \varphi_i(x) - \int_{\partial\Omega} D\ \partial_x u\ \varphi_i(x)\\ &= \sum_p \int_\Omega D\ U_p\ \delta (x - x_p)\ \partial_x \varphi_i(x)\\ &= \sum_p D \ U_p\ \partial_x \varphi_i(x_p)\qquad \forall i \end{align}

Velocity Update

\[ M_i \mathbf{v}_i = (M_i \mathbf{j}^{\mathrm{drift}}_i + M_i \mathbf{j}^{\mathrm{diff}}_i) \]

\[ \mathbf{v}^{\mathrm{drift}}_i = (M_i \mathbf{j}^{\mathrm{drift}}_i) / (M_i u_i) \]

\[ \mathbf{v}^{\mathrm{diff}}_i = (M_i \mathbf{j}^{\mathrm{diff}}_i) / (M_i u_i) \]

We transfer the velocity back to the particles using a Grid-to-Particles (G2P) transfer.*

\[ \mathbf{V}_p = \sum_i \left(\mathbf{v}^{\mathrm{drift}}_i\ + \ \mathbf{v}^{\mathrm{diff}}_i\right)\varphi_i(x_p) \]