Rotation matters - Direct numerical simulations of rectangular particles in suspensions at low to intermediate solid fraction

TL;DR

This paper introduces a novel direct numerical simulation method for rectangular particles in suspensions, enabling detailed analysis of how particle rotation influences suspension dynamics at low to intermediate solid fractions.

Contribution

The study develops a new numerical approach combining Lagrange multipliers and immersed boundary methods for simulating rectangular particles, including an efficient collision scheme.

Findings

Particle rotations significantly influence suspension behavior.

The numerical method accurately captures particle interactions and rotations.

Rotational energy remains small compared to translational energy but still impacts dynamics.

Abstract

Our ability to numerically model and understand the complex flow behavior of solid-bearing suspensions has increased significantly over the last couple of years, partly due to direct numerical simulations that compute flow around individual interfaces and hence resolve unprecedented detail. While most previous studies focus on spherical particles, we develop a direct numerical approach to capture rectangular particles. Our approach uses distributed Lagrange multipliers to enforce rigid-body motion in the solid domain in combination with an immersed boundary method to correctly enforce the no-slip constraint on the solid-fluid interfaces. An important component of our model is an efficient particle collision scheme that prevents overlap between particles of different shapes and allows for the transfer of both translational and angular momentum during particle collision. We verify and validate our numerical method through several benchmark cases. Applied to suspension flow, we test the hypothesize that particle rotations alter the aggregated dynamics of the suspension even if the relative rotational energy of the particles remains small as compared to the translational energy.

Equations

$$\nabla \cdot \mathbf{v}=0,$$

$$-\frac{\nabla p}{\rho }+\frac{\mu }{\rho }{\nabla }^2\mathbf{v}+\mathbf{g}=0,$$

$$\frac{M_{p}d{\mathbf{V}}_p}{dt}$$

$$\frac{d({\mathbf{I}}_p{\omega }_p)}{dt}$$

$$\frac{d{\mathbf{X}}_p}{dt}$$

$${\mathbf{I}}_p=((2l_1{)}^2+(2l_2{)}^2)M_p/12.$$

$${\chi }_q={\mathbf{V}}_p+{\omega }_p\times {\mathbf{r}}_{pq}.$$

$$-r+\sqrt{(x_D-x_O{)}^2+(y_D-y_O{)}^2}<0,$$

$$x_B=\sqrt{r^2-(y_D-y_O{)}^2}+x_O,$$

$$y_C=\sqrt{r^2-(x_D-x_O{)}^2}+y_O,$$

$$\displaystyle A_{p}=A_{BDC}+A_{\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{BOC}}}-A_{BOC}\ ,$$

$$\displaystyle A_{\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{BOC}}}=\frac{\alpha}{2\pi}\times 2\pi r\ ,$$

$$\begin{array}{r}{l}_1^a{x}_D+l_1^b{y}_D+l_1^c<0,\\ l_2^a{x}_D+l_2^b{y}_D+l_2^c>0,\\ l_3^a{x}_D+l_3^b{y}_D+l_3^c>0,\\ l_4^a{x}_D+l_4^b{y}_D+l_4^c<0,\end{array}$$

$$x_B=\frac{-l_3^b{y}_D-l_3^c}{l_3^a},$$

$$y_C=\frac{-l_3^a{x}_D-l_3^c}{l_3^b},$$

$$-\frac{\nabla p}{\rho (\mathbf{x})}+\frac{\mu }{\rho (\mathbf{x})}{\nabla }^2\mathbf{v}+\mathbf{g}+{\mathbf{f}}_{DLM}+{\mathbf{f}}_{IB}=0,$$

$$\nabla \cdot \mathbf{v}=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}=-\frac{\nabla p}{\rho }+\frac{\mu }{\rho }{\nabla }^2\mathbf{v}+\mathbf{g},$$

$$\frac{\hat{\mathbf{v}}-{\mathbf{v}}^n}{\Delta t}=\frac{\mu }{\rho (\mathbf{x})}{\nabla }^2\hat{\mathbf{v}}+\mathbf{g}+{\mathbf{f}}_{IB}^n\text{in}\Omega ,$$

$$\hat{\mathbf{v}}={\mathbf{b}}^n\text{on}\Gamma ,$$

$$\begin{array}{r}\left[1-\frac{\mu \Delta t}{\rho (\mathbf{x})}\frac{{\partial }^2}{\partial x^2}-\frac{\mu \Delta t}{\rho (\mathbf{x})}\frac{{\partial }^2}{\partial y^2}\right]\hat{\mathbf{v}}={\mathbf{v}}^n+\mathbf{g}+{\mathbf{f}}_{IB}^n.\end{array}$$

$$\begin{array}{rl}& \left[1-\frac{\mu \Delta t}{\rho (\mathbf{x})}\frac{{\partial }^2}{\partial x^2}\right]\left[1-\frac{\mu \Delta t}{\rho (\mathbf{x})}\frac{{\partial }^2}{\partial y^2}\right]\hat{\mathbf{v}}={\mathbf{v}}^n+\mathbf{g}+{\mathbf{f}}_{IB}^n,\end{array}$$

$$\hat{\mathbf{v}}={\mathbf{b}}^n+\Delta t\frac{1}{\rho (\mathbf{x})}\nabla p^n.$$

$$\frac{{\mathbf{v}}^{n+1}-\hat{\mathbf{v}}}{\Delta t}=-\frac{1}{\rho (\mathbf{x})}\nabla p^{n+1},$$

$$\nabla \cdot {\mathbf{v}}^{n+1}=0\text{in}\Omega ,$$

$$\mathbf{n}\cdot {\mathbf{v}}^{n+1}=\mathbf{n}\cdot {\mathbf{b}}^{n+1}\text{on}\Gamma ,$$

$$\nabla \cdot (\frac{1}{\rho (\mathbf{x})}\nabla p^{n+1})=\frac{1}{\Delta t}\nabla \cdot \hat{\mathbf{v}}\text{in}\Omega ,$$

$$\nabla p^{n+1}\cdot \mathbf{n}=0\text{on}\Gamma .$$

$$\mathbf{D}=\frac{1}{2}\left[\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^{\text{T}}\right]=0,$$

$${\mathbf{V}}_p=\frac{1}{M_p}{\int }_p{\rho }_p{\mathbf{v}}^{\ast }dA\approx \frac{{\rho }_p}{M_p}\sum _{i,j}{\Phi }_{i,j}{\mathbf{v}}_{i,j}^{\ast },$$

$$\begin{array}{rl}{\omega }_p& =\frac{1}{{\mathbf{I}}_p}{\int }_p\mathbf{r}\times {\rho }_p{\mathbf{v}}^{\ast }dA\\ & \approx \frac{1}{{\mathbf{I}}_p}\sum _{i,j}{\Phi }_{i,j}{\mathbf{r}}_{i,j}\times {\rho }_p{\mathbf{v}}_{i,j}^{\ast }.\end{array}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

The behavior of small particles, drops, or bubbles in a viscous fluid at low Reynolds number is one of the oldest problems in theoretical fluid mechanics dating back as far as Stokes’ analysis of a rigid sphere in an unbounded fluid from 1851 [1]. Nonetheless, the flow dynamics of solid-bearing suspensions remain enigmatic [2], despite decades of research progress as summarized in [3]. One key challenge is that long-range interactions between particles lead to surprising spatial correlations in velocity fluctuations even at low solid fraction [4, 5]. The resulting stochastic behavior is reminiscent of turbulence, even though the Reynolds number is very low [6, 7, 8].

