Extension of Moving Particle Simulation including rotational degrees of freedom for dilute fiber suspension

In the MPS model, the dynamics of fluid velocity obey the continuum Navier-Stokes equation. To incorporate the rotational degrees of freedom into the dynamics model, we employ the micropolar fluid model[18] in which the angular velocity field is incorporated. The conservation laws of linear and angular momentum are written as follows:

	$$\displaystyle\frac{D\bm{u}(\bm{r},t)}{Dt}=-\frac{1}{\rho}\bm{\nabla}P(\bm{r},t)+\nu\bm{\nabla}^{2}\bm{u}(\bm{r},t)+\nu_{r}\bm{\nabla}\times\bm{\Upsilon}(\bm{r},t)+\bm{f}(\bm{r},t),$$		(1)
	$$\displaystyle\mathcal{I}\frac{D\bm{\Omega}(\bm{r},t)}{Dt}=\bm{G}(\bm{r},t)-\nu_{r}\bm{\Upsilon}(\bm{r},t),$$		(2)

where $D/Dt$ is the time material derivative, $\bm{r}$ is the position vector, $t$ is time, $\bm{u}(\bm{r},t)$ is the fluid velocity, $\rho$ is the mass density, $P(\bm{r},t)$ is the pressure, $\nu$ is the kinematic viscosity coefficient, $\nu_{r}$ is the rotational kinematic coefficient, $\bm{\Omega}(\bm{r},t)$ is the angular velocity field, $\bm{\Upsilon}(\bm{r},t)=2\bm{\Omega}(\bm{r},t)-\bm{\nabla}\times\bm{u}(\bm{r},t)$, $\bm{f}(\bm{r},t)$ is the external volume force, $\mathcal{I}$ is the micro-inertia coefficient, and $\bm{G}(\bm{r},t)$ is the torque density due to the external field. According to the second law of thermodynamics, $\nu_{r}$ is a parameter properly chosen in the following range[19]:

$$0\leq\nu_{r}\leq\left(1+\frac{2}{d}\right)\nu.$$

(3)

Here, $d$ is the spatial dimension. For a normal fluid without micropolar degrees of freedom, $\bm{\Omega}$ is given as $\bm{\Omega}=\left(\bm{\nabla}\times\bm{u}\right)/2$ which guarantees that Eq. (1) reduces the standard Navier-Stokes equation[19]. In this work, we simply set $\bm{\Omega}=(\nabla\times\bm{u})/2$ for liquid region.

In this study, we employ the explicit MPS (EMPS) method[20, 21] to discretize Eqs. (1) and (2). The equations for the constituent particle $i$ are as follows:

	$$\displaystyle\frac{d\bm{u}_{i}(t)}{dt}=-\frac{1}{\rho_{i}}\left\llangle\bm{\nabla}P\right\rrangle_{i}\!(t)+\frac{1}{\mathrm{Re}}\left\llangle\bm{\nabla}^{2}\bm{u}\right\rrangle_{i}\!(t)+\frac{1}{\mathrm{Re}_{r}}\left\llangle\bm{\nabla}\times\bm{\Upsilon}\right\rrangle_{i}\!(t)+\bm{f}_{i}(t),$$		(4)
	$$\displaystyle\frac{d\bm{\Omega}_{i}(t)}{dt}=\alpha\bm{G}_{i}(t)-\frac{2\alpha}{\mathrm{Re}_{r}}\bm{\Upsilon}_{i}(t),$$		(5)
	$$\displaystyle\bm{\Upsilon}_{i}(t)=2\bm{\Omega}_{i}(t)-\left\llangle\bm{\nabla}\times\bm{u}\right\rrangle_{i}\!(t),$$		(6)

where $\llangle\rrangle$ indicates the quantity evaluated by the operator model in MPS at the position of particle $i$ mentioned in the next paragraph. The equations are non-dimensionalized using the following quantities: the fluid mass density $\rho_{0}$, the reference kinematic viscosity coefficient $\nu_{0}$, and the size of the fluid particle $l_{0}$. $\nu_{0}$ is a reference value, and $l_{0}$ can be interpreted as the characteristic length scale of the discretized system (which may be interpreted as the grid size in the finite difference scheme). These quantities define units of length, time, and energy, and the quantities discussed below are normalized according to these units. $\mathrm{Re}=\nu_{0}/\nu$ is the Reynolds number, $\mathrm{Re}_{r}=\nu_{0}/\nu_{r}$ is the rotational Reynolds number, and $\alpha$ is defined as $\alpha=l_{0}^{2}/\mathcal{I}$. As the case of the integration of the micropolar fluid model to the SPH model[19], translational and rotational velocities are mapped onto constituent (liquid and solid) particles. In our model, solid particles are micropolar fluid particles, and their motion follows Eqs. (4) and (5). The motion of liquid particles follows the standard Navier-Stokes equation plus the reaction force based on the third term on the right-hand side of Eq. (4) exerted by the surrounding solid particles.

