Computational geometric methods for preferential clustering of particle suspensions

TL;DR

This paper introduces a geometric numerical method for simulating suspensions of particles with improved accuracy and efficiency, preserving key physical properties and better capturing particle distribution dynamics.

Contribution

The paper presents a novel divergence-free radial basis function interpolation combined with a splitting integrator for particle suspension simulation, enhancing accuracy and computational cost-effectiveness.

Findings

More accurate particle distributions compared to conventional methods

Preserves the sum of the Lyapunov spectrum and mimics centrifuge effects

Effective for large suspensions of 10^4 particles in flow fields

Abstract

A geometric numerical method for simulating suspensions of spherical and non-spherical particles with Stokes drag is proposed. The method combines divergence-free matrix-valued radial basis function interpolation of the fluid velocity field with a splitting method integrator that preserves the sum of the Lyapunov spectrum while mimicking the centrifuge effect of the exact solution. We discuss how breaking the divergence-free condition in the interpolation step can erroneously affect how the volume of the particulate phase evolves under numerical methods. The methods are tested on suspensions of $10^4$ particles evolving in discrete cellular flow field. The results are that the proposed geometric methods generate more accurate and cost-effective particle distributions compared to conventional methods.

Tables

Table 1. Table 1 : The terms in the series expansion ( 75 ) for the physical volume Ψ [ n ] superscript Ψ delimited-[] 𝑛 \Psi^{[n]} under various numerical methods. Note that Ψ 1 [ n ] = 0 subscript superscript Ψ delimited-[] 𝑛 1 0 \Psi^{[n]}_{1}=0 for all the methods.

Method	${\Psi }_2^{\left[n\right]}$	${\Psi }_3^{\left[n\right]}$	${\Psi }_4^{\left[n\right]}$
Exact solution	$\frac{\alpha }{2}{\chi }_3$	$\frac{\alpha }{6}\left({\chi }_5-\alpha {\chi }_3\right)$	$\frac{\alpha }{24}({\chi }_1+\alpha {\chi }_2+{\alpha }^2{\chi }_3+3\alpha {\chi }_3^2$ $-2\alpha {\chi }_4-2\alpha {\chi }_5)$
FE1	0	0	0
LT1	$\alpha {\chi }_3$	$\frac{{\alpha }^2}{2}{\chi }_3$	$-\frac{{\alpha }^2}{6}\left(3{\chi }_4-\alpha {\chi }_3-3\chi _3^{}{}_{}{}^2\right)$
CP1	$\frac{\alpha \sqrt{6}}{6}{\chi }_3$	$\frac{\alpha \left(\sqrt{6}-1\right)}{6}{\chi }_5-\frac{{\alpha }^2\sqrt{6}}{12}{\chi }_3$	$\frac{\alpha }{36}(({\alpha }^2{\chi }_3-3\alpha {\chi }_5+\frac{7}{2}{\chi }_1)\sqrt{6}$ $+(3\chi _3^{}{}_{}{}^2-3{\chi }_4+3{\chi }_5)\alpha -6{\chi }_1)$
AB2	$\frac{3\alpha }{16\theta }{\chi }_3$	$\frac{3\alpha }{32}\left({\chi }_5-\alpha {\chi }_3\right)$	$\frac{\alpha }{128}\left(3\theta {\chi }_1+16\alpha {\chi }_3^2-16\alpha {\chi }_4\right)$
RK2	$\frac{\alpha }{2}{\chi }_3$	0	$\frac{{\alpha }^2}{8}\left({\chi }_3^2-{\chi }_4\right)$
CP2	$\frac{\alpha }{2}{\chi }_3$	$\frac{3{\alpha }^2}{16}{\chi }_3+\frac{\alpha }{6}{\chi }_5$	$\frac{\alpha }{576}(33{\alpha }^2{\chi }_3+72\alpha \chi _3^{}{}_{}{}^2+24\alpha {\chi }_2$ $-48\alpha {\chi }_4-60\alpha {\chi }_5+32{\chi }_1)$
SS2	$\frac{\alpha }{2}{\chi }_3$	$\frac{\alpha }{4}\left({\chi }_5-\alpha {\chi }_3\right)$	$\frac{\alpha }{48}\left(3{\chi }_1+4{\alpha }^2{\chi }_3+6\alpha {\chi }_3^2-6\alpha {\chi }_4-6\alpha {\chi }_5\right)$

Table 2. Table 2 : The errors of the physical volume after one time step for the numerical methods under consideration.

	$\left|\Delta {\Psi }^{\left[n\right]}\right|$