Despite the analogy between solid-bearing suspensions at low Reynolds number and turbulent flows, different modeling approaches are required for the two flows. In turbulent fluids, stochasticity emerges from dynamic instability in the fluid phase. Superimposing the motion of small, spherical particles on the fluid motion in a Lagrangian framework, which tracks particles but does not couple these back to the fluid phase, captures at least some of the important components of the aggregate dynamics [e.g., 9, 10]. In contrast, stochasticity in solid-bearing suspensions at low Reynolds number emerges from the coupling between the solid and the fluid phase. A Lagrangian approach is hence questionable and a fully-coupled Eulerian model, where flow transports particles and particles reroute the flow around them and interact, is required.

Direct numerical simulations of particulate flows solve the Stokes or Navier-Stokes equation at the scale of individual particles and offer a compelling, if computationally expensive, opportunity to compute the long-range hydrodynamic interactions without a-priori assumptions about drag forces or settling behavior. The accuracy of the numerical description, however, becomes delicate, which is not surprising for a highly nonlinear, multi-body boundary value problem. It is not sufficient to ensure rigid-body motion for solid particles [e.g., 11, 12]. It is equally important to enforce the no-slip boundary condition on the solid interfaces correctly [e.g., 13] and represent particle-particle collisions in a meaningful way [e.g., 14, 15, 16].

To simplify the problem, most previous approaches have focused on spherical particles [e.g., 11, 12, 17, 18]. While clearly a meaningful first step, analytical descriptions of the motion of single, non-spherical particles through singularity methods [19, 20, 21, 22] or Euler angles [23, 24] have shown that non-spherical particles are more likely to rotate than spheres. The reason is that the area of the particle which is exposed to mean flow depends on the angle of rotation. As a consequence, the hydrodynamic drag force acting on the particle depends on the angle between particle and mean flow, as first noted by Becker [25]. Non-spherical shapes hence introduce an additional physical effect, which is largely absent in the spherical limit.

Here, we hypothesize that the tendency of non-spheroidal particles to rotate affects suspension dynamics sensitively. To test this hypothesis, we develop a numerical methodology to better understand the dynamics of non-spherical particles in suspensions at low to intermediate solid fraction. We use direct numerical simulations to capture different crystal shapes building on [26, 27, 28]. The two main advances of our approach are (1) the usage of an immersed boundary method (IBM) following Qin and Suckale [29] to improve solver accuracy in the vicinity of the solid-fluid interfaces and (2) the development of an efficient collision scheme between particles of arbitrary shape, which is necessary to be able to capture larger particle assemblages of non-spherical particles than Ardekani et al., [27].

The solid-fluid coupling we implement takes advantage of the fact that momentum is conserved in both the solid and fluid phase. Solving the Navier-Stokes equation in the entire domain hence assigns the correct momenta to the phases, but an additional constraint is needed to infer the velocity field in the solid phase from the computed momenta. Following Glowinski et al., [11], we use Lagrange multipliers distributed throughout the solid domain to project the fluid velocity field onto rigid body motion by enforcing that the deformation-rate tensor is zero there. Since the distributed Lagrange multipliers are only defined in the solid domain, they are not ideally suited for enforcing the correct boundary condition in the fluid adjacent to the solid boundary.

To correct the velocity field in the fluid at the solid-liquid interface, we use an immersed boundary method. A number of immersed boundary methods have been developed in the past two decades as reviewed in Mittal and Iaccarino, [30]. In the present work, we implement a direct forcing strategy which is applicable to boundaries of arbitrary shape [e.g., 13] and is more efficient than other approaches such as the feedback forcing strategy [31]. In our verification section, we take advantage of the analytical solution for Stokes flow around a circular cylinder in the vicinity of a moving wall to demonstrate that the proposed immersed boundary method affords higher order accuracy than previous direct forcing schemes and models that only rely on distributed Lagrange Multipliers to correct the solid-fluid boundaries [e.g., 27, 17].

To avoid unphysical overlap of the particles, we develop an efficient collision model for rectangular particles and mixed spherical and rectangular particles. Most previous studies of the behavior of non-spherical particles in suspensions assume spheroidal particles that have an axis of rotational symmetry [e.g., 32, 33, 34]. We focus specifically on rectangular particles, because the development of this numerical technique was motivated specifically by the need to better understand hazardous lava flows. The minerals that form in lava flows during cooling typically have polygonal shapes (see Fig. 1).

The dynamics of lava flows [e.g., 37] and magma reservoirs [e.g., 38] depends on jamming phenomena induced by crystal clustering, which arise not only from the total crystal content but also from the long, elongated shape of some of the crystals such as plagioclase, which contribute a significant portion of the crystal load (e.g., Fig. 1A). Not surprisingly, existing laboratory work suggests that the presence of plagioclase crystals produces jamming at very low crystal fraction [e.g., 39, 40]. The ability to resolve sharp, rectangular shapes is hence critical to advance our ability to model and understand lava flows.

Capturing clustering and jamming of suspended particles depends not only on resolving fluid flow around the solid-fluid boundary, but also requires a shape-adjusted collision or contact model. One commonly adopted collision scheme pioneered by [11] applies a repulsive force shortly before particles come in contact. The repulsive potential is set up such that particles never come in direct physical contact or overlap. The scheme was improved [41] by updating the particle velocity implicitly to avoid numerical instability and an overly restrictive Courant-Friedrichs-Lewy condition.

While the repulsive-potential scheme is common in models of dilute flows [e.g., 42, 43, 18, 44, 45], it becomes questionable in the limit of densely packed sediment beds, where particle-particle contact is ubiquitous [14, 46, 16]. The two main criticisms are that repulsive-potential methods introduce an artificial gap between particles and that they do not resolve the physics of collision and rebound, which is particularly important for understanding the mobilization of dense sediment beds [47, 48].

In the limit of lava flows, the limitations are less consequential than for densely packed sediment beds. Microtextural analyses of crystal-bearing lava flows show that crystals often tend to be separated from each other by a finite melt film (e.g., Fig. 1A). Introducing a small but finite separation between particles into the numerical model is hence not unrealistic. In cases where this separation does not exist or is less clearly pronounced (e.g., Fig. 1B), the main mechanism contributing to clustering is that suspended crystals provide preferred nucleation sites, which leads to inter-growth of crystals. We neglect these geochemical effects here in the interest of simplicity.

The physics of rebound also play a lesser role in magmatic flows than in sedimentary beds, because the shear viscosity of magmatic melts tends to be at least five orders of magnitude more viscous than the shear viscosity of water and often much higher than that. The impact Stokes number, St = $\rho_{p}v_{0}D/(9\mu_{f})$, where $\rho_{p}$ represent the particle density, $D$ the particle diameter, $\mu_{f}$ the pure fluid viscosity, and $v_{0}$ the relative approach speed of the colliding particles, which governs rebound behavior, is hence typically much smaller than zero. For the validation cases that we consider later, St falls into the $10^{-6}$ to $10^{-7}$ range and can hence be neglected. Suppressed rebound is also consistent with the enduring contact of the particles observed in many field samples (e.g., Fig. 1).

We hence argue that in the context of lava-flow modeling, the repulsive potential model proposed by Ardekani and Rangel [41] provides a meaningful starting point. Here, we extend their scheme to particles with both spherical and rectangular shapes. To ensure that the repulsive forces prevent inter-penetration and satisfy the laws of Newtonian dynamics, we apply the analytical model proposed by Baraff [49] to formulate a linear system of inequality and equality constraints on the repulsive forces. The main advantage of this model is that it reduces the computation of the forces between systems of rigid bodies in collision to the solution of a linear programming problem which can be solved through a standard sparsity exploiting linear programming package, here we use the UMFPACK library written by Timothy Davis ([50]).

