Numerical simulation of particulate flows using a hybrid of finite difference and boundary integral methods

A combination of finite difference (FD) and boundary integral (BI) methods is used to formulate an efficient solver for simulating unsteady Stokes flow around particles. The two-dimensional (2D) unsteady Stokes equation is being solved on a Cartesian grid using a second order FD method, while the 2D steady Stokes equation is being solved near the particle using BI method. The two methods are coupled within the viscous boundary layer, a few FD grid cells away from the particle, where solutions from both FD and BI methods are valid. We demonstrate that this hybrid method can be used to accurately solve for the flow around particles with irregular shapes, even though radius of curvature of the particle surface is not resolved by the FD grid. For dilute particle concentrations, we construct a virtual envelope around each particle and solve the BI problem for the flow field located between the envelope and the particle. The BI solver provides velocity boundary condition to the FD solver at “boundary” nodes located on the FD grid, adjacent to the particles, while the FD solver provides the velocity boundary condition to the BI solver at points located on the envelope. The coupling between FD method and BI method is implicit at every time step. This method allows us to formulate an $O\left(N\right)$ scheme for dilute suspensions, where $N$ is the number of particles. For semidilute suspensions, where particles may cluster, an envelope formation method has been formulated and implemented, which enables solving the BI problem for each individual particle cluster, allowing efficient simulation of hydrodynamic interaction between particles even when they are in close proximity. The method has been validated against analytical results for flow around a periodic array of cylinders and for Jeffrey orbit of a moving ellipse in shear flow. Simulation of multiple force-free irregular shaped particles in the presence of shear in a 2D slit flow has been conducted to demonstrate the robustness of the FD-BI scheme. This method can also be used to solve the full Navier-Stokes equations around particles. In this case, it is shown that, below a threshold particle Reynolds number, which in turn depends on the required resolution of the particle surface, the FD-BI scheme will be more computationally efficient than a pure FD scheme solved on a multiblock grid.

Figure 1

(a) Schematic of hybrid FD-BI method, shown for a single particle in a flow domain. For a staggered grid arrangement, the nodes correspond to cell centers of either the $u_1$ grid or $u_2$ grid. (b) Schematic of hybrid FD-BI method, shown for multiple particles in the flow domain. When particle centers are within a distance of $2R_O$, a cluster is formed, with a common outer envelope ${\mathrm{\Gamma }}^{O,1}$.

Figure 2

Schematic for arrangement of pressure nodes and velocity “boundary” nodes around particle. $(\to )$ denotes location of $u_1$ velocity boundary nodes. $(\uparrow )$ denotes location of $u_2$ velocity boundary nodes. $\left(▵\right),(\square )$, and $(◯)$ denote fluid, solid, and boundary pressure nodes, respectively. The shaded region shows the circle of radius $R$ that encloses the particle. The dotted line (circle centered with particle, radius $R+\epsilon$) demarcates solid and boundary velocity nodes. Note that fluid and solid velocity nodes have not been marked here.

Figure 3

(a) Schematic for envelope evolution. The envelopes represent either ${\mathrm{\Gamma }}^O$ or ${\mathrm{\Gamma }}^{\epsilon }$ for each cluster [e.g., Fig. 1]. Envelope shape evolves from $E1$ (right after cluster is formed) to $E3$. Harmonic springs in envelope keep it taut, while repulsive Morse potential from particle ensures that the envelope stays at least a distance $\mathrm{\Delta }{R}_c$ from the particle surface. (b) Schematic showing balance of pressure $\mathrm{\Delta }P$ and tension in envelope ${\sigma }_{\mathrm{avg}}$ (from inset on concave segment of the envelope). The forces on the particle arise purely due to the traction from the fluid and any additional exclusion forces between particles surfaces within a cluster; the envelope does not apply any force on the particle.

Figure 4

(a) Computational setup for convergence test. The two black circles denote the inner and outer boundaries of the computational domain. The red circle shows ${\mathrm{\Gamma }}^O$. The blue arrows denote velocity field in the fluid and boundary velocity nodes. (b) Traction component $f_x$ plotted with respect to angle $\theta =s^P/R$ on cylinder surface for different $\mathrm{\Delta }$ values, for fixed $R^O$ (set A). (c) Traction component $f_x$ versus angle $\theta$ for different $\mathrm{\Delta }$ values, for $R^O=R+\epsilon +3.25\mathrm{\Delta }$ (set B).