To calculate the physical quantities and their differentials at the position of particle $i$, we need the weighting function. We employ the following weighting function:

$$w(r)=\begin{cases}l_{c}/r\,-1&(0<r<l_{c})\\ 0&(r\geq l_{c})\end{cases}.$$

(7)

Here, $l_{c}$ is the cutoff radius. The local density is evaluated by the local number density of the constituent particles defined as

$$n_{i}=\sum_{j\neq i}w\left(\left|\bm{r}_{j}-\bm{r}_{i}\right|\right).$$

(8)

The differential operators in Eqs. (4) and (5) are calculated by the following operator models:

	$\displaystyle\left\llangle\bm{\nabla}\psi\right\rrangle_{i}$	$\displaystyle=\frac{d}{n_{0}}\sum_{j\neq i}\left[\frac{\psi_{i}+\psi_{j}}{|\bm{r}_{j}-\bm{r}_{i}|^{2}}(\bm{r}_{j}-\bm{r}_{i})w\left(|\bm{r}_{j}-\bm{r}_{i}|\right)\right],$		(9)
	$\displaystyle\left\llangle\bm{\nabla}\times\bm{b}\right\rrangle_{i}$	$\displaystyle=\frac{d}{n_{0}}\sum_{j\neq i}\left[\frac{\left(\bm{b}_{j}-\bm{b}_{i}\right)\times\left(\bm{r}_{j}-\bm{r}_{i}\right)}{|\bm{r}_{j}-\bm{r}_{i}|^{2}}(\bm{r}_{j}-\bm{r}_{i})w\left(|\bm{r}_{j}-\bm{r}_{i}|\right)\right],$		(10)
	$\displaystyle\left\llangle\bm{\nabla}^{2}\bm{b}\right\rrangle_{i}$	$\displaystyle=\frac{2d}{\lambda n_{0}}\sum_{j\neq i}\left[(\bm{b}_{j}-\bm{b}_{i})w\left(|\bm{r}_{j}-\bm{r}_{i}|\right)\right],$		(11)
	$\displaystyle\lambda$	$\displaystyle=\frac{\sum_{j\neq i}{\left(\bm{r}_{j}-\bm{r}_{i}\right)}^{2}w\left(\left|\bm{r}_{j}-\bm{r}_{i}\right|\right)}{\sum_{j\neq i}w\left(\left|\bm{r}_{j}-\bm{r}_{i}\right|\right)}.$		(12)

Here, $\psi_{i}$ and $\bm{b}_{i}$ are scalar and vector variables on the particle $i$, $n_{0}$ is the initial particle number density, and $\lambda$ is the parameter defined by Eq. (12)[7]. Note that in Eq. (9) we use $\psi_{i}+\psi_{j}$ instead of $\psi_{j}-\psi_{i}$, as proposed by Oochi et al.[20], for better momentum conservation.

The fiber is modeled as an array of solid particles as shown in Fig. 1. The solid particles are connected with stretching, bending, and torsional potential energies, in a similar manner proposed by Yamamoto and Matsuoka for the other simulation scheme[22]. These potential forces should be a function of the bond vector of neighboring particles and the orientation of each solid particle[23]. To describe the orientation of the solid particles, we introduce two directors $\bm{s}_{i}$ and $\bm{t}_{i}$ on each solid particle $i$. $\bm{s}_{i}$ and $\bm{t}_{i}$ are unit vectors for which directions are parallel and perpendicular to the bond vector, as shown in Fig. 1 The time derivative of directors is related to the angular velocity as follows:

$$\frac{d\bm{s}_{i}(t)}{dt}=\left(\bm{1}-\bm{s}_{i}\bm{s}_{i}\right)\cdot\left(\bm{\Omega}_{i}\times\bm{s}_{i}\right),\quad\frac{d\bm{t}_{i}(t)}{dt}=\left(\bm{1}-\bm{t}_{i}\bm{t}_{i}\right)\cdot\left(\bm{\Omega}_{i}\times\bm{t}_{i}\right).$$

(13)

Here, $\bm{1}$ is the unit tensor. The projection tensors $(\bm{1}-\bm{s}_{i}\bm{s}_{i})$ and $(\bm{1}-\bm{t}_{i}\bm{t}_{i})$ maintain $\bm{s}_{i}\cdot\bm{t}_{i}=0$ within numerical errors.