We have employed five benchmark cases to verify and validate our numerical model. To verify the numerical efficiency and mathematical accuracy of our solid-fluid coupling scheme, we take advantage of analytic solutions including the Wannier flow [51] and the flow field in a rotational viscometer [e.g., 52, 53]. To validate our solid-fluid coupling at finite Reynolds number, we compare our simulations to experimental results for flow past a square cylinder [54]. We then test whether our collision model is consistent with prior numerical results of spherical particles in a Couette device [55, 52] and with experimental results of rectangular particles aggregating while sinking through a stagnant fluid [56].

2 Governing Equations

Our numerical model solves for conservation of mass and momentum in two dimensions. We assume that the pure liquid phase is incompressible. In the liquid domain, the governing equations are hence the incompressibility condition,

[TABLE]

and the Stokes equation

[TABLE]

where $\bf{v}$ represents the velocity, $p$ the pressure, $\bf{g}$ the gravitational acceleration, $\mu$ the dynamic viscosity of the fluid, and $\rho$ the density of the fluid. We assume that the pure liquid phase abides by a Newtonian rheology. All non-Newtonian effects we detect are hence due to the interaction between solid and liquid phases.

We represent the particles as rigid bodies fully immersed in the fluid. In the solid domain, the governing equations are hence Newton’s equation of motion

[TABLE]

where ${\bf{X}}_{p}$ is the position of the particles’ center of mass, $M_{p}$ is the mass of the particle, ${\bf{V}}_{p}$ its velocity, ${\bf{\omega}}_{p}$ its angular velocity, ${\bf{I}}_{p}$ its angular moment of inertia tensor, and ${\bf{F}}_{p}$ and ${\bf{T}}_{p}$ represent the hydrodynamic force and torque exerted onto the particle by the surrounding fluid.

In this paper, we focus specifically on circular and rectangular particles, primarily because of the importance of these shapes for lava flows (see Figure 1). Simplified geometric forms as compared to general shapes are advantageous from a numerical point of view, because the moment of inertia tensor (see Figure 2) and the solid fraction in a given computational cell (see Figure 3) can be integrated analytically for idealized particle geometries. Some more details about these two analytical algorithms are introduced in following paragraphs.

For a circular particle rotating around its center of mass, the angular moment of inertia tensor is ${\bf{I}}_{p}=M_{p}{\bf{r}}_{pq}^{2}$, where we define $||{\bf{r}}_{pq}||$ to be the radius of the circular particle (see Figure 2A). For a rectangular particle (see Figure 2B), ${\bf{I}}_{p}$ is

[TABLE]

While assuming idealized crystal shapes is advantageous because it significantly improves numerical accuracy, our numerical method lends itself to straight-forward extension to general particle shapes through a high-order integration scheme.

In the interest of simplicity, we do not currently consider nucleation, growth, or dissolution of particles, implying that the velocity of point q on the edges of both circular and rectangular particles can be represented as

[TABLE]

This simplification implies the assumption that the time scale of crystal motion is faster than the time scale of crystal nucleation or growth, such that these process will only contribute a relatively negligible component to the overall dynamics of the flow. While clearly not generally true, we use this assumption as a first step towards a more complete understanding of the behavior of cyrstal-bearing magmatic flows.

In 2D, we define the solid fraction in a given computational cell as the ratio of the area occupied by particle to that of the whole cell. For the simplified geometric shapes we use, the solid fraction can be computed analytically by decomposing it into triangles, rectangles, and circular segments. Figure 3 shows two examples for computing the solid fraction in cell, $G$. The first step in the analytical computation of the solid fraction is to identify how many grid vertexes fall into the solid domain. In Figure 3(A), only one vertex, $D$, is inside the solid domain and its coordinates satisfy the following condition,

[TABLE]

where $r$ represents the radius of the circular disk, $x_{D}$ and $y_{D}$ are the x and y- coordinates of $D$, and $x_{O}$ and $y_{O}$ are the coordinates of the center of the disk, $O$. In the second step, we compute the coordinates of the intersection points between the interface and the cell faces, $B$ and $C$, through

[TABLE]

and,

[TABLE]

respectively. Finally, we calculate the solid fraction in cell $G$, $A_{p}$, as the sum of the triangle, $A_{BDC}$, and the circular sector cut off by the chord from B to C, A_{\stackrel{{\scriptstyle\mbox{\large\frown}}}{{BOC}}},

[TABLE]

where $A_{BDC}$($=\frac{1}{2}|BD||DC|$), $A_{BOC}$($=\frac{1}{2}|BO||CO|$) and

[TABLE]

with $\alpha=\text{cos}^{-1}(\frac{OB\cdot OC}{|OB||OC|})$.

Figure 3(B) illustrates an example for computing the solid fraction occupied by a rectangular particle. Similar to the computation for the circular disk, we first detect the number of vertexes inside the solid domain, $D$, by verifying

[TABLE]

where $l_{1-4}^{a}$, $l_{1-4}^{b}$ and $l_{1-4}^{c}$ are the coefficients of four line segments of the rectangle (marked in Figure 3(B)). We then compute $A_{p}$ as the area of the shaded triangle, $A_{BDC}=\frac{1}{2}|BD||DC|$, where the coordinate of the intersections $B$ and $C$ are,

[TABLE]

and,

[TABLE]

respectively. Needless to say, many more cases exist of how the grid might dissect the particles, but all follow the basic strategy of decomposing the solid fraction into simple geometric shapes for which the area can be computed analytically. An important advantage of this strategy is that the accuracy of the integration to obtain solid fraction does not depend on the orientation of the particle with respect to the grid, which can lead to numerical artifacts.

3 Numerical Methodology

Due to the small size of the crystals (see Figure 1) and the high viscosity of the fluid, the Reynolds number is very small in many magmatic flows. As a consequence, the Navier-Stokes equation can often be reduced to the Stokes equation. While nominally easier because of the absence of the nonlinear term and the parabolic time dependence, the Stokes equation is a purely elliptic boundary-value problem. Approximating it numerically depends sensitively on the accuracy with which the numerous interior solid-fluid interfaces are represented. Here, we derive our fluid solver in the Stokes limit. At finite Reynolds number, we rely on the Navier-Stokes we have developed previously [29].

The idea behind our methodology pioneered by Glowinski et al., [11] is to solve the Stokes equation in the entire computational domain at the beginning of each time step. In the first step, both phases are hence treated as a fluid. In the second step, we use the distributed Lagrange multipliers to prevent deformation in the solid domain by projecting fluid onto solid motion, while keeping momentum conserved. In the third step, we enforce a no-slip boundary condition in the immediate vicinity of the solid-fluid interface, $\Gamma$, through our immersed boundary method. We hence generalize the one-phase momentum conservation (Eq. 2) to the following two-phase equation,

[TABLE]

where, ${\bf{f}}_{DLM}$ is the body force solved by the distributed Lagrange multipliers, ${\bf{f}}_{IB}$ is the body force solved by the immersed boundary method. In the following subsections, we discuss the three main components of our numerical approach, namely 1) the projection-iterative solver of the Stokes equation on the Cartesian grid, 2) the distributed Lagrange multipliers; and 3) the immersed boundary method.

3.1 Stokes solver