Figure 5

(a) Computational setup for convergence test for the wavy particle. The red circle shows ${\mathrm{\Gamma }}^O$. The blue arrows denote velocity field in the fluid and boundary velocity nodes. (b) Traction component $f_x$ plotted with respect to angle $\theta ={tan}^{-1}(y/x)$ on particle surface, for different $\mathrm{\Delta }$ values, for fixed $R^O$ (set C). (c) Traction component $f_x$ versus angle $\theta$ for $R^O=R+\epsilon +3.25\mathrm{\Delta }$ (set D) for different values of $\mathrm{\Delta }$.

Figure 6

(a) Relative error (based on $L^2$ norm) in traction ${\mathbf{f}}^P\left(s^P\right)$, between simulation and reference, plotted with respect to $\mathrm{\Delta }$, for sets with $R^O$ fixed w.r.t. $\mathrm{\Delta }$, for particles with smooth (set A) and wavy (set C) surfaces. (c) Relative error (based on $L^2$ norm) in traction ${\mathbf{f}}^P\left(s^P\right)$, between simulation and reference, plotted with respect to $\mathrm{\Delta }$, for sets with $R^O=R+\epsilon +3.25\mathrm{\Delta }$ for particles with smooth (set B) and wavy (set D) surfaces.

Figure 7

Evolution of net force on particle $F_y$ as a function of time, for stationary and moving particle. $R=3$ for all cases. The other parameters for the simulation are (a) $\mathrm{\Delta }=0.5,R^O=5.125$ (b) $\mathrm{\Delta }=2.0,R^O=11.5$.

Figure 8

(a) Vector plot (blue arrows) of flow around an array of cylinders. Red circles indicate ${\mathrm{\Gamma }}^O$ for each cylinder (black circle). (b) Comparison of $F_x$ vs $c$ between current simulations and prior analysis with analytical solutions by [31] (in dilute limit) and [32] (in the dense limit).

Figure 9

(a) Snapshot of velocity field around a force-free elliptical particle (black line) in shear flow, with $b/a=0.25$. The uniform background shear flow has been subtracted from the total flow. (b) Snapshot of velocity field around a force-free elliptical particle (black line) in shear flow, with $b/a=0.75$. The red circle again denotes ${\mathrm{\Gamma }}^O$. (c) Comparison of analytical versus numerical results for angular velocity of ellipse as a function of time, for $b/a=0.75$ (upper) and $b/a=0.25$ (lower).

Figure 10

CPU time versus number of particles (solid line with marker). Curve indicating $N$ scaling has been drawn for reference (dashed line). Size of domain is proportional to number of particles.

Figure 11

(a) Schematic of problem for studying the advantage of clustering. (b) Variation of nondimensional $x$ component of traction over boundary of left particle $\delta =1.4$. (c) Variation of computational time (in seconds) w.r.t. $\delta$.

Figure 12

Breaking of a particle cluster (a) $t=720$, (b) $t=1128$, (c) $t=1200$. Magenta (solid) and green (dotted) curves around the particles denote the outer envelopes ${\mathrm{\Gamma }}^O$ and inner envelopes ${\mathrm{\Gamma }}^{\epsilon }$, respectively. The tagging of the FD $u$-velocity nodes is shown as follows: the blue squares, black asterisks, and red triangles denote the fluid, boundary, and solid $u$-velocity nodes, respectively.

Figure 13

Formation of a particle cluster (a) $t=4032$, (b) $t=4224$, (c) $t=4560$. Other details are same as in Fig. 12.

Figure 14

(a) Flow field and particle configuration at a time instant for force-free particles in confined shear flow. (b) Zoomed in flow field at a time instant for force-free particles in confined shear flow. Background shear flow has been subtracted to highlight flow field caused by the particles.

Figure 15

Variation of effective viscosity w.r.t. nondimensional time for force-free particles in confined shear flow.