Method	Case (a) ($\nabla \cdot 𝐮\left(𝐱\right)\ne 0$)	Case (b) ($\nabla \cdot 𝐮\left(𝐱\right)=0$)	Case (c) ($𝐮\left(𝐱\right)={𝐮}_e\left(𝐱\right)$)
FE1	$\frac{\alpha }{2}{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{12}{h}^4\left|{\chi }_4\right|$	$\frac{{\alpha }^2}{12}{h}^4\left|{\chi }_4\right|$
LT1	$\alpha h^2\left|{\chi }_3\right|$	$\frac{5{\alpha }^2}{12}{h}^4\left|{\chi }_4+\frac{{\delta }_4}{2}\right|$	$\frac{5{\alpha }^2}{12}{h}^4\left|{\chi }_4\right|$
CP1	$\frac{\alpha \sqrt{6}}{6}{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{12}{h}^4\left|{\delta }_4\right|$	$O\left(h^5\right)$
AB2	$\frac{3\alpha }{16\theta }{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{24}{h}^4\left|{\chi }_4+\frac{{\delta }_4}{8}\right|$	$\frac{{\alpha }^2{h}^2}{24}{h}^4\left|{\chi }_4\right|$
RK2	$\frac{\alpha }{2}{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{24}{h}^4\left|{\chi }_4+\frac{{\delta }_4}{8}\right|$	$\frac{{\alpha }^2}{24}{h}^4\left|{\chi }_4\right|$
CP2	$\frac{\alpha }{2}{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{12}{h}^4\left|{\delta }_4\right|$	$O\left(h^5\right)$
SS2	$\frac{\alpha }{2}{h}^2\left|{\chi }_3\right|$	$\frac{{\alpha }^2}{24}{h}^4\left|{\chi }_4+\frac{{\delta }_4}{8}\right|$	$\frac{{\alpha }^2}{24}{h}^4\left|{\chi }_4\right|$

Table 3. Table 3 : Summary of the properties of the integration methods under consideration

	FE1	LT1	CP1	AB2	RK2	SS2	CP2
Order	1	1	1	2	2	2	2
Contractivity-preserving	No	Yes	Yes	No	No	Yes	Yes
Centrifuge-preserving	No	No	Yes	No	No	No	Yes

Table 4. Table 4 : The wall clock times for 10 4 superscript 10 4 10^{4} time steps using different integration methods and exact interpolation.

	FE1	LT1	CP1	AB2	RK2	CP2	SS2
$T_w\left(s\right)$	1.1817	1.3610	1.6577	2.0524	2.0376	2.5001	1.6455

Table 5. Table 6 : The relative entropy E ​ ( P n ) 𝐸 subscript 𝑃 𝑛 E(P_{n}) , first Wasserstein distance W ​ ( P n ) 𝑊 subscript 𝑃 𝑛 W(P_{n}) and Δ ​ 𝐱 n ¯ ¯ Δ subscript 𝐱 𝑛 \overline{\Delta\mathbf{x}_{n}} between the numerical distribution P n subscript 𝑃 𝑛 P_{n} and the reference distribution. The numerical distributions are calculated by various integration methods that use exact interpolation as shown in the second column. The first column contains the Stokes number S ​ t 𝑆 𝑡 St , time step h ℎ h and simulation time T 𝑇 T used in the six simulations. The first row contains the aspect ratio λ 𝜆 \lambda of the particle shape. Values with a − - mean that the numerical solution is unstable.

		$\lambda =1$			$\lambda =1/10$
	$P_n$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$
$St=\frac{1}{5}$	exact+FE1	7.0143	0.4466	0.5534	–	–	–
	exact+LT1	0.4596	0.0070	0.0098	2.3929	0.1793	0.2086
$h=\frac{1}{50}$	exact+CP1	0.1212	0.0049	0.0065	2.3100	0.1750	0.2048
	exact+AB2	4.0952	0.0697	0.0585	–	–	–
$T=4$	exact+RK2	1.6791	0.0234	0.0215	–	–	–
	exact+CP2	0.0581	0.0022	0.0027	0.5339	0.0666	0.0983
	exact+SS2	0.1169	0.0048	0.0062	2.2956	0.1739	0.2039
$St=1$	exact+FE1	6.9916	1.6321	2.5224	0.9086	0.6445	1.7383
	exact+LT1	0.1538	0.0929	1.0000	0.6645	0.4788	1.2536
$h=\frac{1}{20}$	exact+CP1	0.1061	0.1034	1.0084	0.5716	0.4293	1.2554
	exact+AB2	0.4833	0.2715	1.4977	–	–	–
$T=8$	exact+RK2	0.2696	0.1566	1.3083	1.3583	0.4631	1.4790
	exact+CP2	0.0786	0.0558	0.5512	0.2505	0.2892	1.0857
	exact+SS2	0.1435	0.0890	1.0044	0.6374	0.4858	1.2535
$St=10$	exact+FE1	6.8070	2.7326	2.8528	1.3145	1.9214	4.8392
	exact+LT1	0.3531	0.1626	0.6588	0.1002	0.5384	5.1653
$h=\frac{1}{5}$	exact+CP1	0.3014	0.1619	0.6479	0.0658	0.5721	5.1883
	exact+AB2	0.7651	0.3095	0.8671	0.8826	1.9070	4.9994
$T=20$	exact+RK2	0.5804	0.2522	0.9074	0.1593	0.7891	5.2612
	exact+CP2	0.0733	0.0667	0.2919	0.0711	0.2761	4.9678
	exact+SS2	0.3157	0.1567	0.6700	0.0992	0.5433	5.2206

Table 6. Table 7 : The relative entropy E ​ ( P n ) 𝐸 subscript 𝑃 𝑛 E(P_{n}) , first Wasserstein distance W ​ ( P n ) 𝑊 subscript 𝑃 𝑛 W(P_{n}) and average error per particle Δ ​ 𝐱 n ¯ ¯ Δ subscript 𝐱 𝑛 \overline{\Delta\mathbf{x}_{n}} between the numerical distribution P n subscript 𝑃 𝑛 P_{n} and the reference distribution. The numerical distributions are calculated by various interpolation methods that use CP2 integration as shown in the second column. The first column contains the Stokes number S ​ t 𝑆 𝑡 St , time step h ℎ h and simulation time T 𝑇 T used in the six simulations. The first row contains the aspect ratio λ 𝜆 \lambda of the particle shape.

		$\lambda =1$			$\lambda =5$
	$P_n$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$
$St=\frac{1}{10}$	MRBF1+CP2	0.1285	0.0860	0.3063	0.0514	0.0478	0.5389
	MRBF2+CP2	0.1002	0.0388	0.1256	0.0967	0.0657	0.3578
$h=\frac{1}{100}$	MRBF3+CP2	0.0380	0.0164	0.0614	0.0337	0.0242	0.1732
	TP1+CP2	7.9636	0.7427	0.7957	3.1871	0.5657	1.1261
$T=4$	TP2+CP2	3.7166	0.2975	0.3959	1.5779	0.3409	0.8653
	TP3+CP2	0.3349	0.0625	0.1197	0.0778	0.0704	0.3468
$St=1$	MRBF1+CP2	0.9861	0.5802	1.4746	0.0402	0.0965	1.3596
	MRBF2+CP2	0.0444	0.0439	0.5004	0.0351	0.0698	0.7334
$h=\frac{1}{40}$	MRBF3+CP2	0.0367	0.0353	0.2442	0.0258	0.0564	0.3755
	TP1+CP2	1.6501	0.7003	1.5487	0.4281	0.3222	1.8162
$T=8$	TP2+CP2	1.5784	1.2164	1.8284	0.0585	0.1871	1.7269
	TP3+CP2	0.0404	0.0440	0.3332	0.0310	0.0714	0.8375
$St=10$	MRBF1+CP2	1.5452	0.0979	0.1223	0.1401	0.0726	0.1998
	MRBF2+CP2	0.1071	0.0109	0.0154	0.0178	0.0183	0.0548
$h=\frac{1}{10}$	MRBF3+CP2	0.0416	0.0041	0.0065	0.0112	0.0084	0.0216
	TP1+CP2	5.8812	0.1545	0.1674	0.6419	0.1957	0.5617
$T=12$	TP2+CP2	4.0693	0.1540	0.2427	0.0955	0.0537	0.2222
	TP3+CP2	0.7464	0.0135	0.0125	0.0339	0.0242	0.0946

Table 7. Table 8 : The relative entropy E ​ ( P n ) 𝐸 subscript 𝑃 𝑛 E(P_{n}) , first Wasserstein distance W ​ ( P n ) 𝑊 subscript 𝑃 𝑛 W(P_{n}) and average error per particle Δ ​ 𝐱 n ¯ ¯ Δ subscript 𝐱 𝑛 \overline{\Delta\mathbf{x}_{n}} between the numerical distribution P n subscript 𝑃 𝑛 P_{n} and the reference distribution. The numerical distributions are calculated by various combinations of integration and interpolation methods as shown in the second column. The first column contains the Stokes number S ​ t 𝑆 𝑡 St , time step h ℎ h and simulation time T 𝑇 T used in the six simulations. The first row contains the aspect ratio λ 𝜆 \lambda of the particle shape.

		$\lambda =1$			$\lambda =10$
	$P_n$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$	$E\left(P_n\right)$	$W\left(P_n\right)$	$\overline{\Delta {𝐱}_n}$
$St=\frac{1}{10}$	MRBF1+CP1	0.1305	0.1050	0.3351	0.0423	0.0431	0.5453
	TP1+FE1	7.9359	0.6340	0.7866	3.2075	0.7961	1.2645
$h=\frac{1}{100}$	MRBF2+CP2	0.1249	0.0414	0.1242	0.0706	0.0456	0.3513
	TP2+CP2	2.7166	0.2998	0.3990	1.2931	0.3131	0.8785
$T=6$	TP2+AB2	2.6957	0.2836	0.3805	1.2585	0.2852	0.8833
	TP3+AB2	2.9179	0.1339	0.2512	0.1326	0.0804	0.4506
$St=1$	MRBF1+CP1	0.9463	0.5888	1.4858	0.0395	0.1000	1.3961
	TP1+FE1	4.3703	1.4721	2.1060	0.4507	0.5446	1.9263
$h=\frac{1}{40}$	MRBF2+CP2	0.0507	0.0525	0.5055	0.0306	0.0753	0.7845
	TP2+CP2	1.6212	1.2123	1.8251	0.0534	0.1630	1.7856
$T=8$	TP2+AB2	1.1982	1.0943	1.7759	0.0532	0.1595	1.7693
	TP3+AB2	0.0589	0.0532	0.4793	0.0376	0.0865	1.0213
$St=10$	MRBF1+CP1	1.2748	0.3379	0.5122	0.0727	0.1190	0.7322
	TP1+FE1	6.1882	0.7767	0.8750	1.9416	0.6882	1.5169
$h=\frac{1}{10}$	MRBF2+CP2	0.0979	0.0348	0.0663	0.0262	0.0438	0.2688
	TP2+CP2	3.2931	0.4742	0.7083	0.0754	0.1176	0.8004
$T=16$	TP2+AB2	3.2219	0.4063	0.6304	0.0774	0.1194	0.7933
	TP3+AB2	0.2631	0.0575	0.0865	0.0397	0.0594	0.4864

Table 8. Table 9 : The expressions for the constants χ 0 subscript 𝜒 0 \chi_{0} , α 0 subscript 𝛼 0 \alpha_{0} , β 0 subscript 𝛽 0 \beta_{0} and γ 0 subscript 𝛾 0 \gamma_{0} for λ < 1 𝜆 1 \lambda<1 , λ = 1 𝜆 1 \lambda=1 and λ > 1 𝜆 1 \lambda>1 .

		$\lambda =1$	$\lambda >1$
${\chi }_0$	$\frac{{\lambda }^2\left(\pi -{\kappa }_0\right)}{\sqrt{1-{\lambda }^2}}$	$2$	$\frac{-{\kappa }_0\lambda }{\sqrt{{\lambda }^2-1}}$
${\alpha }_0={\beta }_0$	$\frac{-\lambda \left({\kappa }_0-\pi +2\lambda \sqrt{1-{\lambda }^2}\right)}{2{\left(1-{\lambda }^2\right)}^{3/2}}$	$\frac{2}{3}$	$\frac{{\lambda }^2}{{\lambda }^2-1}+\frac{\lambda {\kappa }_0}{2{\left({\lambda }^2-1\right)}^{3/2}}$
${\gamma }_0$	$\frac{\left(\lambda \left({\kappa }_0-\pi \right)+2\sqrt{1-{\lambda }^2}\right)}{{\left(1-{\lambda }^2\right)}^{3/2}}$	$\frac{2}{3}$	$\frac{-2}{{\lambda }^2-1}-\frac{\lambda {\kappa }_0}{{\left({\lambda }^2-1\right)}^{3/2}}$
${\kappa }_0$	$2\arctan \left(\frac{\lambda }{\sqrt{1-{\lambda }^2}}\right)$	$1$	$\ln \left(\frac{\lambda -\sqrt{{\lambda }^2-1}}{\lambda +\sqrt{{\lambda }^2-1}}\right)$

Equations

$$\dot{\mathbf{v}}=$$

$$\dot{\mathbf{x}}=$$

$$J\dot{\mathbf{\omega }}=J\mathbf{\omega }\times \mathbf{\omega }-\mathbf{T},$$

$$\dot{Q}=Q\hat{\mathbf{\omega }},$$

$$\left(\begin{array}[]{c}\omega_{1}\\ \omega_{2}\\ \omega_{3}\\ \end{array}\right)\mapsto\widehat{\boldsymbol{\omega}}=\left(\begin{array}[]{ccc}0&-\omega_{1}&\omega_{2}\\ \omega_{1}&0&-\omega_{3}\\ -\omega_{2}&\omega_{3}&0\\ \end{array}\right),$$

$$\dot{\mathbf{x}}=\mathbf{u}(\mathbf{x})-{\alpha }^{-1}\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)+O({\alpha }^{-2})$$

$$\nabla \cdot \mathbf{v}=\frac{\partial u_i}{\partial x_i}-\frac{1}{\alpha }\left(\frac{\partial }{\partial t}\frac{\partial u_i}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}+u_i\frac{\partial }{\partial x_i}\frac{\partial u_j}{\partial x_j}\right)+O({\alpha }^{-2})$$

$$\nabla \cdot \mathbf{v}=-\frac{1}{\alpha }\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}=-\frac{1}{\alpha }\left(\Vert S{\Vert }_F^2-\Vert \Omega {\Vert }_F^2\right)+O({\alpha }^{-2})$$

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\text{and}{\Omega }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right)$$

$$\dot{\mathbf{y}}=\mathbf{f}(\mathbf{y})-A\mathbf{y}$$

$$Y=\det \left(\frac{\partial \mathbf{y}(t)}{\partial {\mathbf{y}}_0}\right)$$

$$Y=\prod _{i=1}^n{e}^{t{\lambda }_i}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \mathbf{y}(t)}{\partial {\mathbf{y}}_0}=({\mathbf{f}}^{\prime }-A)\frac{\partial \mathbf{y}(t)}{\partial {\mathbf{y}}_0}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\det (M(t))=\det (M(t))\mathrm{Tr}\left(M(t{)}^{-1}\frac{\mathrm{d}}{\mathrm{d}t}M(t)\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}Y=-Y\mathrm{Tr}\left(A\right),$$

$$Y=e^{-t\mathrm{Tr}(A)}.$$

$$\sum _{i=1}^n{\lambda }_i=-\mathrm{Tr}(A).$$

$$\dot{\mathbf{y}}=B\mathbf{y}-A\mathbf{y},$$

$${\Phi }_h^{[RK]}({\mathbf{y}}_0)=R(h(B-A)){\mathbf{y}}_0.$$

$${\Phi }_h^{[RK]}({\mathbf{y}}_0)=\exp \left(h(B-A)\right){\mathbf{y}}_0+O(h^{p+1})$$

$$\det \left(\frac{\partial {\Phi }_h^{[RK]}({\mathbf{y}}_0)}{\partial {\mathbf{y}}_0}\right)=e^{-h\mathrm{Tr}(A)}+O(h^{p+1}).$$

$$\det \left(\frac{\partial {\Phi }_h^{[FE]}({\mathbf{y}}_0)}{\partial {\mathbf{y}}_0}\right)=1-3\alpha h+\alpha h^2\left(3\alpha -\frac{\partial u_i}{\partial x_i}\right)+O(h^3),$$

$$\det \left(\frac{\partial \mathbf{y}(h)}{\partial {\mathbf{y}}_0}\right)=e^{-3\alpha h}.$$

$$h⪅\frac{3}{3\alpha -\frac{\partial u_i}{\partial x_i}}.$$