Our solver for the incompressible Stokes equation is based on the projection method [57] with approximate factorization [58] and fractional time stepping [59] on a staggered grid. Using this solver, we actually replace the governing equation (Eqs. 1-2) by the following equations,

[TABLE]

[TABLE]

and solve them iteratively until the residual term on the left hand side falls below a predefined tolerance level, $\frac{\partial{\bf{v}}}{\partial t}<$ TOLS.

At each iteration n, we apply the projection method which treats the pressure term in Eq. 18 as a projection operator to project the initial velocity guess onto a divergence-free field (Eq. 17). In the framework of the classical three projection steps, n, n+1/2 and n+1 [57], we first obtained an provisional velocity field, $\hat{\bf{v}}$ at n+1/2,

[TABLE]

[TABLE]

where ${\bf{f}}_{IB}^{n}$ and ${\bf{b}}^{n}$ are the boundary condition determined by the immersed boundary method. We treat the viscous term in Eq. 19 as an implicit term, which allows us to eliminate the time-step restriction associated with the viscous term. Hence, we rewrite Eq. 19 in two dimensions as

[TABLE]

Replacing the large sparse coefficient matrix on the left hand side of Eq. 21 with tridiagonal matrices, we solve the following equation,

[TABLE]

instead of Eq. 21, associated with a second-order-accurate boundary condition [59],

[TABLE]

Kim and Moin [59] argued that, as an $O(\Delta t^{3})$ approximation of Eq. 21, Eq. 22 results in a huge reduction in computing cost and memory.

The provisional velocity field $\hat{\bf{v}}$ contributes to the final velocity field ${\bf{v}}^{n+1}$ and the gradient of the pressure later on,

[TABLE]

and

[TABLE]

[TABLE]

where ${\bf{n}}$ represents the normal of the solid-fluid interface.

We then solve the following Poisson equation

[TABLE]

implied by Eqs. 24-25 through an asymmetric multifrontal method [60]. Introducing Eqs. 20 and 26 into Eq. 24, we obtain the following boundary condition for the pressure on a Dirichlet boundary,

[TABLE]

In the present work, we do not impose pressure boundary conditions on $\Gamma$ explicitly. As illustrated later, the momentum equation normal to $\Gamma$ reduces to d$p$/d$n$ = 0 on the boundary points, because we linearize all velocity components in the vicinity of $\Gamma$ in our immersed boundary method.

3.2 Distributed Lagrange multipliers

During the first iteration of the fluid solver, the additional body forces, ${\bf{f}}_{DLM}$ and ${\bf{f}}_{IB}$ are unknown. We first solve the variable-density Navier-Stokes equation without these terms. To ensure that the deformation tensor in the solid domain is zero,

[TABLE]

we use distributed Lagrange multipliers to project the initial fluid velocity field onto rigid-body motion [11]. To estimate the linear and angular velocity of each particle, we take advantage of momentum conservation and integrate the fluid velocity inside the solid boundary to obtain an equivalent linear rigid-body velocity,

[TABLE]

and angular velocity,

[TABLE]

where ${\bf{v}}^{*}$ is the fluid velocity field solved from Stokes equation, which is represented by ${\bf{v}}^{n+1}$ in Section 3.1, $\Phi_{i,j}=\frac{A_{p}}{\Delta x\Delta y}$ represents the solid volume fraction in cell $(i,j)$, and $A_{p}$ is the area inside the solid domain for a given computational cell, which is determined analytically by the scheme introduced in Section 2.

3.3 Immersed boundary method

One drawback of relying on the Navier-Stokes equation in both the fluid and solid domain is that the interface between the domains abides by a free-slip condition at least initially. Since the distributed Langrange multiplier are only defined in the solid domain, they may not correct the interface itself leading to an inaccurate representation of the no-slip condition on the solid-fluid boundary particularly in Stokes flow [e.g., 29]. To improve the representation of the no-slip condition, we use an immersed boundary method to derive the hydrodynamic forces that result from the no-slip condition in the immediate vicinity of the solid interface.

Figure 4 illustrates our approach for the velocity in the solid domain and the velocity on the fluid in the immediate vicinity of the interface. Here, we introduce the x-component of this approach in detail, the y-component is similar. Since our solver is based on a staggered grid, $u$ is stored at the cell faces, such as $(i_{1}+\frac{1}{2},j_{1})$ and $(i_{3}+\frac{1}{2},j_{3})$ in Figure 4. In the solid domain, we define the rigid body motion (see Eq. 7) to obtain a corrected velocity,

[TABLE]

where $p_{y}$ represents the y-coordinate of particle center, thus $(y(j_{1})-p_{y})$ represents the y-component of ${\bf{r}}_{1}$.

In the fluid domain, e.g., at node c $(i_{3}+\frac{1}{2},j_{3})$, we employ a direct forcing method to substitute $u^{*}_{i_{3}+\frac{1}{2},j_{3}}$, the fluid velocity solved from Stokes equation, with the weighted average velocity,

[TABLE]

where $\phi$ is an interpolating factor, $u_{q}$ is the solid velocity at point $q$ that is just on the solid-fluid interface and $u_{f}$ represents the fluid velocity without the influence of the solid-fluid interface. As demonstrated later, a linear interpolation scheme in the normal direction with respect to the interface as proposed by [61] enables a 2nd-order accurate representation of the boundary condition.

Other weighting approaches include Luo et al. [13], who proposed a nonlinear-weighted average scheme depending on the particle Reynolds number Re${}_{p}=\frac{\rho_{f}|u_{q}-u_{f}|D}{\mu_{f}}$. In the framework of this scheme, $\phi$ is defined as follow,

[TABLE]

where $|X|=\frac{|h|}{D}$ is the relative distance, $D$ is the crystal diameter, and $h$ is the distance from the node to the solid-fluid interface. In the Stokes regime, $\text{Re}_{p}\approx 0$, implying $\phi\approx 1$ and the solid velocity is imposed to the fluid in the immediate vicinity of the solid-fluid interface, which is not realistic. Another nonlinear strategy, which sets $u_{f}$ to $u^{*}_{i_{3}+\frac{1}{2},j_{3}}$ and defines $\phi$ as the solid fraction ($\Phi$) in a computational cell, was also widely used in previous studies, but is limited to a 1st order accurate in space [31].

In terms of linear weightings, Qin and Suckale [29] compute $u^{n+1}_{i_{3}+\frac{1}{2},j_{3}}$ through a 1D linear interpolation, defining $\phi=h_{y}/(h_{y}+\Delta y)$ and $u_{f}=u^{*}_{i_{3}+\frac{1}{2},j_{3}+\frac{1}{2}}$, where $h_{y}$ is the distance from the node to the interface in y-coordinate. While this scheme provides 2nd-order accuracy for spherical shapes, it is not suitable for complex shapes. For example, at the corners of the rectangular crystal considered here, it is ambiguous about which unique direction is that the interpolation can be performed over. The scheme proposed by Balaras [61] eliminates this ambiguity by linearly interpolating along the normal direction with respect to the interface. Figure 4 illustrates this scheme, by which $u_{f}$ is defined as the velocity at point $o$, $u_{o}$, and $\phi$ is defined as $\frac{|cq|}{|oq|}$. Note that the point $o$ is not located on a grid node, therefore it is needed to interpolate $u_{o}$ from the surrounding points by

[TABLE]

where $\alpha_{i}$ are the coefficients of a standard bilinear interpolation involving points 1-4. We hence follow the scheme by Balaras [61] here.

4 Collision modeling

