# Angiography Simulation and Planning Using a Multi-Fluid Approach
## 3 Diffusion Simulation of Contrast Media
### Multi-Fluids Mixture Model  
Our multi-fluid flows mixture model of the blood-contrast media is essentially based on a SPH fluid simulation algorithm. The state of each particle \(i\) is computed at every time step by solving following fluid dynamics equations [23]. \[(\frac{\partial\mathbf{u}_{k}}{\partial t}+(\mathbf{u}_{m}\cdot\
abla) \mathbf{a}_{k})_{i}=-(\alpha_{k}\
abla\cdot\mathbf{u}_{m})_{i}-(\
abla\cdot( \alpha_{k}\mathbf{u}_{mk}))_{i} \tag{1}\]  
\[(\frac{\partial\mathbf{u}_{m}}{\partial t}+(\mathbf{u}_{m}\cdot\
abla) \mathbf{u}_{m})_{i}=-(\frac{\
abla p_{m}}{\rho_{m}})_{i}+\mathbf{g}+(\frac{ \
abla\cdot\mathbf{T}_{m}}{\rho_{m}})_{i}+(\frac{\
abla\cdot\mathbf{T}_{m}}{ \rho_{m}})_{i}+(\frac{\mathbf{F}_{c}}{\rho_{m}})_{i} \tag{2}\]  
where, \(\alpha_{k}\) is the volume fraction, \(\mathbf{u}_{m}\) is the mixture velocity, \(\mathbf{u}_{mk}\) is the drift velocity, \(\rho_{m}=\sum_{k}\alpha_{k}\rho_{k}\) is the mixture density, \(p_{m}\) is the mixture pressure, \(\mathbf{g}\) is the gravity, \(\mathbf{T}_{m}\) is the viscous stress tensor of mixture, \(\mathbf{T}_{Dm}\) is the convective momentum transfer between phases, and \(\mathbf{F}_{c}\) is the external force. The left-hand of Equation (1) is the substantial derivative of the volume fraction of phase \(k\) with regard to the mixture velocity \(\mathbf{u}_{m}\), and the left-hand of Equation (2) represents the substantial derivative of the mixture velocity \(\mathbf{u}_{m}\). Therefore, we solve these two equations by computing the right-hand of the equations with an explicit time integration. The drift velocity \(\mathbf{u}_{mk}\) is introduced to describe the change of the volume fraction of each phase, and the phase \(k\) of particle \(i\) of \(\mathbf{u}_{mk}\) expression is:  
\[(\mathbf{u}_{mk})_{i}=\tau(\rho_{k}-\sum_{k^{\prime}}\epsilon_{k^{\prime}i} \rho_{k^{\prime}})\mathbf{a}_{i}-\tau(\
abla p_{ki}-\sum_{k^{\prime}} \epsilon_{k^{\prime}i}\
abla p_{k^{\prime}i})-\sigma(\frac{\
abla\alpha_{ki} }{\alpha_{ki}}-\sum_{k^{\prime}}\epsilon_{k^{\prime}i}\frac{\
abla\alpha_{ki^ {\prime}i}}{\alpha_{k^{\prime}i}}) \tag{3}\]  
where the coefficient \(\tau\) is used to control the inertia and pressure effects, a higher value will cause stronger inertia and pressure effects, \(\sigma\) is set to control the diffusion effects, the higher the value is, the stronger the diffusion effect. \(\rho_{k}\) is the rest density of phase \(k\), \(\epsilon_{ki}=\frac{\alpha_{ki}\rho_{ki}}{\rho_{mi}}\) is the mass fraction of phase \(k\) and the acceleration \(\mathbf{a}_{i}\) is:  
\[\mathbf{a}_{i}=\mathbf{g}-((\mathbf{u}_{m}\cdot\
abla)\mathbf{u}_{m}-\frac{ \partial\mathbf{u}_{m}}{\partial t})_{i} \tag{4}\]  
which represents the difference between the gravity and the substantial derivative of the mixture velocity. To deal with the miscible fluid simulation, the miscible phase can move within each other by different pressures. Thus, the pressure of phase \(k\) of particle \(i\) in Equation (3) is expressed as:  
\[p_{ki}=\alpha_{ki}p_{mi} \tag{5}\]  
\[p_{mi}=\kappa(\overline{\rho}_{i}-\rho_{mi}) \tag{6}\]  
\[\overline{\rho}_{i}=\sum_{j}m_{j}\
abla W_{ij} \tag{7}\]  
where \(\overline{\rho}_{i}\) is the interpolated density over all mixture particles under a smoothing kernel function \(W_{ij}\), we adopt a poly6 kernel function for the density interpolation and a spiky kernel function for all other quantity interpolations in Muller [24], \(\kappa\) is the gas constant, \(m_{j}\) is the mass of particle \(j\). Before solving the governing equation, the drift velocity of the mixture fluid needs to be determined. For each particle \(i\), we have \[(\
abla p_{k})_{i}=\sum_{j}\frac{m_{j}}{\overline{\rho}_{j}}(p_{kj}-p_{ki})\
abla W _{ij} \tag{8}\]  
\[(\
abla\alpha_{k})_{i}=\sum_{j}\frac{m_{j}}{\overline{\rho}_{j}}(\alpha_{kj}- \alpha_{ki})\
abla W_{ij} \tag{9}\]  
where \(\
abla W_{ij}\) is the gradient of the smoothing kernel function. Therefore, the drift velocity \(\mathbf{u}_{mk}\) can be computed from Equation (3) by using \(\
abla p_{k}\) and \(\
abla\alpha_{k}\).