$$\Phi^{[RK]}_{h}(\mathbf{y}_{0})=\mathbf{y}_{0}+h\bigl{(}(1-\tfrac{1}{2\theta})\mathbf{f}(\mathbf{y}_{0})+\tfrac{1}{2\theta}\mathbf{f}(\mathbf{y}_{0}+\theta h\mathbf{f}(\mathbf{y}_{0})),$$

$$\det \left(\frac{\partial {\Phi }_h^{[RK]}({\mathbf{y}}_0)}{\partial {\mathbf{y}}_0}\right)$$

$$+\left(3\alpha (\theta -1)\frac{{\partial }^2{u}_i}{\partial x_i\partial x_j}{v}_j+3{\alpha }^2\frac{\partial u_i}{\partial x_i}-24{\alpha }^3\right)\frac{h^3}{3!}+O(h^4)$$

$${\mathbf{y}}_k$$

$${\mathbf{y}}_i$$

$${\mathbf{y}}_2=$$

$$\det \left(\frac{\partial {\mathbf{y}}_2}{\partial {\mathbf{y}}_0}\right)$$

$$+\left(\alpha \left(3\theta -\frac{21}{2}\right)\frac{{\partial }^2{u}_i}{\partial x_i\partial x_j}{v}_j+\frac{21{\alpha }^2}{2}\frac{\partial u_i}{\partial x_i}-24{\alpha }^3\right)\frac{h^3}{3!}+O(h^4)$$

Full text

1 Introduction

Since the influential work of Maxey and Riley [1] in deriving the equations of motion of an inertial spherical particle immersed in viscous flow, there have been a multitude of studies exploring the collective behavior of suspensions of particles. In particular, the remarkable phenomenon of preferential concentration of inertial particles in turbulence has attracted the attention of many authors. This phenomenon, sometimes referred to as the “centrifuge effect”, is also attributed to Maxey [2] who showed that particles disperse in regions where the fluid velocity strain rate is low compared to the vorticity. The theoretical mechanisms for particle clustering has since been further explored by means of Lyapunov exponent analysis [3], caustics [4] and perturbative methods to name a few. Sophisticated numerical simulations [5] have also advanced and verified our understanding of this phenomenon for a variety of flows and extended such observations to non-spherical particles [6]. As the need for large-scale simulations increase, the demand for cost-effective numerical methods is growing. However, despite the fact that numerical simulations are so well documented, there have been few studies that explore the extent to which the numerical methods used in simulations accurately reproduce the geometric properties that explain the preferential clustering of particles. In this paper we discuss some features of the equations of motion that influence the preferential concentration of particles and determine to what extent these features can be replicated by well designed numerical methods. In doing so, we propose an efficient numerical algorithm that is designed to replicate these features. The method combines matrix-valued radial basis functions for the divergence-free interpolation of the discrete fluid field with a splitting method that is designed specifically for the equations of motion under study.

Interpolation methods are necessary for simulating suspensions of particles as the flow field is usually generated by a direct numerical simulation of the Navier-Stokes equations and is therefore only available at discrete points in space, meaning that it must be approximated at the location of the particle. To achieve this in a simple and efficient manner many authors use a variant of a tri-polynomial interpolant, for example [5, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Previous studies [16, 17] have explored the extent to which these interpolation methods accurately reproduce statistical properties of the turbulent flow field. However, all the interpolation methods considered in the aforementioned references are based on polynomials that create an approximation to the fluid velocity field that is not divergence-free. One major consequence is that the hydrodynamic Stokes force that determines the particle path lines is instead calculated from a non-conservative fluid velocity field. These non-conservative interpolation methods are still used in practice today despite the fact that the theoretical mechanisms that explain the preferential concentration are derived with the assumption of incompresibility. Furthermore, it is often argued (e.g., [11]) that interpolation errors are “averaged out” and it is concluded that one can acheive statistically similar results using a fast low-order interpolation method. This claim is supported by the fact that linear interpolation produces similar statistics to simulations using cubic interpolation [16]. However, neither linear nor cubic interpolation preserves the divergence-free condition of the fluid field and therefore it is not truly understood whether or not errors to the divergence of the fluid field are averaged out in the same way that standard truncation errors are. The implications of these divergence-errors have not been studied in detail, however there is numerical evidence suggesting that breaking this condition can lead to erroneous clustering in PDF methods, first presented in [18] and also in [19]. Divergence free interpolation has been used to good effect in particle-laden flow simulations [20, 21] as well as in other particle simulation problems, such as in geodynamic modelling [22] and magnetospheric physics [23], for example. One of the goals of this study is to explain, from a numerical analysis point of view, the consequences of breaking the divergence-free condition in the flow field. We also show the benefit of divergence-free interpolation using in simulations of suspensions of inertial particles as well as show how these errors affect the numerical time integration.

In addition, we study some numerical integration methods and how their errors affect the preferential concentration of particles. Two popular classes of methods are explicit Runge-Kutta [7, 9] and Adams-Bashforth methods [8, 10, 12, 14, 13]. As the accuracy of the time stepping algorithm is limited by the interpolation error we consider only explicit order one and two methods. Such methods often do a reasonably good job at integrating the ODEs under study as they are efficient and easy to implement. However, the exact solution to the ODEs that govern the dynamics of particles with Stokes drag possess a number of physical features that can be exploited to increase the accuracy of the time stepping methods without increasing its order or cost. Such features include constant contractivity of phase space volume, the centrifuge effect, rigid body motion, linear dissipation and, in some cases, perturbative forces. These features are able to be exploited by a carefully designed splitting method. In this work we propose, as an alternative to Runge-Kutta and Adams-Bashforth methods, a splitting method that is especially designed to reduce the error in the centrifuge effect, which combined with divergence free interpolation techniques allow us to obtain a higher lever of accuracy in the distribution of particles in viscous flows.

1.1 Main contributions and summary of paper

We now highlight the main contributions and give a brief outline of the paper. We begin by outlining the equations of motion and the centrifuge effect in section 2. In section 3 we develop and analyze a contractive splitting method whose flow preserves the sum of the Lyapunov spectrum of the exact solution and show that conventional methods cannot do this. The splitting method is then applied to the equations of motion for spherical particles and the so-called “centrifuge-preserving” methods are presented, which are constructed to minimize the error of the centrifuge effect.

Section 4 presents the use and implementation of matrix-valued radial basis function interpolation to construct a divergence-free interpolation of the discrete flow field. We show that a vector field approximated by matrix-valued radial basis functions are compatible with the Stokes equations due to the fact they they are identical to the method of regularized Stokeslets. This results in a more physically realistic approximation to the underlying Navier-Stokes equations.

In section 5 we focus our attention to how physical volume of the particle phase $\Psi$ evolves over a small time $h$. Upon expanding $\Psi$ in $h$ under the exact solution, we recover the centrifuge effect at $O(h^{4})$. When expanding $\Psi$ under the numerical solution, we find that errors to the divergence of the fluid velocity field appear at $O(h^{2})$, overshadowing the centrifuge effect. However, when a divergence-free interpolation method is used, all the numerical methods under consideration replicate the qualitative behavior of the centrifuge effect. That is, physical volumes of particles will contract in regions where the vorticity is lower than the strain rate and vice versa, however, they do so at a slightly erroneous rate. To account for this error, we show that the centrifuge-preserving methods contract physical volume at the same rate as the exact solution to leading order in $h$, hence also preserving the quantitative behavior of the centrifuge-effect.

Section 6 is dedicated to simulations of particle suspensions evolving in a discrete cellular flow field where we compare the proposed geometric methods against conventional methods. What we observe is that a computationally inexpensive combination of divergence-free interpolation and centrifuge-preserving splitting methods yield far more accurate spatial distributions of particles compared to standard methods of higher cost. We present many examples where our geometric algorithm produced distributions of particles that are more similar to the “exact” distribution despite having higher error per particle than distributions produced by slow conventional methods. The main conclusion here is that numerical solutions that preserve the sum of the Lyapunov spectrum, the contractivity of phase space volume, the divergence-free condition and the centrifuge effect in simulations is of great benefit.

Section 7 is dedicated to conclusions.

2 The equations of motion

The translational dynamics of a small particle immersed in a viscous fluid is governed by the rigid body equations with a Stokes force term

[TABLE]

where $\mathbf{u}(\mathbf{x})$ is the fluid velocity at the particle’s location $\mathbf{x}$, $\mathbf{v}$ the velocity, $K$ is a positive definite resistance tensor and $\alpha=1/St$ is the inverse particle Stokes number, which is a dimensionless measure of particle inertia. Note that unless mentioned we will assume that $\mathbf{u}(\mathbf{x})$ does not explicitly depend on $t$. Doing so improves the readability and presentation of the paper and does not affect the forthcoming results.

For spherical particles, $K=I$ is the identity, the rotational variables are constant and the above ODEs uniquely specify the dynamics of each particle. For non-spherical particles, the resistance tensor $K=Q^{T}K_{b}Q$, where $K_{b}$ is the diagonal positive definite body frame resistance tensor and $Q\in SO(3)$ is a rotation matrix that transforms a vector in the body frame to one in the inertial frame. The angular velocity $\boldsymbol{\omega}$ evolves via

[TABLE]

where $J$ is the diagonal body frame moment of inertia tensor and $\mathbf{T}$ is the hydrodynamic torque. The rotation matrix $Q$ is calculated by solving the matrix ODE

[TABLE]

where $\widehat{\cdot}:\mathbb{R}^{3}\rightarrow\mathfrak{so}(3)$ is defined by

[TABLE]

such that $\widehat{\boldsymbol{\omega}}\mathbf{v}=\boldsymbol{\omega}\times\mathbf{v}$. The expressions for $K_{b}$ and $\mathbf{T}$ for spheroidal particles are given in A.

2.1 The centrifuge effect

Here we will outline the centrifuge effect of the particle equations of motion, which is one of the mechanisms for particle clustering that is referred to throughout the paper. In [2], Maxey assumes $\alpha\gg 1$ and expands the the spherical particle ODEs (1) and (2) in powers of $\alpha^{-1}$ to derive a first-order ODE expression for $\mathbf{x}$

[TABLE]

where we have ignored the effect of gravity. Taking the divergence gives

[TABLE]

where there is an implied summation over repeated indices, which is the convention that is assumed throughout the paper. Assuming that the fluid field is divergence-free, we arrive at the familiar relationship between the fluid field rate of strain, rate of rotation and the divergence of the particle velocity field

[TABLE]

where the rate of strain and rotation tensors $S$ and $\Omega$ are given by

[TABLE]

and $\|\cdot\|_{F}$ is the Frobenius matrix norm. In other words, the divergence of the particle velocity field $\nabla\cdot\mathbf{v}$ is positive when the vorticity is large compared to the strain rate tensor meaning that the particulate phase disperses in these regions. Conversely, particles concentrate in regions where the strain rate is large compared to the vorticity. This phenomenon is the “centrifuge effect” of the exact solution to (1) and (2).

Finally, we remark that while the centrifuge effect was derived for spherical particles, one can make similar observations for non-spherical particles. In this scenario, the resistance tensor can be decomposed into a spherical part and a non-spherical part, e.g., for a spheroidal particle with rotational symmetry (see A) we can write $K_{b}=a\,I+b\,\mathbf{e}_{z}\mathbf{e}_{z}^{T}$, where $\mathbf{e}_{z}=(0,0,1)^{T}$ and $b\rightarrow a$ in the spherical limit. In other words, the centrifuge effect still plays a central role in the preferential clustering of non-spherical particles in addition to the non-spherical effects due to $b\mathbf{e}_{z}\mathbf{e}_{z}^{T}$ term in the resistance tensor.

3 Numerical integration of dissipative vector fields

The dynamics of small inertial particles (both spherical and non-spherical) can be modeled as the flow of a vector field with linear dissipation. Such vector fields arise due to the fact that for low Reynolds number flow the drag forces are linear in the slip velocity, for example the Stokes drag force for small ellipsoids, spheres or rigid slender particles [24]. We begin this section with a discussion of such linearly dissipative vector fields and their contractive properties of phase space volume. We then discuss the application of some conventional explicit methods for integrating such ODEs. In particular, we show that conventional methods cannot preserve the contractivity of phase space volume. A splitting scheme is then shown to preserve the exact contractivity of phase space volume. The section concludes with the application of the splitting scheme to spherical particles.

3.1 Linearly dissipative vector fields and contractivity of phase space volume

A linearly dissipative vector field in $n$ dimensions is given in general by

[TABLE]

where $A$ is a positive definite matrix and $\mathbf{f}(\mathbf{y})$ is volume preserving, that is, it satisfies $\nabla\cdot\mathbf{f}(\mathbf{y})=0$ (e.g., any Hamiltonian vector field). Note that the ODEs of both non-spherical and spherical particles can be cast in this form, where $\mathbf{f}(\mathbf{y})$ represents the free rigid-body vector field plus the conservative part of the Stokes force and $-A\mathbf{y}$ represents the dissipative part of the Stokes force.

The quantitative behavior of particle clustering can be explained in part by analyzing the Lyapunov exponents $\lambda_{i}$ of the ODE, see for example [3]. It is therefore desirable that the numerical solution of the ODE reproduces similar Lyapunov exponent characteristics. Whilst there do not currently exist numerical methods that preserve individual Lyapunov exponents we can however construct a numerical method that preserves the sum of the Lyapunov spectrum $\sum_{i=1}^{n}\lambda_{i}$. From a backward error analysis point of view, a numerical method that preserves the Lyapunov spectrum is one that is the exact solution to an ODE with the same Lyapunov spectrum sum as the ODE being solved. For the equations of motion for spherical particles (equations (1) and (2)), the sum of the first three Lyapunov exponents characterizes the divergence of the velocity field and the sum of the spatial Lyapunov exponents characterizes the rate at which particle clouds contract or expand [20]. Generally speaking, the sum of the Lyapunov spectrum describes the rate at which phase space volume exponentially contracts or expands [25]. That is, by letting $\mathbf{y}(t)$ denote the exact solution of (10) with initial conditions $\mathbf{y}(0)=\mathbf{y}_{0}$ and

[TABLE]

the $n$ dimensional phase space volume, then

[TABLE]

It is also known that linearly dissipative systems contract phase space volume at a constant rate, as we will now show. Taking the Jacobian of $\mathbf{y}(t)$ with respect to $\mathbf{y}_{0}$ gives

[TABLE]

We now recall Jacobi’s formula, which relates the derivative of the determinant of a square matrix $M(t)$ by the following

[TABLE]

Differentiating (11) with respect to time and applying (14) gives

[TABLE]

as $\mathbf{f}^{\prime}$ has zero trace. This is solved by

[TABLE]

By equating this with (12) we obtain the relation

[TABLE]

As the trace of $A$ is by definition positive, the phase space volume $Y$ is strictly monotonically contracting in time. Equation (17) implies that a numerical integration method that preserves phase space volume $Y$ also preserves the sum of the Lyapunov exponents of the underlying ODE. It is therefore logical to imply that a numerical flow of (10) that preserves the contractivity of phase space volume will better reproduce the clustering properties of the exact solution than one that doesn’t. The rest of this section is dedicated to analysing to what extent some common numerical integration methods for particle dynamics can preserve this constant contractivity of phase space volume.

3.2 Preservation of the contractivity of phase space volume by numerical methods

Denote by $\Phi_{h}$ a numerical method for solving (10) such that ${\Phi_{h}(\mathbf{y}_{0})\approx\mathbf{y}(h)}$ for time step $h\ll 1$. For $\Phi_{h}$ to be called contractivity preserving when applied to (10), we require that $\det\left(\frac{\partial\Phi_{h}(y)}{\partial y}\right)=e^{-h\,\mathrm{Tr}(A)}$ [26]. It is known that no standard methods (e.g., one with a B-series [27]) can preserve phase space volume for divergence free vector fields [28]. The same is also expected when it comes to preserving contractivity of phase space volume for dissipative vector fields [29]. Instead, a weaker requirement is that they contract phase space volume when the exact solution does so, that is, $\det\left(\frac{\partial\Phi_{h}(y)}{\partial y}\right)<1$ when applied to (10). Such a numerical method that possesses this property is called contractive.

We now consider some popular numerical methods for particle dynamics. We consider explicit methods used in the literature, namely, explicit Runge-Kutta methods, Adams-Bashforth methods and a splitting method scheme.

3.2.1 Runge-Kutta methods and phase volume contractivity

Take an order-$p$ Runge-Kutta method $\Phi^{[RK]}_{h}(\mathbf{y}_{0})$ with stability function $R(z)$ applied to an ODE of the form (10) with linear $\mathbf{f}(\mathbf{y})$, say

[TABLE]

where $B$ is a square and traceless matrix. Then the numerical solution by the Runge-Kutta method is given by

[TABLE]

Recall that $R(z)$ is an order-$p$ Padé approximation to the exponential function. We therefore have

[TABLE]

which means that over one time-step, the phase space volume contracts via

[TABLE]

That is, for such a linear system, a Runge-Kutta method will only preserve the phase volume contractivity up to the order of the method. So for a non-linear dissipative ODE of the form (10), one can hardly expect a Runge-Kutta method to preserve phase space volume exactly. In fact, due to this error explicit Runge-Kutta methods usually have a time-step restriction on $h$ to even be contractive at all [26]. This is illustrated by the following examples of some low order explicit Runge-Kutta methods applied to equations (1) and (2).

Example 1**.**

We apply the forward Euler method $\Phi^{[F\!E]}_{h}$ to the ODEs (1) and (2). Note that these ODEs are non-linear due to $\mathbf{u}(\mathbf{x})$. Setting $\mathbf{y}:=(\mathbf{v},\mathbf{x})$, we have for the contractivity of phase space volume under the forward Euler method

[TABLE]

which is an order one approximation to the exact contractivity

[TABLE]

The Forward Euler method must therefore satisfy the following time-step restriction for it to be contractive

[TABLE]

Violating this restriction means that the forward Euler method will expand phase space volume despite the ODE dictating that it is always contracting. Furthermore, it can be seen that large values of $|\frac{\partial u_{i}}{\partial x_{i}}|$ will place further restrictions on the size of $h$.

Example 2**.**

Consider the following second order explicit Runge-Kutta method

[TABLE]

where $\theta=\frac{1}{2},\frac{2}{3}$ and $1$ correspond to the explicit midpoint method, Ralston’s method and Heun’s method, respectively. We apply this to the ODEs (1) and (2). Setting $\mathbf{y}:=(\mathbf{v},\mathbf{x})$, we have for the contractivity of phase space volume under $\Phi^{[RK]}_{h}$

[TABLE]

which is an order two approximation to the exact contractivity (23), which is expected for an order two method. The time-step $h$ must be chosen small enough such that the $O(h^{3})$ error term does not violate the contractivity condition. Violating this restriction means that the method will expand phase space volume despite the ODE dictating that it is always contracting. Furthermore, it can be seen that large values of $|\frac{\partial u_{i}}{\partial x_{i}}|$ and $v_{i}$ will place further restrictions on the size of $h$.

We remark that we can make similar observations for the above methods applied to the ODEs for non-spherical particles. That is, the phase space volume is conserved only to the order of the method.

3.2.2 Multi-step methods and phase volume contractivity

Another popular numerical method used for particle dynamics are multi-step methods. Consider the explicit $k$-step Adams-Bashforth methods. Such methods are of global order-$k$ and can be seen as a map $\Phi_{h}^{[AB]}:(\mathbf{y}_{0},...,\mathbf{y}_{k-1})\rightarrow(\mathbf{y}_{1},...,\mathbf{y}_{k})$ such that

[TABLE]

where the coefficients $b_{i}$ satisfy $\sum_{i=0}^{k-1}b_{i}=1$. That is, $\Phi_{h}^{[AB]}$ takes a point in a $kn$ dimensional phase space to another point in the same space. Due to this and the fact that the initial vectors in the domain $\mathbf{y}_{i}$ for $i=0,...,k-1$ are independent of one another, it’s less clear how to define the notion of numerical phase space volume that relates to that of the underlying ODE. However, in practice these initial vectors $\mathbf{y}_{k}$ for $i=0,...,k-1$, are usually computed by an order-$k$ one step method, for example a Runge-Kutta method. Therefore, the vectors $\mathbf{y}_{i}$ depend on the vectors $\mathbf{y}_{j}$ for $j<i$. This implies that each $\mathbf{y}_{i}$ has the same series expansion as the exact solution up to $O(h^{k})$. While a detailed analysis of multi-step methods and the preservation of phase volume lie outside the scope of the paper, we can illustrate this concept by the following example, which considers the phase volume properties of a second order Adams-Bashforth using a second order Runge-Kutta method to compute the initial vectors.

Example 3**.**

Consider the second-order Adams-Bashforth method $\Phi_{h}^{[AB]}:(\mathbf{y}_{1},\mathbf{y}_{0})\rightarrow(\mathbf{y}_{2},\mathbf{y}_{1})$ where

[TABLE]

Now define $\mathbf{y}_{1}=\Phi^{[RK]}_{h}(\mathbf{y}_{0})$ in the domain by the second order Runge-Kutta method (25). Applying $\Phi_{h}^{[AB]}$ to the ODEs (1) and (2) and taking the Jacobian determinant of $\mathbf{y}_{2}$ with respect to $\mathbf{y}_{0}$ gives

[TABLE]

which is an $O(h^{2})$ approximation to the exact contractivity (23). Note that here we have taken the the contractivity over two time-steps $2h$. Like the previous example, we observe that the contractivity is affected by non-zero values of $\frac{\partial u_{i}}{\partial x_{i}}$ and $v_{i}$.

3.3 A splitting scheme that preserves the contractivity of phase space volume

We now analyze the splitting method based on the following splitting of equation (10)

[TABLE]

where $\mathbf{b}(\mathbf{y})$ is any vector that is constant along the flow of the second vector field. A similar splitting was proposed in [30] for non-spherical particle dynamics. Denote their exact flow operators by $\psi^{[1]}_{h}$ and $\psi^{[2]}_{h}$, respectively. In the context of small particles immersed in a viscous fluid, the first vector field represents the free rigid body equations and the second is due to the Stokes viscous drag forces. The free rigid body vector field can be solved exactly. That is, by a forward Euler step for the spherical case or otherwise using trigonometric or Jacobi elliptic functions depending whether or not the body is axially symmetric [31]. Due to the existence of an exact solution, we immediately have volume preservation

[TABLE]

The second vector field is solved by the variation of parameters formula

[TABLE]

Taking the Jacobian determinant gives

[TABLE]

which is consistent with the exact solution (16). As the Jacobian of the composition of two or more maps is the product of the Jacobians of the maps, any splitting method based on the alternating compositions of the flows $\psi^{[1]}_{h}$ and $\psi^{[2]}_{h}$ will be contractivity preserving.

In forthcoming numerical experiments, we will consider only order one and two methods including the order one Lie-Trotter method

[TABLE]

and the order two Strang method

[TABLE]

Here, we denote by $\Phi^{[LT]*}_{h}=\psi^{[2]}_{h}\circ\psi^{[1]}_{h}$ the conjugate of the $\Phi^{[LT]}_{h}$.

3.4 Application to spherical particle dynamics and the centrifuge-preserving methods

To construct a contractivity preserving splitting method for spherical particles, we split the ODEs (1) and (2) in to the following two vector fields

[TABLE]

Their exact flow operators are

[TABLE]

Indeed, letting $\mathbf{y}_{0}=(\mathbf{v}_{0}^{T},\mathbf{x}_{0}^{T})^{T}$ we see that

[TABLE]

hence any composition of the above flows will preserve contractivity. For the construction of the splitting method for non-spherical particle dynamics, we refer the reader to [30].

A draw back of splitting methods is that composing methods of order higher than two requires the use of negative time steps. As the ODEs in question are dissipative, such higher order methods would therefore require strict time step restrictions, which is preferably avoided. The idea behind geometric numerical integrators is that preserving relevant properties of the exact solution upon discretisation can lead to better qualitative and long time numerical solutions. With this in mind, instead of improving the accuracy of the method in a conventional sense by increasing the order, we propose as an alternative the following composition methods

[TABLE]

which are order one and order two methods, respectively. We propose that the above splitting methods are particularly well suited to calculation of particle dynamics as their numerical solution preserves the centrifuge effect of the exact solution when considering the contraction of physical volume of the particle field. We will therefore refer to the methods (42) and (43) as the “centrifuge-preserving” methods. This favorable property is discussed in more detail in section 5.2. In section 6.4 we show through numerical simulations that integrators possessing this property predict more accurately the spatial distribution of particles (both spherical and non-spherical) compared to methods without this property.

4 Divergence-free interpolation with matrix-valued radial basis functions

To construct a divergence-free approximation to the discrete fluid field, we propose using matrix-valued radial basis functions (MRBFs). In this section, we will give a brief outline on their use and implementation. We then further motivate their use by showing that the interpolated vector field generated by MRBFs is a solution to the Stokes equation.

The interpolation problem is as follows. Given a set of vector-valued data $\{\mathbf{u}_{i},\mathbf{x}_{i}\}_{i=1}^{n^{3}}$ generated by an accurate direct numerical simulation to the Navier-Stokes equations, construct a divergence-free vector field that locally interpolates the data. In our context, $\mathbf{u}_{i}=\mathbf{u}(\mathbf{x}_{i})$ is the fluid velocity vector at the grid node located at $\mathbf{x}_{i}=(x_{i},y_{i},z_{i})^{\mathrm{T}}$. When implementing a polynomial interpolation method, one usually chooses the $n\times n\times n$ cube of data points neighboring the particle, where $n=2,3$ or $4$. This is because polynomial interpolation of degree $n-1$ requires $n$ data points in each dimension to specify a unique interpolating polynomial. MRBFs are not restricted by this particular choice of data points, however to keep the interpolation methods comparable we will adopt this convention. The MRBF interpolating vector field $\mathbf{s}(\mathbf{x})$ is then constructed by

[TABLE]

where

[TABLE]

is called an MRBF, $\theta(r_{i}(\mathbf{x}))$ is a (scalar-valued) radial basis function, $r_{i}(\mathbf{x})=\|\mathbf{x}_{i}-\mathbf{x}\|$ is the distance from the point $\mathbf{x}_{i}$ and $I$ is the identity matrix in three dimensions. The $n^{3}$ vector-valued coefficients $\mathbf{c}_{i}\in\mathbb{R}^{3}$ are chosen such that $\mathbf{s}(\mathbf{x}_{i})=\mathbf{u}(\mathbf{x}_{i})$, which amounts to solving the following $3n^{3}$ dimensional linear system

[TABLE]

The particular RBF we use in the forthcoming experiments is the Gaussian $\theta(r)=\exp(-\epsilon^{2}r^{2})$, where $\epsilon$ is some user defined parameter that controls the flatness of the RBF. In general, one should choose $\epsilon$ as small as possible as this leads to less interpolation error, although more ill conditioned systems.

It can be easily seen that $\mathbf{s}(\mathbf{x})$ is divergence free. Using the double curl identity in $\mathbb{R}^{3}$ we have

[TABLE]

Finally, we list some advantages of MRBF interpolation over standard tri-polynomial interpolation: (1) they work equally well on scattered data points, meaning that they are just as well suited to interpolate data generated by a direct numerical simulation involving complex geometries on unstructured grids; (2) they have faster convergence of their derivatives [32, 33], compared to tri-polynomial interpolation [34]; and (3), they are compatible with the Stokes equations, meaning that they construct a more physically realistic fluid field for fluid simulations. We will discuss point (3) in the next section.

4.1 MRBFs as regularised Stokeslet solutions to the Stokes equations

In addition to the fact that the underlying flow field should be divergence-free, we are given extra knowledge that can be exploited; namely that the data is a numerical solution to the incompressible Navier-Stokes equations

[TABLE]

We are only interpolating in space and hence approximating steady-state solutions to (48) and as the grid-spacing $\Delta x$ is comparable to the smallest length scales of the flow (e.g., the Kolmogorov scale for turbulent flows), the Reynolds number is small and the non-linear terms of equations (48) can be ignored. Under the above assumptions, a good approximation for the local flow in a grid-cell can be given by the steady Stokes equations, which reads

[TABLE]

where we have set $\mu=1$. Cortez [35] presents what’s called the regularised Stokeslet solution to the Stokes equation, which is an approximation to Green’s function of the Stokes equation for body force $\mathbf{F}=\phi_{\epsilon}(\mathbf{x})\,\mathbf{f}_{0}$, where $\mathbf{f}_{0}\in\mathbb{R}^{3}$ is constant. Here, $\phi_{\epsilon}(\mathbf{x})$ is the so-called “blob” function, which is a radially symmetric smooth approximation of the Dirac delta function $\delta(\mathbf{x})$ that decays to zero at infinity whilst satisfying

[TABLE]

Now define the functions $G_{\epsilon}(\mathbf{x})$ and $B_{\epsilon}(\mathbf{x})$ as the solutions to

[TABLE]

which are smooth approximations to Green’s function and the biharmonic equation $\nabla^{4}B(\mathbf{x})=\delta(\mathbf{x})$, respectively. Then Cortez’s regularised Stokeslet solution reads

[TABLE]

with pressure term

[TABLE]

Using the definition for $G_{\epsilon}(\mathbf{x})$, we can rewrite the regularised Stokeslet (52) as

[TABLE]

which is identical to an MRBF element if we can identify $B_{\epsilon}(\mathbf{x})$ with a positive-definite RBF $\psi(||\mathbf{x}||)$ (e.g., the Gaussian $\psi(||\mathbf{x}||)=\exp\left(-\epsilon^{2}||\mathbf{x}||^{2}\right)$) and the force vectors are identified with the interpolation coefficient vectors $\mathbf{c}_{i}$ from equation (44). This means that a vector field that is constructed from a linear combination of MRBFs, (i.e., equation (44)) corresponds to a linear combination of regularised Stokeslet solutions, with force $\mathbf{f}_{i}=\mu\mathbf{c}_{i}$. This leads to the following solution to the Stokes equation, now written in terms of MRBFs

[TABLE]

One implication of this is that the interpolated background fluid field is related to the gradient of a scalar pressure field, when MRBF interpolation is used. The benefit of this can be illustrated by inserting equation (48) into equation (6) to derive an expression for $\nabla\cdot\mathbf{v}(\mathbf{x})$ in terms of the pressure field [36]

[TABLE]

This equation tells us that the pressure field is also related to the preferential concentration of particles. Moreover, this suggests that particles cluster in regions of maximum pressure ($\nabla^{2}p_{\epsilon}(\mathbf{x})<0$) [36]. Indeed, in [37], numerical evidence is found to support the correlation between the Laplacian of the pressure field $\nabla^{2}p_{\epsilon}(\mathbf{x})$ and the spatial distribution of the particles. However, if the background fluid field is interpolated by a standard polynomial method, then there is no background scalar pressure field, which could erroneously influence the particle path lines.

5 Numerical errors and preferential concentration of spherical particles

In what follows we will consider how volumes of inertial spherical particles evolve under the flow of the ODEs (1) and (2). The goal of this section is to relate the numerical interpolation and integration errors to the clustering mechanisms of the exact solution. This is done by first expanding the exact solution into its elementary differentials. We then discuss the effects of integration errors and interpolation errors on the evolution of volumes of particles. In what follows, we will initially assume that $\mathbf{u}(\mathbf{x})$ is an arbitrary vector field that is not necessarily divergence-free until explicitly mentioned.

5.1 Expanding the exact solution

We start this section by defining the notion of volume of the particle suspension. Given an open and bounded set $D_{t}\subset\mathbb{R}^{3}$ at time $t$, then its volume at $t=0$ is given by

[TABLE]

where $\mathbf{x}(0)=\mathbf{x}_{0}$. This can be thought of as the volume occupied by a suspension of inertial particles confined to the region $D_{0}$. The idea is to consider how this volume expands or contracts in time. Consider now the same volume after evolving under the ODEs (1) and (2) for time $t$

[TABLE]

Hence, the quantity

[TABLE]

determines how volumes of particles contract or expand over time. That is, given a volume of particles, if $\Psi>1$ the volume is expanding, $\Psi<1$ the volume is contracting and $\Psi=1$ the volume is preserved. These three cases correspond to the particulate phase dispersing, concentrating or remaining a constant density, respectively. Note that we will refer to $\Psi$ as the physical volume, to distinguish between phase space volume.

To illustrate the connection between $\nabla\cdot\mathbf{v}$ and $\Psi$, we can take the Jacobian of equation (2) with respect $\mathbf{x}_{0}$ [38]

[TABLE]

Applying Jacobi’s formula (14) yields a differential equation for $\Psi$

[TABLE]

It is clear that if $\nabla\cdot\mathbf{v}<1$, then $\Psi$ is decreasing and if $\nabla\cdot\mathbf{v}>1$ then $\Psi$ is increasing.

We will now show this more concretely, by expanding the exact solution into its elementary differentials, which we now recall. Denote by $\mathbf{y}(t)$ the exact solution of an ODE

[TABLE]

For some small time $0<h\ll 1$, $y_{i}(h)$ has the following elementary differential expansion [27]

[TABLE]

for $i=1,...,n$. We note that the expansion is convergent if $h$ is small compared to $\|\mathbf{f}\|$.

The elementary differentials of the ODEs (1) and (2) are calculated and the terms up to $O(h^{3})$ are presented

[TABLE]

where the variable appearing on the right hand side are evaluated at $t=0$. Here we assumed nothing about the size of $\alpha$, but instead take $h\ll\alpha$ such that the series converges. Taking the determinant of the Jacobian of $\mathbf{x}(h)$ from equation (68) with respect to $\mathbf{x}_{0}$ yields an expansion for $\Psi$ with repect to $h$

[TABLE]

where

[TABLE]

The elementary differentials $\chi_{i}$ are given by

[TABLE]

This can be verified by expansion of (61) into its Taylor series. If we insist that the fluid field is divergence-free then the $\Psi_{i}$ and $\chi_{i}$ all vanish except for $\Psi_{4}$ and $\chi_{4}$. We are then left with

[TABLE]

which relates the fluid rate of strain and rotation with the contractivity of physical volume in the same way as the centrifuge effect (8). That is, if the rate of vorticity is greater than the rate of strain, physical volumes of particles will contract and vice versa.

5.2 Expanding the numerical solution and numerical errors

In this section, we perform a similar analysis to that of section 5.1 but instead of the exact flow of the ODE, we consider now how physical volumes of particles evolve under the numerical flow. That is, we will look at how errors to the divergence of the fluid field affect the evolution of volumes of particles under the numerical solution to the equations of motion. We do so by expanding the numerical methods into their elementary differentials and comparing the expansions with the exact solution. We then look at how errors to the divergence of the fluid field affect the evolution of physical volume under the numerical flow. The results are that if a divergence-free interpolation method is used, the numerical methods preserve the same qualitative behavior of the centrifuge effect. Moreover, we show here that the centrifuge-preserving methods replicate the centrifuge-effect from equation (74) up to the accuracy of the interpolation method when the fluid field is divergence-free.

Consider the map $\Phi^{[n]}_{h}:(\mathbf{v}_{0},\mathbf{x}_{0})\rightarrow(\mathbf{v}_{1},\mathbf{x}_{1})$, where the superscript $[n]$ denotes the numerical method in consideration. We calculate $\Psi^{[n]}=\det\left(\frac{\partial\mathbf{x}_{1}}{\partial\mathbf{x}_{0}}\right)$ and expand the solution in $h$ yielding an expression of the form

[TABLE]

The values of $\Psi^{[n]}_{i}$ for $i=2,3,4$ for the Forward Euler (FE1), Lie-Trotter (LT1), order one centrifuge-preserving (CP1), Ralston (RK2), Adams-Bashforth two-step (AB2), order two centrifuge-preserving (CP2) and Strang splitting (SS2) methods are presented in table 1. Note that $\Psi^{[n]}_{1}=0$ for all the above methods.

We make a number of observations from this table. First, the divergence of the fluid field affects $\Psi^{[n]}$ at $O(h^{2})$ for each method. The one exception to this is FE1, which satisfies $\Psi^{[F\!E]}=1$ and therefore erroneously preserves physical volume. When the divergence of the fluid field is zero, all the $\chi_{i}=0$ except $\chi_{4}$. For example, setting $\nabla\cdot\mathbf{u}=0$ gives for the SS2 method

[TABLE]

This means that the numerical solution generated by the Strang splitting method (38) reproduces the qualitative nature of the centrifuge effect, in the sense that \Psi^{[SS]}\big{|}_{\nabla\cdot\mathbf{u}=0}>1 when $\left.\|\Omega\|_{F}>\|S\|_{F}\right.$. This qualitative centrifuge effect is seen by all the methods (other than FE1) by setting $\nabla\cdot\mathbf{u}=0$ in table 1. However, we note here that the coefficient of the $O(h^{4})$ term in equation (76) is different to that of the exact solution (74). Meaning that, while the method contracts physical volume when the exact solution does, it does so at an erroneous rate. This issue is circumvented by the centrifuge-preserving methods (CP1 and CP2), where we have chosen the time-step coefficients in such a way such that they yield the exact same expansion as (74) up to $O(h^{4})$ and hence contracts physical volume at the same rate as the exact solution to leading order.

We now discuss the effect of interpolation errors in simulations of spherical particles in numerically calculated flows. Say that $\mathbf{u}_{e}(\mathbf{x})$ is the true solution to the underlying Navier-Stokes equations that satisfies $\nabla\cdot\mathbf{u}_{e}(\mathbf{x})=0$. As this exact solution is generally not available, we consider the following three cases

Case (a): the fluid field has interpolation errors $\boldsymbol{\delta}(\mathbf{x})$ that are not divergence free $\mathbf{u}(\mathbf{x})=\mathbf{u}_{e}(\mathbf{x})+\boldsymbol{\delta}(\mathbf{x})$, where $\nabla\cdot\boldsymbol{\delta}(\mathbf{x})\neq 0$ (e.g., using standard polynomial interpolation) 2. 2.

Case (b): the fluid field has interpolation errors $\boldsymbol{\delta}(\mathbf{x})$ and is divergence free $\mathbf{u}(\mathbf{x})=\mathbf{u}_{e}(\mathbf{x})+\boldsymbol{\delta}(\mathbf{x})$, where $\nabla\cdot\boldsymbol{\delta}(\mathbf{x})=0$ (e.g., using MRBF interpolation) 3. 3.

Case (c): the fluid field is free of errors $\mathbf{u}(\mathbf{x})=\mathbf{u}_{e}(\mathbf{x})$ and $\nabla\cdot\mathbf{u}(\mathbf{x})=0$ (e.g., when the velocity field is available in closed form, also referred to as “exact” interpolation.)

We pay particular attention to how errors resulting in $\nabla\cdot\boldsymbol{\delta}(\mathbf{x})\neq 0$ affect how the numerical methods evolve physical volume. To quantify this we define the physical volume error by

[TABLE]

Here, $\Psi$ is used to denote the physical volume over time $h$ of the exact solution with fluid field corresponding to case (c), that is, the physical volume of the true solution in the absence of any errors, whereas $\Psi^{[n]}$ is the physical volume of the numerical solution with fluid field corresponding to one of the three cases given below. The results are presented in table 2. We see here that the physical volume errors in case (a) are $O(h^{2})$ and proportional to $\nabla\cdot\boldsymbol{\delta}(\mathbf{x})$ for all the methods. In case (b), the physical volume errors are proportional to $O(h^{4})$ except for the centrifuge-preserving methods, which have physical volume error proportional to $O(h^{4}\delta_{4})$, where $\delta_{4}=\left|\chi_{4}(\mathbf{u})-\chi_{4}(\mathbf{u}+{\boldsymbol{\delta}})\right|=O(|\boldsymbol{\delta}|)\ll|\chi_{4}(\mathbf{u})|$ is the error of $\chi_{4}$ from the interpolation method, which we assume is small. In case (c), $\delta_{4}=0$ and the centrifuge-preserving methods have physical volume error proportional to $O(h^{5})$, while the other methods are $O(h^{4})$. It is due to this behavior that we expect all the methods to more accurately evolve physical volume, when a divergence-free interpolation method is implemented such as with MRBFs. In this case, we expect the centrifuge-preserving methods to perform especially well due to table 2.

In addition to the erroneous contraction of physical volume, we note from examples 1 - 3 that large divergence errors impose more stringent restrictions on the time step for the numerical methods to be contractive.

6 Numerical simulations

In this section we test our numerical methods for simulating suspensions of particles in viscous flows. The section begins by outlining the flow field and summarizing the methods and numerical parameters. We then outline the computational cost and verify the convergence of the methods. Next we simulate suspensions of $10^{4}$ particles in Taylor-Green vortices. This is the most important part of the section and is comprised of three experiments. The first compares the integration methods with exact evaluation of the fluid field. The second compares the effect of different interpolation errors with the CP2 integration. The third and final experiment explores how a combination of the proposed interpolation and integration methods can be used to generate cost-effective accurate particle distributions compared to conventional methods.

6.1 Preliminaries

Here we will briefly outline the numerical methods that are under consideration in the forthcoming numerical experiments, the fluid field and finally the particle models.

The integration methods under consideration and their properties are summarized in table 3.

We will abbreviate the divergence-free MRBF interpolation with the nearest $(n+1)\times(n+1)\times(n+1)$ data points by MRBF$n$ and the non-divergence-free order $n$ tripolynomial interpolation by TP$n$. The MRBF shape parameters are set to $\epsilon_{1}=0.31$, $\epsilon_{2}=0.23$ and $\epsilon_{3}=0.16$ corresponding to the MRBF1, MRBF2 and MRBF3 schemes, respectively, and are chosen empirically. We will compare the methods against a reference solution that uses exact evaluation of the analytic fluid field and the classical fourth order Runge-Kutta method for time integration with a time step that is 10 times smaller then that of the other methods. Note that such a reference solution is only available in the case that the flow field is known in closed form.

The discrete fluid field is generated by evaluating a closed form solution to the Navier Stokes equation on a regularly spaced grid with uniform sampling in each direction $\Delta x=\Delta y=\Delta z=1/10$. We use a stationary Taylor-Green vortex solution that was proposed in [39] and has been used by other authors to study the behaviour of particles in cellular flow fields [40, 41, 42, 43]. The particular Taylor-Green flow field used in the experiments is given by $\mathbf{u}(\mathbf{x})=(u(\mathbf{x}),v(\mathbf{x}),w(\mathbf{x}))^{T}$ where

[TABLE]

We will perform experiments on both spherical and non-spherical particles. Denoting by $\lambda$ the aspect ratio of the particle, then $\lambda=1$ corresponds to spherical particles, $\lambda>1$ corresponds to a prolate spheroid and $\lambda<1$ corresponds to an oblate spheroid. For $\lambda=1$, the equations of motion are given by equations (1) and (2), while for $\lambda\neq 1$ the equations of motion are (1), (2), (3) and (4). For details about the moment of inertia tensor $J$, torque term $\mathbf{T}$, resistance tensor $K$ for the $\lambda\neq 1$ cases we refer to A. Finally, we note that for all of the following experiments, the particles are given a random initial location within in a box of width $0.01$ centered at the point $x_{0}=(1/3,1/5,1/7)^{T}$ in the domain and a random initial orientation for non-spherical particles.

6.2 Computational cost

Here, we outline the main computational costs associated with the methods. The two main steps in the algorithm are the interpolation step and the time integration step, which we examine separately. The wall clock times $T_{w}$ for $10^{4}$ time steps of the considered integration methods using exact evaluation of the fluid field are measured and presented in table 5 and $T_{w}$ for $10^{4}$ time steps of the various interpolation methods using the FE1 method are presented in table 5.

We note that the centrifuge-preserving methods are slightly more costly due to extra evaluations of the $\Phi_{ah}^{[LT]}$ operator. However, we note that one could speed up many of these splitting methods by observing that they are conjugate to a lower stage faster method, for example

[TABLE]

hence repeated evaluations of the operator $\Phi_{h}^{[SS]}$ when implemented in this way effectively has the same cost as $\Phi_{h}^{[LT]*}$. Similar observations are made for the centrifuge preserving methods.

For the interpolation step, there are two main calculations that contribute the most to the computational cost. The first being the solution of a linear system of size $3n^{3}\times 3n^{3}$ to find the interpolation coefficients. Guassian elimination is used for this purpose due to simplicity and the fact that the systems are not so large (at most $192\times 192$ for the MRBF3 and TP3 methods). However, we note the existence of the exact matrix inverses for the coefficient matrices of these linear systems. This can be found in [44] for the TP method and [45] for the MRBF method, the latter being due to the fact that the coefficient matrix has a block toeplitz structure for MRBF interpolation on Cartesian grids. The next most significant cost is evaluation of the sums of basis functions, that is, the sum in (44) and a similar equation for the TP methods. MRBF interpolation involves evaluation of more complex basis functions (i.e., matrix-vector products containing exponentials of polynomials), which is more costly than evaluating sums of monomials for the TP interpolation. This cost is more of a burden for the MRBF2 and MRBF3 methods as seen in table 5.

We see here that the MRBF1 and TP1 methods are roughly equal in cost. The MRBF2 method is about double that of TP1 and MRBF1 and is more expensive than the TP2 method. The MRBF3 method is double the cost of the TP3 method. We recall that we are not constrained to these three choices of MRBF methods and one is free to use any number of data points to achieve an optimum balance of accuracy and cost. This freedom is due to the fact that MRBF interpolation was designed for interpolation on scattered data points [33]. This option is not available for the TP$n$ methods, where $n+1$ grid points in each dimension are required to ensure the existence of a unique degree $n$ interpolating polynomial.

6.3 Convergence

In this section, we will verify the convergence of the integration methods first with exact interpolation then with various combinations of the interpolation methods for spherical ($\lambda=1$) and non-spherical ($\lambda=10$). In these experiments, we set $St=1$ and compute the particles’ dynamics for time $T=1$.

The convergence of the error, measured in the 2-norm, of the integration methods are presented in figures 1(a) and 1(b). We observe here that the FE1 and LT1 methods have similar accuracy as do the RK2 and SS2 methods. It is noted that the benefits of preserving contractivity in the various splitting methods are expected to be seen after longer times. One remarkable observation here is that for this Stokes number the first order CP1 method is competitive with the second order RK2 and AB2 methods at large time steps, furthermore, the CP2 method is the most accurate by a factor of about 5 in both scenarios.

Figures 1(c) and 1(d) show the convergence of the CP2 and RK2 methods with different interpolation methods. We see that the methods initially converge at their expected order, but as $h$ goes to zero we see that the integration error becomes overshadowed by the $h$-independent interpolation error. We observe that the MRBF solutions yield more accurate solutions than the TP solutions using the equal number of data points. However, the TP3 solution is expected to perform better for longer simulations where particles cross grid cells. This is because the piece-wise fluid field constructed from the TP3 method is globally $C^{1}(\mathbb{R}^{3})$, meaning that the spatial derivatives of the fluid velocity are everywhere continuous. This is not true for the other methods.

Finally, we remark that the centrifuge-preserving methods perform equally well for non-spherical particles.

6.4 Simulating suspensions of particles

Up until now we have mainly focused on the average error in the positions of individual particles. It is well known that standard methods such as polynomial interpolation and Adams-Bashforth integration do a good job at minimizing this truncation error in some norm. However, while it is indeed important that this conventional measure of error is kept at a minimum, in practice one is usually more interested in properties of distributions of many indistinguishable particles meaning that the individual error of each particle is less important. Due to this, it is more desirable that an algorithm reproduces accurately the spatial statistical properties of many particles rather than minimizing the absolute errors of each individual particle. With this in mind, the main goal here is to test to what extent the aforementioned errors affect suspensions of particles when viewed as a single discrete probability distribution. We will show that distributions of particles calculated by our proposed geometric methods will more closely resemble that of the exact solution, despite sometimes having higher average error per particle in the conventional sense.

In our context, a distribution $P=\{(\mathbf{x}_{i},w_{i})\}_{i=1}^{n_{c}}$ is a set of $n_{c}$ non-empty equally sized cells, where $\mathbf{x}_{i}$ is the location of the cell center and $w_{i}$ is a weight that is equal to the number of particles in that cell. We let $P_{n}$ denote a distribution where the particle locations are calculated by a numerical method, $P_{\mathrm{ref}}$ refers to the distribution obtained by the reference solution and we use $300\leq n_{c}\leq 400$ depending on the spread of particles. We will determine the accuracy of $P_{n}$ using three measures which we now outline.

The first is the first Wasserstein distance, which is a natural metric to compare the distance between two discrete probability distributions of equal size (also known as the Earth Mover Distance). The first Wasserstein distance between two probability distributions is denoted by $W_{1}(P_{1},P_{2})$ and is a measure of the cost of transporting the distribution $P_{1}$ into $P_{2}$ in the cheapest way possible. The cost is measured as the distance between cell centers, measured in the 2-norm and weighted by the number of particles being transported. For mathematical details about the first Wasserstein distance, we refer the reader to [46] and the numerical computation of the first Wasserstein distances are computed using a publicly available MATLAB code [47]. We denote by $W_{1}(P_{n})=W_{1}(P_{n},P_{\mathrm{ref}})$ the first Wasserstein distance between $P_{n}$ and $P_{\mathrm{ref}}$.

The second is the relative entropy (also known as the Kullback-Leibler divergence) [48], which is a measure of how much information is lost from a reference distribution $P_{2}$ when an approximate distribution $P_{1}$ is used. The relative entropy is calculated by

[TABLE]

where $P(\mathbf{x}_{i})=w_{i}$ is the number of particles in the cell at $\mathbf{x}_{i}$ and $\Omega_{P}$ is the support of the two distributions. If there is an empty cell in one distribution and not the other, say at $\mathbf{x}=\mathbf{x}_{0}$ we use $P(\mathbf{x}_{0})=10^{-1}$, to avoid singularities. This modestly penalizes the approximate solution for predicting a non-zero probability of having a particle in a cell that should have zero particles according to the reference distribution. We denote by $E(P_{n})=E(P_{n},P_{\mathrm{ref}})/n_{p}$ the relative entropy between $P_{n}$ and $P_{\mathrm{ref}}$ scaled by the number of particles $n_{p}=10^{4}$.

Finally, the third means of determining the accuracy of the distribution is by the average error of the particle positions $\overline{\Delta\mathbf{x}_{n}}$. This conventional measure of error is calculated by taking the difference between the final position of the numerical and reference solution starting from the same initial conditions and averaging over all the $n_{p}=10^{4}$ particles, that is

[TABLE]

where $\mathbf{x}_{n,i}$ is the $i$th particle calculated by the numerical method and $\mathbf{x}_{\mathrm{ref},i}$ is the $i$th particle under the reference solution. As the rotational variables are strongly coupled with the translational variables, errors in the rotational dynamics will also influence the final positions of the particles, hence this is a reasonable measure of the error of the algorithms’ overall accuracy in computing the dynamics of a single particle. We recall that this error is that which the conventional methods are designed to reduce when referring the the global order of accuracy of a method.

In the forthcoming experiments, we will use various combinations of integration and interpolation methods to compute the paths of $10^{4}$ particles in the discrete Taylor-Green vortices starting with random positions and orientations within a cube of width $1/100$ centered at the point $(1/3,1/5,1/7)^{T}$. We perform three experiments. The first compares how the various numerical integration methods and their errors affect the spatial distribution of suspensions of particles in the absence of interpolation errors. The second experiment investigates how interpolation errors affect the spatial distribution of particles using the CP2 method. Finally, we look at how a combination of MRBF interpolation and centrifuge-preserving integration can be used to calculate fast and accurate suspensions of particles compared to the conventional AB2+TP$n$ methods, similar to the methods used in [8, 10, 12, 13, 14], for example.

6.4.1 Comparison of integration methods

In this experiment we use the seven integration methods outlined in table 3 to simulate a suspension of particles evolving in the Taylor-Green flow with exact interpolation. The methods are each tested in six separate simulations, three with Stokes numbers of $St=1/5,1,10$ for spherical particles ($\lambda=1$) and three with the same Stokes numbers for non-spherical particles ($\lambda=1/10$). At the end of the simulation the relative entropy $E(P_{n})$, first Wasserstein distance $W(P_{n})$ and the average spatial error $\overline{\Delta\mathbf{x}_{n}}$ between the numerical distribution and the reference distribution are calculated and presented in table 6. The time step $h$ and total simulation time $T$ are also presented in this table. We start by discussing some qualitative features of the final distributions, examples of which are given in figure 2. We then discuss the results of table 6 in detail.

Figure 2 depicts the final distribution of the particles for the various integration methods. The particle positions are plotted modulo 2 for presentation purposes and represented by black dots, while the reference solution is plotted using green dots. Figures 2(a) to 2(f) correspond to the $St=10$, $\lambda=1$ simulation and is viewed along the $y$ direction. We see here that the CP2 solution is able to predict the correct clustering in all the regions that are predicted by the green reference solution. The LT1, CP1 and SS2 solutions are visually similar to each other, however do not correctly predict clustering of particles in some regions, given by regions of green dots that are void of black dots. The RK2 and AB2 solutions do a worse job as seen, again, by even more regions with a higher concentration of green dots compared to black dots. Similar observations are again seen in figures 2(g) to 2(l), which correspond to the $St=1$, $\lambda=1/10$ simulation, viewed along the $z$ direction. In this simulation, the particles more closely follow the streamlines of the fluid field and more quickly concentrate in regions of high strain as seen by the regions of dense green dots. In these figures, it is even more easily seen that the four contractivity preserving methods do a good job at correctly clustering particles in regions where the reference solution does, while we see with the FE1 and RK2 solutions multiple regions exhibiting an erroneous concentration of black dots that are void of green dots. The AB2 solution is unstable for these parameters.

To quantify the above observations, which have up until now been visual, we turn our attention to table 6. We start by outlining some general observations that are common to all six simulations. Looking first at the order one methods, we observe that the LT1 and CP1 methods, which are contractivity preserving, outperform the FE1 method in almost all measures in each simulation despite the fact that their order of accuracy is the same. What is striking here is that in most simulations the LT1 and CP1 methods generally have a lower $\overline{\Delta\mathbf{x}_{n}}$, $W(P_{n})$ and $E(P_{n})$ compared to the conventional RK2 and AB2 methods despite being of lower order and computational cost. Similar observations are made if we turn our attention towards the order two methods. That is, the SS2 method in most cases has lower $\overline{\Delta\mathbf{x}_{n}}$, $W(P_{n})$ and $E(P_{n})$ than the RK2 and AB2 methods, and better still is the CP2 method. The advantage of the CP2 method over the SS2 is more pronounced than the advantage of the CP1 method over the LT1. The CP2 method has the lowest $\overline{\Delta\mathbf{x}_{n}}$, $W(P_{n})$ and $E(P_{n})$ in all six simulations and is clearly the best method here in all three metrics.

For the $St=1$, $\lambda=1/10$ simulation, the CP1 solution has a larger $\overline{\Delta\mathbf{x}_{n}}$ than the SS2 solution, but a lower $W(P_{n})$ and $E(P_{n})$ which suggests that for these simulation parameters, the centrifuge-preserving property is more advantageous for producing more accurate distributions than simply reducing the accuracy of the method in the conventional sense. In this simulation, the CP1 method is the second best in all measures, the best being the CP2 method. It is also noteworthy that in other simulations the CP1 method has roughly equal, and sometimes lower $\overline{\Delta\mathbf{x}_{n}}$, $W(P_{n})$ and $E(P_{n})$ than the SS2 solution, which further suggests that the centrifuge-preserving property is advantageous.

One of the most remarkable observations is made for the $St=10$, $\lambda=1/10$ experiment. Here, the values of $\overline{\Delta\mathbf{x}_{n}}$ are quite severe and roughly the same for all methods, due to the fact that the time step is quite large and the non-spherical ODEs are more stiff. Despite this, the contractivity preserving methods have a much lower $W(P_{n})$ and $E(P_{n})$ and the centrifuge preserving methods are better still. This highlights the fact that preserving the aforementioned physical features in the numerical solution results in spatial distributions that more closely resemble the reference solution, despite having the same $\overline{\Delta\mathbf{x}_{n}}$.

Finally, we mention that all the splitting schemes have better stability properties and are still able to produce accurate clusters of particles for low Stokes numbers and reasonably large time steps as noted by the $\lambda=1/10$ simulations for $St=1/5$ and $St=1$ where we begin to see some of the conventional methods losing stability.

6.4.2 Comparison of interpolation methods

In this experiment we compare the MRBF and TP interpolation methods in combination with the CP2 method to simulate a suspension of particles evolving in the Taylor-Green flow. Six separate simulations are performed, three with Stokes numbers of $St=1/10,1,10$ for spherical particles ($\lambda=1$) and three with the same Stokes numbers for non-spherical particles ($\lambda=5$). At the end of each simulation the average spatial error $\overline{\Delta\mathbf{x}_{n}}$, relative entropy $E(P_{n})$ and the first Wasserstein distance $W(P_{n})$ between the numerical distribution and the reference distribution are calculated and the results are presented in table 7. The time step $h$ and total simulation time $T$ are also presented in this table. Some spatial distributions produced by the different interpolation methods are presented in figure 3.

Directing our attention towards figures 3(a) to 3(f), which show the final distribution of the $St=1/10$, $\lambda=5$ simulation looking down the $z$-axis. It can be seen here that the three MRBF solutions look visually very similar to the reference solution, as does the TP3 solution. If we look towards the corresponding part of table 7, we see that the TP3 solution has a $\overline{\Delta\mathbf{x}_{n}}$ of 0.3468, which is lower than the MRBF1 solution, which has a $\overline{\Delta\mathbf{x}_{n}}$ of 0.5389. Despite this, the MRBF1 solution, which we note performs exceptionally well here, has a lower $E(P_{n})$ and $W(P_{n})$ meaning that the final distribution is more similar to the reference distribution even though the $\overline{\Delta\mathbf{x}_{n}}$ is greater.

Figures 3(g) to 3(l) show the final distribution of the $St=1$, $\lambda=1$ simulation looking down the $x$-axis. Here we see that the TP1, TP2 and MRBF1 solutions differ visually from the reference solution, whilst the TP3, MRBF2 and MRBF3 solutions look quite similar. From the corresponding section of table 7, we see that the MRBF2 solution’s $\overline{\Delta\mathbf{x}_{n}}$ is 0.5004 compared to the TP3 solution, which is 0.3332, but both have a similar $E(P_{n})$ and $W(P_{n})$. Furthermore, there are many examples here of the MRBF solutions having higher $\overline{\Delta\mathbf{x}_{n}}$, but lower $E(P_{n})$ and $W(P_{n})$. This can be seen in all three $\lambda=1$ simulations, where the MRBF2 solution has larger $\overline{\Delta\mathbf{x}_{n}}$ than the TP3 solution but similar or lower $E(P_{n})$ and $W(P_{n})$. For the $\lambda=5$ simulations and for $St=1$ and $St=10$, the MRBF2 solution outperforms the TP3 solution in all three measures. Such examples indicate that preserving the divergence-free condition is important to acheive accurate spatial distributions.

We now make some general observations about table 7. We see that in all but one simulation, the MRBF1 solutions outperform the TP2 solution in all three measures. It is noteworthy that the MRBF1 solution is as fast as TP1 interpolation where both require only eight data points for the interpolation as opposed to 27 data points for the TP2 interpolation, which is a slower method. Additionally, in all six simulations the MRBF3 interpolation method outperforms all the TP solutions in all measures.

To summarize, we have seen many examples of the MRBF solutions producing distributions that are more similar to the reference solution than the TP solutions, despite having worse $\overline{\Delta\mathbf{x}_{n}}$. These observations are consistent with the fact that the CP2 method, among others, loses accuracy in $\Delta\Psi^{[n]}$ when evolving particles in a non-divergence-free flow field. That is, the physical volume $\Psi^{[n]}$ is more strongly affected when the divergence-free condition is broken, despite the fact that the order of accuracy of the method remains unaffected.

6.5 Comparison of interpolation and integration methods

Our final experiment explores the benefit that is gained by combining MRBF interpolation with the centrifuge- and contractivity-preserving methods compared to the standard methods used in the literature. We will compare the methods TP1+FE1, MRBF1+CP1, TP2+AB2, TP3+AB2, MRBF2+CP2 and TP2+CP2. The first two methods are the cheapest and are of roughly equal cost. The method TP2+AB$n$ are used in, for example [8, 10, 12, 13, 14] and subsequent studies. We include TP3+AB2 to test whether increasing the interpolation accuracy is worthwhile use of computational resources. We also consider the MRBF2+CP2 solution, which is an accurate and economical combination of our proposed geometric methods. Finally, the TP2+CP2 method is considered to emphasize the negative implications of using a non-divergence-free interpolation method with the CP2 method.

Six simulations are performed, three with Stokes numbers of $St=1/10,1,10$ for spherical particles ($\lambda=1$) and three with the same Stokes numbers for non-spherical particles ($\lambda=10$). At the end of the simulation, the average spatial error $\overline{\Delta\mathbf{x}_{n}}$, relative entropy $E(P_{n})$ and the first Wasserstein distance $W(P_{n})$ between the numerical distribution and the reference distribution are computed and presented in table 7 along with the time step and total simulation times used.

Figures 4(a) to 4(f) show the final distribution of the $St=1/10$, $\lambda=10$ simulation looking down the $z$-axis. We see that the only methods that look similar to the green reference solution are the MRBF1+CP1, MRBF2+CP2 and TP3+AB2 solutions. That is, both of our geometric methods and the most costly conventional method. It is worth noting that the TP2+CP2 and the TP2+AB2 solutions look very similar, despite one being generated by the CP2 method. This is likely due to the fact that the CP2 method loses its centrifuge-preserving properties when the fluid field is not divergence-free as seen in table 2 as well as the interpolation method being the dominant source of error. Turning our attention towards the corresponding part of table 8 (i.e., the top-right) we make a few remarks. The most striking one here is that the MRBF1+CP1 solution has a larger $\overline{\Delta\mathbf{x}_{n}}$ than the TP3+AB2 solution, but its $E(P_{n})$ and $W(P_{n})$ are both far lower. This is in agreement with the fact that the TP3+AB2 method is better in reducing error in the conventional sense, (i.e., the average 2-norm of the particle position errors $\overline{\Delta\mathbf{x}_{n}}$), but does a poorer job at reproducing the mechanisms that are responsible for the preferential distribution of particles (i.e., the centrifuge effect and the sum of the Lyapunov spectrum). We can make similar remarks for the $St=1/10$, $\lambda=1$ experiment where, despite having higher $\overline{\Delta\mathbf{x}_{n}}$, the MRBF1+CP1 method has a lower $E(P_{n})$ and $W(P_{n})$ than the TP3+AB2 method. These significant observations are prevalent in both $St<1$ experiments, which corresponds to a more stiff vector field with greater relative influence of fluid inertia.

We finish with some general observations. First, the MRBF1+CP2 solutions outperform the TP1+FE1, TP2+AB2 and TP2+CP2 methods in all three measures. Second, the TP2+AB2 and TP2+CP2 methods perform about the same in each simulation, suggesting there is little advantage from using the CP2 method in conjunction with the TP2 solution. Finally, we note that the MRBF2+CP2 method is more accurate than the TP2+AB2 and the TP3+AB3 in all cases. The only exception is the $St=1$, $\lambda=1$ simulation where the MRBF2+CP2 method has worse $\overline{\Delta\mathbf{x}_{n}}$ than the TP3+AB3 solution.

7 Conclusions

A novel combination of geometric numerical methods for calculating accurate distributions of inertial particles in viscous flows is proposed. The algorithm consists of MRBFs to construct a divergence-free approximation of the background flow field and a geometric splitting method for the time integration.

The splitting method is shown to preserve the sum of the Lyapunov spectrum and hence the contractivity of phase space volume. By expanding the exact solution we derive an expression for how a physical volume of particles change over a small time step $h$, with which we recover the centrifuge effect at $O(h^{4})$. We show that when a divergence-free interpolation method is used, one can implement a so-called centrifuge-preserving splitting method that preserves not only the qualitative but also the quantitative behavior of this centrifuge effect. Moreover, it is shown that errors to the divergence of the fluid field can overshadow this effect when a conventional polynomial interpolation method is used, for example.

It is shown through numerical experiments that MRBF interpolation yields particle distributions that are more similar to the exact solution than standard TP interpolation. In many examples, this is observed even when the MRBF solution has higher error per particle. This is, in part, explained by the fact that: (1) MRBF interpolation is divergence-free meaning that the numerical time integration methods mimic the qualitative centrifuge effect; and (2) MRBF interpolation produces a vector field that solves the Stokes equations, meaning that the background flow field more physically resembles that of the exact solution (e.g., the flow field is related to the gradient of a scalar pressure function).

Furthermore, we see that the proposed centrifuge-preserving methods are superior to the standard methods in terms of error per particle and how closely the particle distribution resembles the exact distribution. This is true, remarkably so, even when comparing the order-one CP1 method to the order-two RK2 and AB2 methods. In particular, for experiments with low particle inertia, the MRBF1+CP1 method produces more accurate distributions of particles than the expensive TP3+AB2 solutions, despite having slightly worse error and being an order one integration method that uses far less data points for the interpolation step.

These observations strongly suggest that preserving certain physical features of ODEs under study in the numerical solution is of importance when simulating inertial particles in discrete flow fields. Of particular interest for future studies would be to implement the proposed methods in a physically realistic flow fields generated by a direct numerical simulation of homogeneous isotropic turbulence or turbulent channel flow, for example.

8 Acknowledgments

This work has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement (No. 691070). B. K. Tapley, E. Celledoni and B. Owren would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Geometry, compatibility and structure preservation in computational differential equations (2019) where part of the work on this paper was undertaken.

Appendix A Non-spherical particle model

Here, we give details of the specific rigid spheroid model that is used in the numerical experiments. The surface of a spheroid is defined by the equation

[TABLE]

where $a$ and $c$ are the distinct semi-axis lengths. The particle shape is characterised by the dimensionless aspect ratio $\lambda=c/a>0$, which distinguishes between spherical ($\lambda=1$), prolate ($\lambda>1$) and oblate ($\lambda<1$) particles (the latter two shapes are also called as rods and disks).

An inertial particle immersed in a fluid will experience forces on its surface that have magnitude governed by many parameters such as the particles density $\rho_{p}$, length $a$, fluid density $\rho_{f}$, kinematic viscosity $\nu$ and fluid time scale $\tau_{f}$. The particle Stokes number is formally defined as the ratio of the particle and fluid time scales $St=\tau_{p}/\tau_{f}$. For a spherical particle the Stokes number is

[TABLE]

where $D=\rho_{p}/\rho_{f}$ is the particle-fluid density ratio. Note that this definition only depends on the particle size and inertia. For spheroidal particles, the following shape dependent Stokes numbers are used, which are derived by Shapiro and Goldenberg [49] and Zhao, et al. [50]

[TABLE]

where $k_{0}=\log((\lambda-\sqrt{\lambda^{2}-1})/(\lambda+\sqrt{\lambda^{2}-1}))$. Note that $St\rightarrow St_{0}$ as $\lambda\rightarrow 1$ from above or below. All the following equations are implemented in their non-dimensional form and all parameters have dimension equal to $1$.

The particle experiences a hydrodynamic drag force due to Brenner [51],

[TABLE]

where $\mathbf{u}=\mathbf{u}(\mathbf{x},t)$ is the fluid velocity evaluated at the particle center of mass $\mathbf{x}$ and $\mathbf{v}=\mathbf{p}/m$ is the particle velocity. The body frame resistance tensor $K_{b}$ was calculated by Oberbeck [52], is diagonal, positive definite and given by

[TABLE]

where the constants $\chi_{0}$, $\alpha_{0}$, $\beta_{0}$ and $\gamma_{0}$ were calculated for ellipsoidal particles by Siewert et al. [53] and are presented in table 9

The torque vector $\mathbf{T}$ depends on the particle shape and the local fluid velocity derivatives, and is given in non-dimensional form by [54]

[TABLE]