To model the collision and sustained contact between particles of different shapes, we implement a repulsive-potential model. More specifically, we build on the formulation by Ardekani and Rangel [41] by updating the particle velocity implicitly to avoid numerical instability and an overly restrictive Courant-Friedrichs-Lewy condition. One consequence is that our particles never come directly in contact, which is consistent with field observations of lava flows showing a small but finite separation between crystals (see Figure 1).

4.1 Detecting a Collision Between Two Particles

We define that a collision occurs when the closest distance between two particles, $d$, is less than or equal to a threshold distance $h$, where $h$ is selected to be twice of the grid spacing, $h=2\Delta x$, to ensure that the particles never overlap. Depending on the shapes of the colliding particles, we apply three different strategies to compute the closest distance $d$ and to identify the collision point (see Figure 5).

When two circular particles collide, the collision point always lies on the line connecting the two centers (see Figure 5A). We compute the closest distance $d$ by subtracting the sum of two radius from the distance between two centers, such that $d=||{\bf{x}}_{1}-{\bf{x}}_{0}||-(||{\bf{r}}_{\alpha}||+||{\bf{r}}_{\beta}||)$. The collision between two rectangles is more complex. We determine the closest distance $d$ by considering each particle as a collection of four points, the vertices of the rectangle, and four line segments, the edges of the rectangle. First, we compute the distance between each line segment of particle $\alpha$ and each vertex of particle $\beta$. The shortest distance over all of these computations can be denoted as $d_{\alpha\beta}$. Then, we compute the distance between each vertex of particle $\alpha$ and each line segment of particle $\beta$, and find the shortest distance over all of these computations, denoted as $d_{\beta\alpha}$. In the end, the distance $d$ is assigned to the shorter one of $d_{\alpha\beta}$ and $d_{\beta\alpha}$.

The distance between a point ${\bf{x}}_{0}$ and a line segment defined by points ${\bf{x}}_{1}$ and ${\bf{x}}_{2}$, is computed based on whether the collision involves two vertices, a vertex and a edge (see Figure 5B) or two edges (see Figure 5C).

Collision between two vertices: This is a special case of Figure 5B. If $(\rm{min}(x_{1},x_{2})>x_{0}\ \rm{or}\ \rm{max}(x_{1},x_{2})<x_{0})$ and $(\rm{min}(y_{1},y_{2})>y_{0}\ \rm{or}\ \rm{max}(y_{1},y_{2})<y_{0})$ where $x_{0,1,2}$ and $y_{0,1,2}$ are the components of ${\bf{x}}_{0,1,2}$, the distance is $\rm{min}(||{\bf{x}}_{1}-{\bf{x}}_{0}||,||{\bf{x}}_{2}-{\bf{x}}_{0}||)$. 2. 2.

Collision between a vertex and a edge: If (1) is not the case, then let

[TABLE]

where ${\bf{t}}$ is a normalized vector pointing from ${\bf{x}}_{1}$ to ${\bf{x}}_{2}$, ${\bf{\hat{n}}}$ perpendicular to ${\bf{t}}$ is a vector whose length is the distance represented by $d$. 3. 3.

Collision between two edges: The colliding particles are approximately parallel, as shown in Figure 5C. We keep using above strategy but replace ${\bf{x}}_{0}$ by $\bar{{\bf{x}}}_{0}$ which the center point of ${\bf{x}}_{0}$ and $s$.

Finally, it is possible for a circular and a rectangular particle to collide, as illustrated in Figure 5D. This case arises in mixed shape simulations of when a circular particle collides with one of the domain walls. In this case, the line connecting the point of collision and the center of circular particle is perpendicular to the rectangular particle’s closest edge to the center of the circular particle. The closest distance $d$ is then computed in a similar way to the one for two rectangles but using the radius of the circular object, $||{\bf{r}}_{\alpha}||$.

If $(\rm{min}(x_{1},x_{2})>x_{0}\ \rm{or}\ \rm{max}(x_{1},x_{2})<x_{0})$ or $(\rm{min}(y_{1},y_{2})>y_{0}\ \rm{or}\ \rm{max}(y_{1},y_{2})<y_{0})$, then the distance is $\rm{min}(||{\bf{x}}_{1}-{\bf{x}}_{0}||,||{\bf{x}}_{2}-{\bf{x}}_{0}||)-||{\bf{r}}_{\alpha}||$. 2. 2.

If (1) is not the case, then let

[TABLE]

4.2 Computing particle velocities after collision

Suppose rectangles $\alpha$ and $\beta$ are colliding at points $q$ and ${\bf{x}}_{0}$, as illustrated in Figure 6. The initial linear velocities at each particles’ center of mass are denoted with a superscript $I$ and subscript letter specifying the particle, such as ${\bf{V}}_{\alpha}^{I}$. The final central linear velocities, computed through our collision model, are denoted with a superscript $F$ and subscript letter specifying the particle, such as ${\bf{V}}_{\alpha}^{F}$. At the collision points, the velocities of each particle are

[TABLE]

Since there always is a thin liquid film between the colliding particles, the particles are free to move along the linear space spanned by the tangential vector, $\bf{t}$, without being affected by friction,

[TABLE]

The velocities parallel to ${\bf{n}}$ are constrained by the conservation of linear momentum, such that

[TABLE]

where $c_{R}$ is the coefficient of restitution. For a purely inelastic collision, $c_{R}=0$, and for a purely elastic collision $c_{R}=1$. For all of the simulations in this paper, $c_{R}=0$ since the viscosities are too high and the particles to small to be able to rebound. We compute the impact Stokes numbers for all validation cases involving multiple particles later to quantify this point better.

Since the rectangles are rigid bodies, the change in the central velocity that occurs from the collision happens instantaneously and can be represented with a pair of impulses, $f{\bf{n}}$ and $f(-{\bf{n}})$,

[TABLE]

Combining Eqs. 36 and 39 then giving them into Eq. 38 produces an equation that can be solved for $f$

[TABLE]

In the end, the corrected velocity of the colliding particles can be computed by giving $f$ and the constraint in the tangential direction, Eq.37, into Eq.39.

4.3 Multiple-particle collisions

As illustrated in Figure 1, lava flows are characterized by long-lived clustering between multiple particles. To capture multiple-particle collisions, we generalize our collision model to capture force chain between systems of $N$ colliding particles, as shown in Figure 7. In each cluster, one particle might collide with several other particles and hence experience a repulsive forces at the momentum of collision from multiple sources (e.g., particle $\alpha$ suffers the force from both $\beta$ and $\gamma$ in Figure 7). Assuming the collision between the system of particles creates $H$ repulsive forces, we define the final linear and angular velocities at the center of a single particle through a generalized form of Eq. 39,

[TABLE]

where $\kappa=1$ when particle $\xi$ collides with particle $\eta$ while $\kappa=0$ when $\xi$ and $\eta$ do not collide, $f_{\epsilon}$ is the potential repulse force affecting between $\xi$ and $\eta$, ${\bf{n}}_{\xi\eta}$ is the normal vector from $\xi$ to $\eta$, ${\bf{r}}_{\xi q_{\xi\eta}}$ is the vector from the center of $\xi$ to the collision point, $q_{\xi\eta}$. Note that $q_{\xi\eta}\neq q_{\eta\xi}$ because of the thin separating film introduced in above sections, like the black points versus the red points shown in Figure 7.

Applying above definition of the final velocities of each single particle (Eq. 41) into the non-overlapping constraint (similar to Eq. 38) for each colliding pairs, we formulate a $H\times H$ linear system to compute the net forces as follow,

[TABLE]