The stretching potential $U_{s}$, bending potential $U_{b}$, and torsional potential $U_{t}$ are defined as

	$\displaystyle U_{s}\left(\left\{\bm{r}_{i}\right\}\right)$	$\displaystyle=\sum_{\left\langle i,j\right\rangle}\frac{k_{s}}{2}{\left(\left|\bm{r}_{j}-\bm{r}_{i}\right|-1\right)}^{2},$		(14)
	$\displaystyle U_{b}\left(\left\{\bm{r}_{i}\right\},\left\{\bm{s}_{i}\right\}\right)$	$\displaystyle=\sum_{\langle i,j\rangle}\left[\frac{k_{b}}{2}{\left(\bm{s}_{j}-\bm{s}_{i}\right)}^{2}-\frac{k_{r}}{2}\left\{{\left(\bm{s}_{i}\cdot\frac{\bm{r}_{j}-\bm{r}_{i}}{|\bm{r}_{j}-\bm{r}_{i}|}\right)}^{2}+{\left(\bm{s}_{j}\cdot\frac{\bm{r}_{i}-\bm{r}_{j}}{|\bm{r}_{i}-\bm{r}_{j}|}\right)}^{2}\right\}\right],$		(15)
	$\displaystyle U_{t}\left(\left\{\bm{t}_{i}\right\}\right)$	$\displaystyle=\sum_{\langle i,j\rangle}\frac{k_{t}}{2}{\left(\bm{t}_{j}-\bm{t}_{i}\right)}^{2}.$		(16)

Here, $k_{s},k_{b},k_{r},k_{t}$ are the spring constants and $\left\langle i,j\right\rangle$ represents a pair of two adjacent solid particles. The potential force $\bm{f}_{i}$ and torque $\bm{G}_{i}$ are calculated as

$$\bm{f}_{i}=-\frac{\partial\left(U_{s}+U_{b}\right)}{\partial\bm{r}_{i}},\quad\bm{G}_{i}=\bm{s}_{i}\times\left(-\frac{\partial U_{b}}{\partial\bm{s}_{i}}\right)+\bm{t}_{i}\times\left(-\frac{\partial U_{t}}{\partial\bm{t}_{i}}\right).$$

(17)

According to Eq. (14), if $k_{s}$ is large sufficiently, the fiber length $L$ corresponds to the number of solid particles in the fiber. Since the unit length of the system is the size of the fluid particle, the aspect ratio of the fiber $r_{p}$ corresponds to $L$.

In the EMPS method, the fractional step algorithm is applied for time integration as in the original semi-implicit scheme for MPS. Each integration step is divided into prediction and correction steps. In the prediction step, predicted velocity $\bm{u}_{i}^{\ast}$ is calculated by using terms other than the pressure gradient term in Eq. (4), and the angular velocity of the solid particles is also updated according to Eq. (5) as follows:

	$\displaystyle\bm{u}_{i}^{\ast}$	$\displaystyle=\bm{u}^{k}_{i}+\Delta t\left[\frac{1}{\mathrm{Re}}\left\llangle\bm{\nabla}^{2}\bm{u}\right\rrangle_{i}^{k}+\frac{1}{\mathrm{Re}}_{r}\left\llangle\bm{\nabla}\times\bm{\Upsilon}\right\rrangle_{i}^{k}+\bm{f}_{i}^{k}\right],$		(18)
	$\displaystyle\bm{\Omega}_{i}^{k+1}$	$\displaystyle=\bm{\Omega}_{i}^{k}+\Delta t\alpha\left[\bm{G}_{i}^{k}-\frac{2}{\mathrm{Re}_{r}}\bm{\Upsilon}_{i}^{k}\right],$
	$\displaystyle\bm{\Upsilon}_{i}^{k}$	$\displaystyle=2\bm{\Omega}_{i}^{k}-\left\llangle\bm{\nabla}\times\bm{u}\right\rrangle_{i}^{k}.$

Here, $\Delta t$ is the step size, and the upper indexes $k$ represent the step number: $\bm{b}_{i}^{k}=\bm{b}_{i}(t=k\Delta t)$. The predicted position $\bm{r}_{i}^{\ast}$ and directors are updated as

	$\displaystyle\bm{r}_{i}^{\ast}$	$\displaystyle=\bm{r}_{i}^{k}+\Delta t\bm{u}_{i}^{\ast},$		(19)
	$\displaystyle\bm{s}_{i}^{k+1}$	$\displaystyle=\bm{s}_{i}^{k}+\Delta t\left(\bm{1}-\bm{s}_{i}^{k}\bm{s}_{i}^{k}\right)\cdot\left(\bm{\Omega}_{i}^{k+1}\times\bm{s}_{i}^{k}\right),$
	$\displaystyle\bm{t}_{i}^{\ast}$	$\displaystyle=\bm{t}_{i}^{k}+\Delta t\left(\bm{1}-\bm{t}_{i}^{k}\bm{t}_{i}^{k}\right)\cdot\left(\bm{\Omega}_{i}^{k+1}\times\bm{t}_{i}^{k}\right),$