which can be seen as an extension of Eq. 40. In above linear system, $\bf{A}$ represents the coefficient matrix, where the rows corresponds with the particles involved in each collision pair between particles $k$ and $l$, and the columns corresponds with the repulse forces affecting on the particles $i$ and $j$, so

[TABLE]

where $i<j$, $k<l$. Each element in $\mathbf{A}$ is given by

[TABLE]

where

[TABLE]

We solve Eq. 42 by a sparsity exploiting linear programming package, UMFPACK (Unsymmetric MultiFrontal method and direct sparse LU factorization).

5 Verification

We are not aware of any analytical solutions of flow around rectangular bodies that we can use for verification purposes. We hence verify the accuracy of our coupled fluid-solid solver for two circular geometries. In the first case, we verify our numerical technique by reproducing the case of Stokes flow around a fixed cylinder in the vicinity of a moving plate. In the second case, we compute the velocity profile in a rotational viscometer to test whether our model correctly enforces the no-slip boundary condition on the moving solid-fluid interfaces.

5.1 Wannier flow

Stokes flow around a circular cylinder in the vicinity of a moving wall is named Wannier flow honoring an analytical solution by Wannier [51]. Using the analytical solution, we can quantify the numerical error of our coupled solid-fluid solver. Figure 8(A) shows the flow configuration and the computed streamlines on the grid of 64$\times$64. The Reynolds number for this case is $\sim$10*-3*. In the computation, we assign a moving wall boundary condition to the bottom boundary and the analytical solution on all other boundaries.

To test whether our simulation method correctly enforces the no-slip condition, we plot evolution of the maximum normal derivative of the pressure on the solid-liquid interface in the infinity norm (see Figure 8B). The normal derivative of the pressure reduces to 1% of the initial value during the second iteration and decreases to less than 10*-4* after a few iterations. We can hence capture the no-slip boundary condition iteratively, although we did not explicitly impose a pressure boundary condition.

Figure 8(C) shows the $L_{2}$ norm of the error between the numerical and the analytical solution of the velocity components in x direction. We implement and compare several immersed boundary methods here to compare their convergence properties. We find that the linear-weighting scheme by Balaras [61] and our implementation afford second-order convergence rates while the nonlinear-weighting schemes [31, 13] are limited to approximately first order.

5.2 Rotational viscometer

Rotational viscometers are routinely used in laboratory studies to quantify the rheological properties of solid-bearing suspensions [e.g., 53, 62]. Similar to Wannier flow, the velocity profile is sensitively dependent on enforcing the no-slip condition correctly. We hence use this case to demonstrate the value of the immersed boundary method. Figure 9(A) shows our 2D representation of a rotational viscometer that consists of a fixed outer cylinder (gray circle) and an inner cylinder (black line) rotating at constant angular velocity, $\omega$. We set the radii of the outer and inner cylinder to be same as the rotational viscometer used in [53], where $R_{1}$ = 0.003 m and $R_{2}$ = 0.015 m. For Stokes flow of a Newtonian fluid, the velocity profile, $V_{\theta}$, can be solved for analytically as,

[TABLE]

where $r$ represents the radial coordinate.

In the numerical simulation, we treat the cross sections of both cylinders as circular disks and enforce a finite angular velocity, $\omega=1$, on the inner cylinder and a zero angular velocity on the outer cylinder. The velocity of the center is zero for both disks. In the gap between two disks, we set the density and the viscosity of the fluid to, $\rho_{f}=1000$ kg/m3 and $\mu_{f}=1000$ Pa s, respectively. Because of the small length scale and the high viscosity, the Reynolds number is approximately $10^{-5}$. It is hence reasonable to consider the fluid as a Stokes fluid.

Figure 9 (B) shows the difference between the numerical radial velocity averaged over all $\theta$ on a grid of 128$\times$128 and the theoretical solution given by Eq. 46. We also plot the numerical result based only on distributed Lagrange Multipliers without enforcing the immersed boundary method and observe a significant drop in numerical accuracy. The drop of the numerical error quantified through the $L_{\infty}$ shows that we obtain second order convergence (see Figure 9 C).

6 Validation

Over the last couple of years, an increasing number of experimental studies have considered rectangular shapes that we can use for validation purposes. Here, we investigate flow past a rectangular cylinder at intermediate Reynolds number [54], joint settling of rectangular particles in a stagnant fluid [56] and jamming of spherical particles in a shear suspension [52] as test cases.

6.1 Flow past a square cylinder

To assess that our numerical scheme correctly resolves the impact of geometry on flow in the immediate vicinity of the boundary, we simulate the two-dimensional flow over a fixed square cylinder at $\mathrm{Re}=410$ and compare our numerical results to the analogue experiments by Dutta et al., [54]. Similar to the experimental setup, we fix a square cylinder of 3.4 mm edge (labeled ”D”) on the centerline of the computational domain whose width, W, is 4.8 cm. The boundary conditions imposed on the flow field are no-slip, ${\bf{v}}=0$, along the sidewalls. The inflow boundary is constant velocity, $u_{\text{inf}}$, and the outflow boundary, defined on the opposite side, is solved by a convective flow condition [63],

[TABLE]

where $\bf{n}$ is the normal vector on the boundary.

Figure 10 compares the contours of the numerical streamfuction to the experimental streak lines for four different cylinder orientation ($\theta=0,22.5,30,45^{\circ}$), and aspect ratios (16 and 28 for experiments and $\infty$ for the numerical simulation), where the aspect ratio is defined as the ratio of the height of the cylinder to the edge of its cross-section. The size of the corresponding computational domain is W$\times$2W, and the simulation is based on a resolution of $200\times 400$. In both the numerical simulation and the experimental results, the flow around the cylinder separates from the surface of the square and forms two recirculation vortex in the near wake. The unequal vortices presented in panels b, c, f, g, i and j demonstrate that the flow is asymmetric when the cylinder orientation is 22.5 or $30^{\circ}$. Comparing the width and length of the vortex between the cases with AR = 16 and AR = 28, we find that the vortices shrink at higher aspect ratio for all cylinder orientations. Moreover, the distance between the vortex and the square decreases at higher aspect ratio. These differences between higher and lower aspect ratio are likely a consequence of the three dimensionality of the flow field [54]. The 2D numerical simulation, which by its lower dimensionality implicitly assumes that the aspect ratio is infinite large, shows similar trends. Comparing panels i-l to e-h, we observe a smaller vortex that is much closer to the square.

Figure 11 compares the centerline recovery of the streamwise velocity and decay of transverse velocity for experiments and numerical simulation with various cylinder orientations. We plot the nondimensional streamwise velocity, $u_{\text{CL}}/u_{\text{inf}}$ along the x axis at the centerline and the nondimensional transverse velocity, $v/u_{\text{inf}}$, along the x axis at a offset location from the centerline (y=1). In agreement with the experiments, the time-average streamwise velocity in our simulations is zero on the square and is negative in the recirculation zone. Later, it increases and reaches an asymptotic value, which is in the range of 0.6-0.65 for all aspect ratios and all cylinder orientations. We obtain similar agreement between experimental and numerical results for the transverse velocity that switches from the maximum positive value to a negative value in the recirculation zone and later decays to approximately zero. Comparing the velocity profiles between the cases with different aspect ratios, we observe that both centerline recovery and the decay of transverse velocity are faster at the higher aspect ratio. The asymptotic limit of u${}_{\text{CL}}$ and the zero value of v are reached at around x = 15, 7, and 5, respectively for AR = 16, 28 and $\infty$. This trend is in accordance with the larger recirculation vortex at the lower aspect ratio, as shown in Figure 10.

Figure 12 shows the contours (represented by capital letters) and the profiles (represented by lower-case letters) of the time-average nondimensional vorticity for three aspect ratios and four cylinder orientations. In the numerical simulation, we compute the 2D nondimensional vorticity through direct differentiation of the velocity field,

[TABLE]

where $u$ and $v$ represent the streamwise velocity and the traverse velocity, respectively. In agreement with the experimental record, the maximum value of the numerical vorticity is observed at the cylinder corner (see pictures I-L). A comparison between the numerical and experimental vorticity profiles demonstrates that the vorticity magnitudes in the immediate vicinity of the cylinder are quite similar for different aspect ratios. Consistent with the observation on the recirculation vortex in Figure 10, the smaller spreading of the vortices occurs at higher aspect ratios, and the vortices are asymmetric when the cylinder orientation is 22.5 or $30^{\circ}$.

Figure 13A compares the numerical Strouhal number to the experimental values. We find good agreement between the numerical and experimental Strouhal number at the highest aspect ratio available experimentally, $\mathrm{AR}=60$. The numerical simulation correctly reproduces the qualitative trend, observed in experiments for all aspect ratios, that the Strouhal number increases from the $0^{\circ}$ cylinder orientation angle to $22.5^{\circ}$ but drops subsequently with a further increasing angle. The maximum Strouhal number arises around $22.5^{\circ}$. Dutta et al., [54] explained the angular variations in the Strouhal number through an angular dependence in the projected dimension of the cylinder with respect to the incoming flow.

Figure 13B shows the variation of the time-average drag coefficient with respect to the cylinder orientations. In 2D numerical simulations, the hydrodynamic drag force, $F_{D}$, acting on the cylinder cross-section in streamwise direction, is given by

[TABLE]

where $S$ represents the circumference of the cylinder cross-section. Hence, we compute the drag coefficients as,

[TABLE]

In agreement with the experimental results, our simulation shows that the minimum drag coefficient occurs at $22.5^{\circ}$ angle. Contrary to the Strouhal number, the drag coefficient decreases from $0^{\circ}$ to $22.5^{\circ}$ and increases again for larger angles.

6.2 Jamming in shear suspensions

Earlier studies, including numerical [52] and experimental studies [55], observed that circular particles in a sheared suspension are prone to clustering, which alters the rheology of the suspension significantly. Assuming spherical particles, Brady and Bossis [52] simulated suspension behavior in a Couette device. They observed that particle clusters at intermediate particle fraction begin to span the gap between the two plates, resulting in a boundary-dominated flow where the suspension moves approximately as a plug with an average speed of approximately half of the imposed shear speed. Karnis et al. [55] experimentally observed similar behavior in both shear flow in a cylindrical Couette device and pressure-driven flow in tubes.

In our numerical model, we use two elongated, rectangular particles to represents the parallel planes of a Couette cell. We fix the upper particle, $\alpha$, and move the lower particle, $\beta$, with a constant speed $U$. We set the width of the gap between plates, $L$, to $L$ = 0.012 m. Since the motivation for our numerical method is to simulate crystal-bearing lava flows, our parameters are inspired by a magmatic context with a fluid viscosity of $\mu_{f}=1000$ Pa s and a density of $\rho_{f}=3000$ kg/m3. With these material properties, the Reynolds number is approximately $\sim 10^{-5}$. We then introduce 120 neutrally buoyant, circular particles with radii of $r=1/20L$.

Figure 14 (A)-(C) show three snapshots in time of the instantaneous crystal configurations at a solid area fraction of $\psi\approx 0.20$. The color scale represents the velocity in the flow direction. As a comparison, three crystal configuration are reproduced from Brady [52]. Initially, the particles are uniformly distributed in the domain. Soon, they start forming clusters spanning from the lower to the upper plate (Figure 14A). The number of particles in the cluster increases with time (Figure 14B) until most of the particles are clustered (Figure 14C). At this point, the particles translate more or less as a single entity or plug. We plot the average nondimensional velocity, $<\bar{u}>$, of the suspension in the flow direction versus the transverse coordinate $y$ in Figure 14 (D). The result by Brady [52] is shown in (E).

Following Brady [52], we define $<\bar{u}>$ as,

[TABLE]

where $\bar{u}$ represents the average dimensional velocity in the flow direction, $U$ is the speed of the lower plane and 0 is the speed of the upper plane. For comparison, we plot $<\bar{u}>$ for the pure fluid flow as a dashed line in (D). We obtain good agreement with [52] and [55] in reproducing that the suspension, including the crystal and the fluid, moves at roughly half the imposed boundary speed in the plug-flow limit. As a consequence, two regions of rapid shear form in the immediate vicinity of the bounding planes.

6.3 Settling and aggregation of rectangular particles

We validate the ability of our numerical method to capture the dynamic interactions between moving rectangular particles by reproducing the laboratory experiments by Schwindinger [56]. In this suite of experiments, prism-shaped particles are dropped into cold Karo syrup to observe settling behavior. Our simulations use the experimental parameters with crystal sizes of $l_{1}\times l_{2}=0.5\times 0.25$ cm, crystal densities of $\rho_{p}$=1.372-1.423 ($\pm$0.004) g/cm3, a constant fluid density of $\rho_{f}=1.37$ g/cm3, and a viscosity of $\mu_{f}=35$ Pa s. The simulations are carried out in a 2D rectangular computational domain which has the same dimeter as the cylinder used in the experiments. We apply a no-slip boundary condition to the wall of the computational domain. In the simulation and the analogue experiments, the particle Reynolds number is $\sim$ 10*-6*.

Figure 15 1-6 shows an experiment recording the settling-induced interactions between nine particles. Initially, the particles are arranged in three rows (Figure 15 1a). Shortly after the onset of motion, they aggregate into clusters (Figure 15 1b). One of the crystal clusters separates from the larger aggregate (Figure 15 1c, d) and settles more rapidly. Later, other clusters of two to three particles separate from the other particles and settle (Figure 15 1e,f). Different from spheres, which tend to form polygonal arrangements [56], the non-spherical particles are drawn into the center of the domain and create pairs or triplets. Our numerical simulations show similar behavior as the experiments.

7 Particle rotation may alter suspension behavior

The behavior observed in the Schwindinger experiments [56] is reminiscent of drafting, kissing, and tumbling [64], where the low-pressure wake behind a leading particle drafts in a trailing particle of comparable density until the two came in contact and then tumble apart. At very low Reynolds ($\mathrm{Re}\ll 1$), however, this mechanism no longer applies, because of the absence of low-pressure wakes. Nonetheless, Schwindinger observed that a tailing, lighter particle can catch up with a leading, denser particle even at approximately zero Reynolds number. Figure 16 1-10 shows an experiment where three particles is dropped into the Karo syrup. In this particular experimental setup, the tailing and leading particles ($\rho$ = 1.423 g/cm3) are slightly denser than the middle particle ($\rho$ = 1.372 g/cm3). Schwindinger reports similar behavior even for equal densities.

In agreement with the analogue experiment, our simulations in Figure 16 a-j show that the middle particle accelerates through the meniscus of seemingly disturbed fluid between the tailing and leading particles and catches up with the leading one. The reason for this behavior, however, is not related to spatial pressure variations but to the long-range interactions of the particles. As shown in Figure 16 a-j, the middle particle experiences a mean flow field that is already downward oriented. Its motion is a superposition of mean-field transport and individual settling and hence faster than buoyancy-driven settling alone.