where $\bm{t}_{i}^{\ast}$ is the predicted torsional director. To maintain the relation $\bm{s}_{i}\cdot\bm{t}_{i}=0$, we adjust $\bm{t}$ as follows:

$$\bm{t}_{i}^{k+1}=(\bm{1}-\bm{s}_{i}^{k+1}\bm{s}_{i}^{k+1})\cdot\bm{t}_{i}^{\ast}.$$

(20)

In the correction step, the velocity and position are calculated as

	$\displaystyle\bm{u}_{i}^{k+1}$	$\displaystyle=\bm{u}_{i}^{\ast}-\frac{\Delta t}{\rho_{i}}\left\llangle\bm{\nabla}P\right\rrangle_{i}^{k+1},$		(21)
	$\displaystyle\bm{r}_{i}^{k+1}$	$\displaystyle=\bm{r}_{i}^{\ast}+\left(\bm{u}_{i}^{k+1}-\bm{u}_{i}^{\ast}\right)\Delta t.$

In the EMPS, the pressure is calculated by the following explicit form[20]:

$$P_{i}^{k+1}=\frac{\rho_{i}c_{s}}{n_{0}}\left(n_{i}^{\ast}-n_{0}\right).$$

(22)

Here, $c_{s}$ is the sound speed, and $n_{i}^{\ast}$ is the number density at $\bm{r}^{\ast}$. This $c_{s}$ is optimized concerning reasonable incompressibility and numerical stability.

We apply shear flows in the following boundary conditions. Hereafter, we refer to flow, shear gradient, and vorticity directions as $x$, $y$, and $z$ directions. We employ periodic boundary conditions for $x$ and $z$ directions, whereas we place solid walls at $y=0$ and $h$ perpendicular to the $y$ direction. These walls consist of three layers of liquid particles, which are fixed on a squared lattice. Following the earlier study[17], we move the walls toward the $x$ direction with the speed of $u_{\mathrm{wall}}=\pm\dot{\gamma}h/2$, where $\dot{\gamma}$ is the apparent shear rate. We have confirmed that the actual shear rate is equal to $\dot{\gamma}$ and uniform throughout the system within a numerical error, in simulations without fibers, as shown in Appendix A.

Simulations of a single fiber in a simple shear flow were carried out, and the rotational motion of the fiber was observed. To describe the fiber motion, we use the orientation angles $\phi$ and $\theta$ as shown in Fig. 2. The number of MPS particles was $N=64000$ in total including those for walls and the fiber. The simulation box dimension was $40\times 40\times 40$ in $x$-$y$-$z$ directions, respectively, and the distance between the walls was $35$. The kinematic viscosity coefficient $\nu$ and the strain rate $\dot{\gamma}$ were chosen so that the fiber-based Reynolds number was $\mathrm{Re}_{f}=L^{2}\dot{\gamma}/\nu=0.1$ to realize a viscous dominant condition. The sound speed of the fluid $c_{s}$ was set so that the Mach number became $\mathrm{Ma}=0.5h\dot{\gamma}/c_{s}=0.03$. The numerical step size $\Delta t$ was chosen to $0.01$ according to the Courant condition, the viscous constraint, and the relation to the spring constant. Other model parameters were set as $l_{c}=3.1$, $\nu_{r}=1.5\nu$, $\mathcal{I}=0.8$ unless otherwise noted. The mass density of the solid particles is the same as that of liquid particles. The aspect ratio of the fiber $r_{p}$ was varied in the range from 2 to 20. The spring constants were chosen at $k_{s}=1000$ and $k_{b}=k_{t}=k_{r}=200$. These values realized a rigid fiber, for which the effect of fiber deformation is negligible in the result as shown later. We performed the simulations with a house-made code.

In the initial condition, we placed the fiber at the center of the simulation box to overlap the center of mass of the fiber and the box. The initial fiber orientation angle to $x$-direction, $\phi_{0}$, was fixed at $\pi/2$, whereas the initial angle to $z$-direction, $\theta_{0}$, was chosen at $\pi/6$, $\pi/3$, or $\pi/2$. Surrounding liquid particles were randomly arranged by the particle packing algorithms proposed by Colagrossi et al. [24].