An indirect indication of this mechanism is the differing rotational behavior observed by Schwindinger [56]. Figure 17 A shows the observed ratio of settling speeds for the tailing, $V_{tr}$, and the leading, $V_{ld}$ particle. In all of the shown experiments, the tailing particles approach the leading particles ($V_{tr}/V_{ld}>1$), but at different rates. High approach rates (i.e., experiments 21 and 14) are associated with more pronounced rotation of the particles as visually evident in the plotted interface positions and quantified more systematically in Figure 17 B. We compute the rotational energy by $E_{r}=1/2I_{p}\omega_{p}^{2}$ and compute the translational energy by $E_{T}=1/2M_{p}V_{p}^{2}$, where $I_{p}$ is particle’s angular moment of inertia tensor, $\omega_{p}$ is its angular velocity, $M_{p}$ is its particle mass and $V_{p}$ is its velocity. Rotations are indicative of unsteadiness in the mean flow since isolated non-spherical particles do not rotate when immersed in a steady-state flow at zero Reynolds number [25]. Rotations are hence suggestive of variability in the mean flow field, which in this particular case, leads to a relative speed-up of the tailing particle.

While rotations are a proxy for existing variability in the mean flow, they also generate variability in the mean flow. This effect is particularly pronounced in experiment 21 and simulation 1 from Figure 17, where both particles experience significant rotation due to their interactions, and most suppressed in experiment 16 and simulation 4, where the crystals are in close contact from the beginning and aligned in a way that makes them less susceptible to rotation. Contrary to the well-known drafting, kissing, and tumbling behavior at finite Reynolds number [64], this zero-Reynolds variant of drafting, kissing, and tumbling (e.g., experiments 21 and simulation 1) is hence highly dependent on crystal shape and the initial position of the crystals. Nonetheless, it is another instance of inertia-like behavior in solid-bearing suspensions at zero Reynolds number and highlights the nonlinearity and stochasticity introduced into these flows through long-range particle interactions analogous to the effect of turbulence at high Reynolds number [6, 7, 8].

One could argue that in the Schwindinger experiments [56], rotation is mostly a diagnostic of spatial variability in the mean flow. It is helpful for understanding the tendency of particles to aggregate and differences in settling speed, but it is not clear whether it is an essential component of the dynamics that a large-scale model of particle settling would need to capture. Potentially, the behavior of the clay particles in [56] could be approximated by considering an assemblage of spheres with variable size to represent clustering. At low solid fraction, that argument has some merit. It is not clear whether it generalizes to larger solid fractions, however, where jamming becomes important, requiring stable, long-term contact between particles to build force chains.

To test the importance of rotation in this regime, we revisit the simulation of a suspension sheared in a Couette device from Section 6.2 using square crystals. In the interest of easier comparability, we keep the overall solid fraction the same ($\psi=0.2$) and define our square particles such that a single square has approximately the same area as one of circular particles used previously. Figure 18 shows that the plug-like flow observed in Figure 14 for circular particles no longer occurs. In fact, the flow field in the fluid phases changes rather profoundly from the linear profile typical of Couette flow (Figure 14A) to a partially channelized flow field with a maximum flow speed in the center of the domain, not unlike Poiseuille flow (Figure 18C) if unsteady.

Figure 18D compares the profile of the velocity in the flow direction averaged over the transverse direction, $\bar{u}$, for circular and square particles. The velocity profile in the suspension of square crystals has two peaks at approximately $y=0.002$ m and $y=0.008$ m. The reason for this behavior is that the square particles do not distributed equally throughout the domain as highlighted in Figure 18E, where we plot average solid fraction averaged in the flow direction for the two suspensions. The circular crystal suspension reaches maximum solid fraction in the two rapid shear regions in the immediate vicinity of the bounding plates. In contrast, the square particles accumulate predominantly in the lower portion of the domain, thereby leaving space for almost solid-free and hence comparatively rapid flow in the upper third of the domain. In the example shown here, the flow speed reaches its peak value around $y=0.008$ m in the approximately particle-free channel.

In this case, the particle shape has altered flow behavior more fundamentally than in the Schwindinger experiments [56]. Aggregation within the solid phase has forced a shift from Couette-type flow to Poiseuille-type flow in the fluid phase. Our simulations suggest that the increased tendency of square particles to rotate is essential for creating this shift. Figure 18F shows the ratio of rotational to translational energy for square as compared to circular particles. While the ratio of kinetic energy stored in rotations is low as compared to translation, rotations have a strong impact on the dynamics, because they destabilize force chains. As a consequence, force chains are much more long-lived for circular as opposed to square particles. Instead of forming force chains throughout the domain, square particles aggregate and form clusters at the bottom of the domain, aided by their much higher random packing fraction (e.g., Figure 14C). The flow is diverted around the crystal-rich zone and channelizes, which in turn makes it more difficult for a force chain to build across it.

Evidently, this effect is dependent on the solid fraction of the suspension as illustrated in Figure 19. At the lowest solid fraction, $\psi=0.05$, the Couette flow is only slightly disturbed by the motion of the crystals (Figure 19 A), but starts to collapse to an unsteady Poiseuille-type profile at $\psi=0.1$ as cluster form at the two bounding plates. For $0.1\leq\psi\leq 0.2$, the Poiseuille flow profile occupies most of the computational domain, but it increasing localizes around $\psi=0.3$. At this stage, several large particle clusters form that tend to break up the flow in some locations. Nonetheless, a Poiseuille-type flow persists in spatially and temporally averaged sense even at $\psi=0.3$ (see Figure 19D). We compare the ratio of rotational to translational energy for the particles at different solid areal fractions to find that rotations remain pronounced until $\psi=0.3$ (Figure 19F). The rotational to translational energy ratio is lowest at $\psi$=0.05 because particles are, on average, separated the furthest in that case and interactions are less pronounced.

8 Conclusion

The goal of this study is to better understand how rotations of non-spherical particles alter suspension behavior. Since the stochasticity of solid-bearing suspensions emerges from the long-range interactions between particles, we develop a direct numerical technique for fully resolving solid-fluid coupling at the scale of individual interfaces. The broader motivation for this work is to better understand and model hazardous, crystal-bearing lava flows and we hence focus specifically on rectangular and square particles commonly found in field samples. Our numerical method relies on distributed Lagrange Multipliers to enforce rigid-body motion of the particles, coupled with an immersed interface method to correctly enforce the no-slip constraint on the solid interfaces. To prevent overlap between particles, we build a collision model that applies to both circular and rectangular shapes and correctly computes the translational and rotational momentum transfer during collision. We have verified and validated the efficiency and accuracy of our method at low and intermediate Reynolds number. Our simulations show that particle interactions tend to amplify rotational motion of square and rectangular crystals both at low and intermediate solid fraction. The consequences range from introducing unexpected dynamics such as drafting, kissing and tumbling behavior at zero Reynolds number to completely altering the flow field in the fluid. While the current model is 2D, these effects are likely to persist in 3D, since the rotational energy of non-spherical particles would tend to be greater and hence potentially even more important dynamically in higher dimensions.

Acknowledgement

This research was supported by the US Army Research Office through grant W911NF-18-1-0092. We would like to thanks Paul Wallace for encouraging us to look at crystal synneusis and Katherine V. Cashman for permission to use her thin-section images in Figure 1.

Author contributions

Z. Q. developed the code, performed the numerical simulations and created the figures. K. A. developed the collision model for rectangular particles. J. S. contributed to code development as well as study conception and analysis. All authors contributed to the text.

Reference