A.1. Oblique coordinate system

To impose simple shear flow on the system, a time-dependent oblique coordinate $\hat {\boldsymbol {r}}$ evolving with mean shear velocity $\boldsymbol {U}=\dot {\gamma }r^2\boldsymbol {e}_1$ is introduced as

(A1)\begin{gather} \hat{r}^1=r^1-\dot{\gamma}tr^2, \end{gather}

(A2)\begin{gather}\hat{r}^2=r^2, \end{gather}

(A3)\begin{gather}\hat{r}^3=r^3, \end{gather}

(A4)\begin{gather}\hat{t}=t, \end{gather}

where the quantities with a caret $(\,\hat {\cdot }\,)$ represent variables observed in the oblique coordinate system, and the upper indices 1, 2 and 3 represent the shear flow, velocity gradient and vorticity direction, respectively. By introducing an oblique coordinate system, advection by the mean flow, whose term explicitly depends on $r^2$, i.e. $(\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {\nabla })= \dot {\gamma }r^2\partial /\partial r^1$, is eliminated from the shear-enforced hydrodynamic equations. This enables the use of the periodic boundary conditions (Rogallo Reference Rogallo1981; Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011; Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016). From the coordinate transformation, the covariant and contravariant transformation matrices $[\boldsymbol {\varLambda }]_{\nu \mu }=\varLambda ^{\nu }{}_{\mu }=\partial r^{\nu }/\partial \hat {r}^{\mu }$ and $[\boldsymbol {\varLambda }']_{\mu \nu }=\varLambda '^{\mu }{}_{\nu }=\partial \hat {r}^{\mu }/\partial r^{\nu }$ are derived as

(A5a,b)\begin{equation} \boldsymbol{\varLambda}= \left( \begin{array}{@{}ccc@{}} 1 & \gamma(t) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right), \quad \boldsymbol{\varLambda}'= \left(\begin{array}{@{}ccc@{}} 1 & -\gamma(t) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right), \end{equation}

respectively, where $\boldsymbol {\varLambda }\boldsymbol {\cdot }\boldsymbol {\varLambda '}= \boldsymbol {\varLambda }'\boldsymbol {\cdot }\boldsymbol {\varLambda }=\boldsymbol{\mathsf{I}}$ by definition. Einstein's summation rule is applied hereafter.

By using transformation matrices, the covariant and contravariant basis vectors $\hat {\boldsymbol {E}}_{\mu }$ and $\hat {\boldsymbol {E}}^{\mu }$, respectively, and the corresponding components of the position vectors $\boldsymbol {r}=r_{\mu }\boldsymbol {e}^{\mu }=r^{\mu }\boldsymbol {e}_{\mu }=\hat {r}^{\mu } \hat {\boldsymbol {E}}_{\mu }=\hat {r}_{\mu }\hat {\boldsymbol {E}}^{\mu }$, are represented as

(A6a,b)\begin{gather} \hat{\boldsymbol{E}}_{\mu}=\varLambda{^{\nu}}{}_{\mu}\,\boldsymbol{e}_{\nu},\quad \hat{\boldsymbol{E}}^{\mu}=\varLambda{^{\mu}}{}_{\nu}\,\boldsymbol{e}^{\nu}, \end{gather}

(A7a,b)\begin{gather}\hat{r}_{\mu}=\varLambda{^{\nu}}{}_{\mu}\,r_{\nu},\quad \hat{r}^{\mu}=\varLambda{^{\mu}}{}_{\nu}\,r^{\nu}. \end{gather}

Since the oblique coordinate system is not an orthogonal system, covariant and contravariant bases are used, where $\hat {\boldsymbol {E}}_{\mu }\boldsymbol {\cdot }\hat {\boldsymbol {E}}^{\nu } =\delta _{\mu }{}^{\nu }$ holds. The lower and upper indices ($\mu ,\nu =1, 2, 3$) of the tensor variables represent the covariant and contravariant components of the tensor, respectively.

Figure 18 shows a schematic diagram of this transformation; a 2-D diagram on the shear plane is used for the sake of explanation. At $t=0$, the basis vectors of the oblique coordinates $\hat {\boldsymbol {E}}_1$ and $\hat {\boldsymbol {E}}_2$ coincide with those of the static Cartesian coordinates $\boldsymbol {e}_1$ and $\boldsymbol {e}_2$. At $t>0$, the second basis vector of the oblique coordinate changes with time. The contravariant metric tensor for the oblique coordinate is defined as

(A8)\begin{equation} G^{\mu\nu}\equiv\hat{\boldsymbol{E}}^{\mu}\boldsymbol{\cdot}\hat{\boldsymbol{E}}^{\nu} = \left( \begin{array}{@{}ccc@{}} 1+\gamma^2(t) & -\gamma(t) & 0 \\ -\gamma(t) & 1 & 0 \\ 0 & 0 & 1 \end{array} \right). \end{equation}

Note that, in the static Cartesian coordinate system, there is no distinction between the covariant and contravariant expressions, i.e. $r_{\mu }=r^{\mu }$ and $\boldsymbol {e}_{\mu }=\boldsymbol {e}^{\mu }$, and the metric tensor is identical to the unit tensor, i.e. $\boldsymbol{\mathsf{G}}=\boldsymbol{\mathsf{I}}$. In the oblique system in figure 18, where the coordinates are non-orthogonal but linear and spatially homogeneous, the metric tensor is time-varying and spatially constant. In this situation, the Christoffel term in the covariant differentiation is zero, and the covariant differentiation is represented by usual partial differentiation: $\hat {\boldsymbol {\nabla }}_{\mu }=\partial /\partial \hat {r}^{\mu }$.

Figure 18. Schematic diagram of the oblique coordinate system. Here $\hat {\boldsymbol {E}}_1$ and $\hat {\boldsymbol {E}}_2$ and $\boldsymbol {e}_1$ and $\boldsymbol {e}_2$ are the basis vectors in the oblique system and those in the static Cartesian system, respectively; and $r^1$ and $r^2$ are the components of the position vector in the Cartesian system. The Cartesian system and the initial oblique system coincide (blue square). At $t>0$, the oblique system is sheared (pink parallelogram) by the strain of $\dot {\gamma }t$, where $\dot {\gamma }$ is the applied shear rate. The region of the oblique system outside the initial square (right triangle) can be periodically transformed back into the square (left triangle).

The periodicity in the governing equations can also be achieved by only the coordinate transformation (A1)–(A4) with the orthogonal basis system (Rogallo Reference Rogallo1981; Onuki Reference Onuki1997) without using the oblique dual-basis system. This single-basis formalism has the advantage that the tensorial representation for an equation is expressed uniquely but has the disadvantage that a spatial differential operator includes the cross oblique term explicitly. In contrast, the dual-basis formalism adopted in our method has the advantage that the forms of differential operators and governing equations in the oblique coordinate system are almost the same as that in the orthogonal system, as explained in § 2, although these forms have dual (covariant and contravariant) expressions. This simple expression of the governing equations in the coordinate system is preferable for a convenient implementation of the practical simulation code.

A.2. Governing equations in the oblique coordinate system

The tensorial component representation of the fluid momentum equation on the general coordinate system (Luo & Bewley Reference Luo and Bewley2004; Venturi Reference Venturi2009; Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016) is

(A9)\begin{equation} \frac{\delta\hat{u}^{\mu}}{\delta \hat{t}}=\rho^{{-}1}\hat{\boldsymbol{\nabla}}_\nu\hat{\sigma}^{\nu\mu}+\hat{\phi}\hat{f}_p^{\mu}. \end{equation}

The left-hand side of (A9) is the intrinsic time derivative in the general coordinate system,

(A10)\begin{equation} \frac{\delta \hat{A}^{\mu}}{\delta \hat{t}}\equiv\frac{\partial \hat{A}^{\mu}}{\partial \hat{t}} +(\hat{u}^{\nu}-\hat{U}^{\nu})\hat{\boldsymbol{\nabla}}_{\nu}\hat{A}^{\mu} +\hat{A}^{\nu}\hat{\boldsymbol{\nabla}}_{\nu}\hat{U}^{\mu}, \end{equation}

where $\hat {U}^{\mu }\equiv -\partial \hat {r}^{\mu }/\partial t$ is the moving velocity of the coordinate, and for simple shear flow $\boldsymbol {U}=\dot {\gamma }(t)r^2\boldsymbol {e}_1= \dot {\gamma }(t)\hat {r}^2\hat {\boldsymbol {E}}_1$. Since $\hat {A}^{\mu }$ is defined in the moving system, the advection in the second term in (A10) is by the relative velocity to the coordinate flow. The last term in (A10) arises from the affine deformation caused by the coordinate flow. Introducing the relative fluid velocity to the coordinate flow $\boldsymbol {\xi }=\boldsymbol {u}-\boldsymbol {U}$, (A9) becomes

(A11)\begin{equation} (\hat{\partial}_{\hat{t}}+\hat{\xi}^{\nu}\hat{\boldsymbol{\nabla}}_{\nu}) \hat{\xi}^{\mu}=\rho^{{-}1}\hat{\boldsymbol{\nabla}}_{\nu}\hat{\sigma}^{\nu\mu}+ \hat{\phi}\hat{f}_p^{\mu}-2\dot{\gamma}(t)\hat{\xi}^2\delta^{\mu,1}, \end{equation}

with the incompressibility condition $\hat {\boldsymbol {\nabla }}_{\mu }\hat {\xi }^{\mu }=0$, where $\hat {\partial }_{\hat {t}}\equiv \partial /\partial \hat {t}|_{\hat {r}^{\mu }}=\partial /\partial t|_{r^{\mu }}+\dot {\gamma }(t)r^2\partial /\partial r^1$. The last term in (A11) arises from the spatial gradient of the coordinate flow. Since this equation does not explicitly depend on the coordinate components $\hat {r}^{\mu }$, periodic boundary conditions can be assigned to (A11), and hence (A11) can be solved by a spectral method (Rogallo Reference Rogallo1981; Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988). The stress tensor gradient in a Newtonian fluid is obtained as

(A12)\begin{equation} \hat{\boldsymbol{\nabla}}_{\nu}\hat{\sigma}_n^{\nu\mu}={-}G^{\nu\mu} \hat{\boldsymbol{\nabla}}_{\nu}\hat{p}+\eta_sG^{\nu\gamma} \hat{\boldsymbol{\nabla}}_{\nu}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{\xi}^{\mu}. \end{equation}

In a viscoelastic fluid, the polymer stress gradient term $\hat {\boldsymbol {\nabla }}_{\nu }\hat {\sigma }_p^{\nu \mu }$ is considered in addition to (A12). In our method, the tensorial expression for the constitutive equation of the polymer stress is additionally introduced in a manner consistent with the previous Newtonian formulation (Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016).

The intrinsic time derivative for conformation tensor $\boldsymbol{\mathsf{C}}= \hat {C}^{\mu \nu }\hat {\boldsymbol {E}}_{\mu }\hat {\boldsymbol {E}}_{\nu }$, which is represented by its second-rank contravariant tensor, is expressed as (Venturi Reference Venturi2009)

(A13)\begin{equation} \frac{\delta \hat{C}^{\mu\nu}}{\delta \hat{t}}\equiv\frac{\partial \hat{C}^{\mu\nu}}{\partial\hat{t}} +(\hat{u}^{\gamma}-\hat{U}^{\gamma})\hat{\boldsymbol{\nabla}}_{\gamma}\hat{C}^{\mu\nu} +\hat{C}^{\mu\gamma}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{U}^{\nu} +\hat{C}^{\gamma\nu}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{U}^{\mu}, \end{equation}

and the upper-convected time derivative is expressed by

(A14)\begin{equation} \frac{\textrm{d}_c\hat{C}^{\mu\nu}}{\textrm{d}\hat{t}}\equiv\frac{\delta \hat{C}^{\mu\nu}}{\delta\hat{t}}-\hat{C}^{\gamma\nu}\hat{\boldsymbol{\nabla}}_{\gamma} \hat{u}^{\mu}-\hat{C}^{\mu\gamma}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{u}^{\nu}. \end{equation}

Substituting (A13) into (A14), one obtains

(A15)\begin{equation} \frac{\textrm{d}_c\hat{C}^{\mu\nu}}{\textrm{d}\hat{t}}=\frac{\partial \hat{C}^{\mu\nu}}{\partial\hat{t}} +\hat{\xi}^{\gamma}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{C}^{\mu\nu} -\hat{C}^{\mu\gamma}\hat{\boldsymbol{\nabla}}_{\gamma} \hat{\xi}^{\nu}-\hat{C}^{\gamma\nu}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{\xi}^{\mu}. \end{equation}

By using (A15), the single-mode Oldroyd-B constitutive equation in the general coordinate system is represented as

(A16)\begin{gather} \frac{\partial \hat{C}^{\mu\nu}}{\partial\hat{t}} +\hat{\xi}^{\gamma}\hat{\boldsymbol{\nabla}}_{\gamma}\hat{C}^{\mu\nu} =\hat{C}^{\mu\gamma} \hat{\boldsymbol{\nabla}}_{\gamma}\hat{\xi}^{\nu} +\hat{C}^{\gamma\nu}\hat{\boldsymbol{\nabla}}_{\gamma} \hat{\xi}^{\mu}-\frac{1}{\lambda} (\hat{C}^{\mu\nu}-G^{\mu\nu} ), \end{gather}

(A17)\begin{gather}\hat{\sigma}_p^{\mu\nu}=\frac{\eta_p}{\lambda} (\hat{C}^{\mu\nu}-G^{\mu\nu} ). \end{gather}

Here, again, (A16) is independent of the coordinate components and has the same form as that in the orthogonal coordinate system. Therefore, periodic boundary conditions can be assigned to (A16).

B.1. Time-stepping algorithm for the coupling between fluid and particles

A flow chart showing the calculation procedure for one time step is shown in figure 19. The couplings between the flow and conformation and between the fluid and particles are established in the following explicit fractional step approach. Throughout the evolution process, field variables are converted from real space to wavenumber space and vice versa as necessary. In this section, continuum variables in the Fourier space are denoted by the subscript $\boldsymbol {k}$, where $\boldsymbol {k}$ represents the wavenumber vector. The discretized $n$th time step is indicated by the superscript of a variable as $(\cdot )^{n}$. Here, the constitutive equation is a single-mode Oldroyd-B model to ease explanation, and the extension to the multi-mode constitutive equations is straightforward. The calculation proceeds according to the following procedure.

1. (i) Initialization of variables (M-1). Starting with the coordinate strain $\gamma =0$, the field variables are initialized as $\boldsymbol {u}=\boldsymbol {\xi }=\boldsymbol {0}$ and $\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{I}}$ at $t=0$ over the entire domain. Correspondingly, the translational and angular velocities of the particles are set to zero. For a many-particle system, the positions of the particles are randomly generated to keep the distance between particle surfaces at least $2\varDelta$, where $\varDelta$ is the grid size.

   

   Figure 19. Flow chart of the main calculation procedure over one time step.

2. (ii) Update of the conformation tensor field (C-1). The conformation tensor $\hat {\boldsymbol{\mathsf{C}}}$ is updated to the next time step by integrating (A16) or (B31) over time to obtain the polymer stress field. As mentioned in § 2.2, the small error of $\hat {\boldsymbol{\mathsf{C}}}$ accumulates in the inner particle region according to the time evolution. This error can be eliminated by resetting $\hat {\boldsymbol{\mathsf{C}}}=\boldsymbol{\mathsf{G}}$ over the $\hat {\phi }=1$ region if necessary. In this study, this reset operation is safely omitted because the error is sufficiently small.

3. (iii) Update of the intermediate velocity field (V-1). Equation (A11) without the $\hat {\phi }\hat {\boldsymbol {f}}_p$ term is time-integrated to obtain an intermediate velocity field $\hat {\boldsymbol {\xi }}^*$. In this step, the fluid stress, i.e. the solvent and polymer stresses, are considered, and the solid–fluid coupling is not considered.

4. (iv) Update of the shear strain (M-2). After the field calculations, the shear strain is updated: $\gamma ^{n+1}=\gamma ^n+\dot {\gamma }\Delta t$, where $\Delta t$ is the time increment. This corresponds to the deformation of the oblique coordinate system. In this step, if the apparent shear strain equals the threshold value $\gamma _{{th}}$, the remeshing process is conducted. In this study, $\gamma _{{th}}=1$.

5. (v) Remeshing (M-3). Practically, as time evolves, the oblique mesh is gradually distorted, which can lead to a decrease in accuracy. To continue the simulation as the strain increases infinitely while maintaining accuracy, the strained oblique coordinates should be reset to a less strained state or to the static Cartesian coordinates at some finite shear strain (Rogallo Reference Rogallo1981). In this study, the oblique coordinate system is reset to the orthogonal Cartesian coordinate system when $\gamma$ reaches $\gamma _{{th}}=1$. First, the shear strain of the oblique system is reset as $\gamma \leftarrow \gamma -\gamma _{{th}}$. Then, the field variables $\hat {\boldsymbol {\xi }}^*$ and $\hat {\boldsymbol{\mathsf{C}}}$ on the oblique grid outside of the initial orthogonal grid are remapped through the periodic boundary in the flow direction. Simultaneously, the components of the variables in the oblique coordinate system are transformed to those in the reset coordinate system by the transformation matrix. Correspondingly, the metric tensor (A8) and the norm of the wavenumber vector in the spectral scheme, $\hat {\boldsymbol {k}}\boldsymbol {\cdot }\hat {\boldsymbol {k}}$, are updated. The norm of wavenumber vector in the wavenumber space corresponds to the Laplacian operator in real space, i.e. $G^{\mu \nu }\hat {\boldsymbol {\nabla }}_{\mu }\hat {\boldsymbol {\nabla }}_{\nu }\ \Longleftrightarrow \ -G^{\mu \nu }\hat {k}_{\mu }\hat {k}_{\nu }= -\hat {\boldsymbol {k}}\boldsymbol {\cdot }\hat {\boldsymbol {k}}$, and then

   (B1)\begin{equation} \hat{\boldsymbol{k}}\boldsymbol{\cdot}\hat{\boldsymbol{k}} =\hat{k}_1^2+(\hat{k}_2-\gamma\hat{k}_1)^2+\hat{k}_3^2, \end{equation}
   where $\hat {k}_i$ is the $i$th component of the covariant wavenumber vector in the oblique coordinate system.

6. (vi) Interpolation and transformation of the intermediate velocity field (M-4). To simplify the reconstruction of the $\phi$ field based on particle positions, the coupling between the fluid and particle is treated on the usual static orthogonal coordinate system. The grid points in the oblique coordinate system do not always coincide with those in the static orthogonal coordinate system. Therefore, the intermediate velocity field $\hat {\boldsymbol {\xi }}^*$ on the oblique grids should be interpolated to the static orthogonal grids. This is done by using a periodic cubic spline interpolation (Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016). After the interpolation, the oblique-basis components of $\hat {\boldsymbol {\xi }}^*$ are transformed to those in Cartesian basis. In this step, the absolute velocity field $\boldsymbol {u}^*$ is constructed using the transformed $\boldsymbol {\xi }^*$ and the base flow $\dot {\gamma }r^2\boldsymbol {e}_1$:

   (B2)\begin{equation} \boldsymbol{u}^*=(\dot{\gamma}r^2\delta^{\mu,1} +\varLambda^{\mu}_{\ \nu}\hat{\xi}^{\nu*})\boldsymbol{e}_{\mu}. \end{equation}

7. (vii) Update of the particle position (P-1). Hereafter, the calculation is conducted in the orthogonal coordinate system (blue block in figure 19). Using the particle velocity at the previous time step $\boldsymbol {V}^n$, the position of the $i$th particle is updated:

   (B3)\begin{equation} \boldsymbol{R}_i^{n+1}=\boldsymbol{R}_i^n +\int_{t^n}^{t^{n+1}}\boldsymbol{V}_i^n\,\textrm{d}t. \end{equation}
   In this step, if the updated particle position crosses the top and bottom boundaries, the position and velocity of the particle are modified according to the Lees–Edwards boundary conditions (Lees & Edwards Reference Lees and Edwards1972; Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011; Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016). In this time, the $\phi$ field is also updated by using the new particle positions consistent with Lees–Edwards boundary conditions. Then, the intermediate particle velocity field $\boldsymbol {u}_p^*$ is calculated:
   (B4)\begin{equation} \phi^{n+1}\boldsymbol{u}_p^*=\sum_i\phi_i^{n+1}[\boldsymbol{V}_i^n +\boldsymbol{\varOmega}_i^n\times\boldsymbol{r}_i^{n+1}], \end{equation}
   where $\boldsymbol {r}_i=\boldsymbol {r}-\boldsymbol {R_i}$. This corresponds to the mapping of the Lagrangian particle velocity on the Euler velocity field.

8. (viii) Calculation of hydrodynamic forces acting on particles (M-5). The hydrodynamic force and torque exerted on the particles $\boldsymbol {F}^H$ and $\boldsymbol {N}^H$ are calculated by the change in momentum in the particle domain:

   (B5)\begin{gather} \int_{t^n}^{t^{n+1}}\boldsymbol{F}_i^H\,\textrm{d}t= \int\rho\phi_i^{n+1}(\boldsymbol{u}^*-\boldsymbol{u}_p^*)\,\textrm{d}\boldsymbol{r}, \end{gather}
   (B6)\begin{gather}\int_{t^n}^{t^{n+1}}\boldsymbol{N}_i^H\,\textrm{d}t= \int\boldsymbol{r}_i^{n+1}\times\rho\phi_i^{n+1} (\boldsymbol{u}^*-\boldsymbol{u}_p^*)\,\textrm{d}\boldsymbol{r}. \end{gather}

9. (ix) Update of the particle velocity and angular velocity (P-2). Using (B5) and (B6), the particle velocities are updated as

   (B7)\begin{gather} \boldsymbol{V}_i^{n+1}=\boldsymbol{V}_i^n+\frac{1}{M_i}\int_{t^n}^{t^{n+1}} [\boldsymbol{F}_i^H+\boldsymbol{F}_i^C ]\,\textrm{d}t, \end{gather}
   (B8)\begin{gather}\boldsymbol{\varOmega}_i^{n+1}=\boldsymbol{\varOmega}_i^n+ \boldsymbol{\mathsf{I}}_{p,i}^{{-}1} \boldsymbol{\cdot}\int_{t^n}^{t^{n+1}}\boldsymbol{N}_i^H\,\textrm{d}t. \end{gather}
   In this study, for the inter-particle force $\boldsymbol {F}^C$, the soft-core (truncated Lennard-Jones) potential, which produces the short-range repulsive force, is adopted:
   (B9)\begin{gather} \boldsymbol{F}^C_i(\boldsymbol{R})={-}\sum_{j\neq i}^N F_{{soft}}(r_{ij})\frac{\boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij}|}, \end{gather}
   (B10)\begin{gather}F_{{soft}}(r_{ij})={-}\left(\frac{\partial U_{{soft}}}{\partial r}\right)_{r=r_{ij}}, \end{gather}
   (B11)\begin{gather}U_{{soft}}(r)=\begin{cases} 4\epsilon\left[\left(\dfrac{2a}{r}\right)^{36}-\left(\dfrac{2a}{r}\right)^{18}\right] +\epsilon & (r < r_c), \\ 0 & (r\geqslant r_c), \end{cases} \end{gather}
   where $\boldsymbol {r}_{ij}=\boldsymbol {R}_j-\boldsymbol {R}_i$ is the distance vector from the $i$th particle to the $j$th particle and $r_c=2^{1/18}(2a)$. Vector $\boldsymbol {r}_{ij}$ is modified according to periodic boundary conditions if necessary. This potential force is simply applied to avoid particle overlap. The force parameter $\epsilon$, which tunes the interaction strength, is set at $\epsilon /(\eta _0\dot {\gamma }a^3)=0.561$ in all many-particle calculations in this study. Under denser particle concentration conditions, where the particle collisions and/or friction and its contribution to the total stress can become significant, more realistic modelling of inter-particle force may be required.

10. (x) Calculation of the coupling term for fluid (M-6). Now that both the positions and velocities of particles have been updated, the final particle velocity field $\boldsymbol {u}_p$ is obtained as

    (B12)\begin{equation} \phi^{n+1}\boldsymbol{u}_p^{n+1}=\sum_i\phi_i^{n+1} [\boldsymbol{V}_i^{n+1}+\boldsymbol{\varOmega}_i^{n+1}\times\boldsymbol{r}_i^{n+1} ]. \end{equation}
    Then, the body force $\phi \boldsymbol {f}_p$ is calculated as
    (B13)\begin{equation} \int_{t^n}^{t^{n+1}}\phi\boldsymbol{f}_p(\boldsymbol{x},t)\,\textrm{d}t= \phi^{n+1}(\boldsymbol{u}_p^{n+1}-\boldsymbol{u}^*). \end{equation}
    To calculate the stresslet (2.19), (B13) is further transformed as
    (B14)\begin{equation} \int_{t^n}^{t^{n+1}}\phi\boldsymbol{f}_p\,\textrm{d}t= \phi^{n+1}(\boldsymbol{u}_p^{n+1}-\boldsymbol{u}^*_p) -\phi^{n+1}(\boldsymbol{u}^*-\boldsymbol{u}_p^*). \end{equation}
    The first term on the right-hand side is expressed by the changes in particle velocity $\Delta \boldsymbol {V}_i=\boldsymbol {V}_i^{n+1}-\boldsymbol {V}_i^n$ from (B7) and in angular velocity $\Delta \boldsymbol {\varOmega }_i=\boldsymbol {\varOmega }_i^{n+1} -\boldsymbol {\varOmega }_i^n$ from (B8), as
    (B15)\begin{equation} \phi^{n+1}(\boldsymbol{u}_p^{n+1}-\boldsymbol{u}_p^*)=\sum_i\phi_i^{n+1} [\Delta\boldsymbol{V}_i^H+\Delta\boldsymbol{\varOmega}_i^H\times\boldsymbol{r}_i^{n+1}] +\sum_i\phi_i^{n+1}\Delta\boldsymbol{V}_i^C, \end{equation}
    where $\Delta \boldsymbol {V}_i^H$, $\Delta \boldsymbol {\varOmega }_i^H$ and $\Delta \boldsymbol {V}_i^C$ are the updates by the hydrodynamic forces $\boldsymbol {F}_i^H$ and $\boldsymbol {N}_i^H$ and the inter-particle force $\boldsymbol {F}_i^C$, respectively. Therefore, $\phi \boldsymbol {f}_p$ can be decomposed into the individual contributions from the hydrodynamic interactions $\phi \boldsymbol {f}_p^H$ and direct inter-particle interactions $\phi \boldsymbol {f}_p^C$, as $\phi \boldsymbol {f}_p=\phi \boldsymbol {f}_p^H+\phi \boldsymbol {f}_p^C$. The individual contributions of the body force are expressed to first order in time as
    (B16)\begin{gather} \phi\boldsymbol{f}_p^H\Delta t=\sum_i\phi_i^{n+1}[\Delta\boldsymbol{V}_i^H+ \Delta\boldsymbol{\varOmega}_i^H\times\boldsymbol{r}_i^{n+1}] -\phi^{n+1}(\boldsymbol{u}^*-\boldsymbol{u}_p^*), \end{gather}
    (B17)\begin{gather}\phi\boldsymbol{f}_p^C\Delta t= \sum_i\phi_i^{n+1}\Delta\boldsymbol{V}_i^C. \end{gather}
    In the calculation of the stresslet contribution from $\phi \boldsymbol {f}_p^C$, the direct virial expression was used instead of (B17) for computational efficiency (Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016):
    (B18)\begin{equation} \boldsymbol{\mathsf{S}}^C={-}\frac{1}{N}\int_{D_V} \boldsymbol{r}\rho\phi\boldsymbol{f}_p^C\textrm{d}\boldsymbol{r} ={-}\frac{1}{N}\sum_{i < j}\boldsymbol{r}_{ij}\boldsymbol{F}_{ij}^C, \end{equation}
    where $\boldsymbol {F}_{ij}^C$ is the inter-particle force on the $i$th particle due to the $j$th particle. In this study, for conditions up to $\phi _p=0.1$ and ${Wi}=2.5$, the contribution of ${\mathsf{S}}_{12}^C$ to the total stresslet ${\mathsf{S}}_{12}$ is small compared to the hydrodynamic contributions ($\langle {\mathsf{S}}_{12}^C\rangle /\langle {\mathsf{S}}_{12}\rangle <0.05$ at $\phi _p=0.1$).

11. (xi) Update of the velocity field (V-2). Finally, the integrated body force $\phi \boldsymbol {f}_p$ is remapped and transformed from the orthogonal coordinate system to the oblique coordinate system and added to the intermediate velocity field $\hat {\boldsymbol {\xi }}^*$:

    (B19)\begin{equation} \hat{\boldsymbol{\xi}}^{n+1}=\hat{\boldsymbol{\xi}}^* + \widehat{\left[\int_{t^n}^{t^{n+1}}\phi\boldsymbol{f}_p\,\textrm{d}s\right]}. \end{equation}
    At this stage, incompressibility is assigned in the Fourier space,
    (B20)\begin{equation} \hat{\boldsymbol{\xi}}_{\boldsymbol{k}}\leftarrow \hat{\boldsymbol{\xi}}_{\boldsymbol{k}}- \frac{(\hat{\xi}_{\boldsymbol{k}}\boldsymbol{\cdot}\hat{\boldsymbol{k}})} {\hat{\boldsymbol{k}}\boldsymbol{\cdot}\hat{\boldsymbol{k}}}\hat{\boldsymbol{k}}. \end{equation}
    This solenoidal projection is also adopted after calculating $\hat {\boldsymbol {\xi }}^*$ (V-1).

12. The described fractional steps are repeated until the calculated time reaches the target final time. Further information about the time-stepping algorithm is detailed in previous work (Nakayama et al. Reference Nakayama, Kim and Yamamoto2008; Molina et al. Reference Molina, Otomura, Shiba, Kobayashi, Sano and Yamamoto2016).

B.2. Spatial discretization and time integral scheme

Since the periodic boundary conditions are assigned in each direction, the continuum variables such as $\hat {\boldsymbol {\xi }}$ and $\hat {\boldsymbol{\mathsf{C}}}$ are Fourier-transformed. In real space, the continuum variables are collocated on the equispaced mesh point with spacing $\varDelta$. Spatial derivatives are calculated in Fourier space while the second-order terms like the advection term in (A11) and (A16) are calculated by a transformation method (Orszag Reference Orszag1969).

For the integration of $\hat {\boldsymbol {\xi }}_{\boldsymbol {k}}$ over time, the exact linear part method, which is preferred for solving stiff equations (Beylkin, Keiser & Vozovoi Reference Beylkin, Keiser and Vozovoi1998), is adopted, where the nonlinear part is discretized by the Euler method.

For the polymer constitutive equation, the explicit Euler method (first-order) is adopted. In this study, to evaluate a single-particle system that corresponds to very dilute suspensions ($\phi _p\sim 0.001$), the discretized equation (A16) is solved directly. This naive implementation has been stable and accurate in such dilute conditions. However, at high $\phi _p$ and ${Wi}$ conditions, the large growth rate of the polymer stress around the particles violates the positive-definiteness of the conformation tensor, thus resulting in an inaccurate solution or divergence. Therefore, to evaluate a many-particle system that corresponds to semidilute suspensions ($0.025\leqslant \phi _p\leqslant 0.1$), the log-conformation formalism is used in which the time evolution equation of $\log \boldsymbol{\mathsf{C}}$ rather than $\boldsymbol{\mathsf{C}}$ is solved to guarantee the positive-definiteness of $\boldsymbol{\mathsf{C}}$ (Fattal & Kupferman Reference Fattal and Kupferman2004; Hulsen, Fattal & Kupferman Reference Hulsen, Fattal and Kupferman2005). A detailed description of the log-conformation formalism is provided in appendix B.3.

The particle position is updated by discretizing (B3) by the Euler method (first-order) at the first time step and the second-order Adams–Bashforth scheme later. In the update of the particle velocity (B7) and (B8), the impulsive hydrodynamic force and torque are calculated by applying (B5) and (B6), respectively, and the potential force $\boldsymbol {F}^C$ in (B7) is discretized by the second-order Heun scheme because both particle positions at $t^n$ and $t^{n+1}$ are already obtained in that stage:

(B21)\begin{equation} \frac{1}{M_i}\int_{t^n}^{t^{n+1}}\boldsymbol{F}_i^C\textrm{d}s =\frac{\Delta t}{2M_i} [\boldsymbol{F}_i^C(\boldsymbol{R}^{n+1})+\boldsymbol{F}_i^C(\boldsymbol{R}^n) ]. \end{equation}

The time increment is determined based on the stability given by the fluid momentum diffusion: $\Delta t=\rho /\eta _0K_{max}^2$ ($K_{max}$ is the largest wavenumber in the spectral scheme). As proven in the code validations in § 3.1 and appendix C, this choice is reasonable considering the conditions in this study.

B.3. Log-conformation-based constitutive equation for Oldroyd-B model

Because the conformation tensor $\boldsymbol{\mathsf{C}}$ is real-symmetric and positive-definite, $\boldsymbol{\mathsf{C}}$ can be diagonalized as

(B22)\begin{equation} \boldsymbol{\mathsf{C}}= \boldsymbol{\mathsf{R}}\boldsymbol{\cdot}\boldsymbol{\varLambda} \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^\textrm{T}, \end{equation}

where $\boldsymbol {\varLambda }=\textrm {{diag}}(\lambda _1,\lambda _2,\lambda _3)$ and $\lambda _i>0\ (i=1,2,3)$ are the eigenvalues of $\boldsymbol{\mathsf{C}}$, and $\boldsymbol{\mathsf{R}}$ is the rotation matrix composed of the eigenvectors of $\boldsymbol{\mathsf{C}}$. Here, the new tensor variable $\boldsymbol {\varPsi }$ is introduced (Fattal & Kupferman Reference Fattal and Kupferman2004) as

(B23)\begin{equation} \boldsymbol{\varPsi}= \boldsymbol{\mathsf{R}}\boldsymbol{\cdot} \boldsymbol{\varLambda}_{\varPsi}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^{\rm T}, \end{equation}

where $\boldsymbol {\varLambda }_{\varPsi }=\textrm {{diag}}(\ln \lambda _1,\ln \lambda _2,\ln \lambda _3) =\ln \boldsymbol {\varLambda }$. From (B22) and (B23),

(B24)\begin{equation} \boldsymbol{\mathsf{C}}= \boldsymbol{\mathsf{R}}\boldsymbol{\cdot}\exp (\boldsymbol{\varLambda}_{\varPsi})\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^\textrm{T}. \end{equation}

Note that, when $\boldsymbol{\mathsf{C}}$ is obtained through (B24), $\boldsymbol{\mathsf{C}}$ is strictly positive-definite by definition. Furthermore, utilizing the time evolution of $\boldsymbol {\varPsi }$ instead of (2.4), the exponential growth in $\boldsymbol{\mathsf{C}}$ is translated to the linear growth of $\boldsymbol {\varPsi }$, which enables numerical stability in the time evolution. Specifically, the stretching in the principal axes of $\boldsymbol{\mathsf{C}}$ by the velocity gradient tensor $\boldsymbol {\nabla }\boldsymbol {u}$ is extracted as

(B25)\begin{gather} \boldsymbol{\nabla}\boldsymbol{u}= \boldsymbol{\mathsf{R}}\boldsymbol{\cdot} \boldsymbol{\mathsf{M}}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^\textrm{T}, \end{gather}

(B26)\begin{gather}\boldsymbol{\mathsf{B}}= \boldsymbol{\mathsf{R}}\boldsymbol{\cdot}\textrm{diag}(M_{11},M_{22},M_{33})\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^\textrm{T}, \end{gather}

where $\boldsymbol{\mathsf{B}}$ is symmetric and commutes with $\boldsymbol{\mathsf{C}}$ by definition.

The residual component $\boldsymbol {\nabla }\boldsymbol {u}-\boldsymbol{\mathsf{B}}$ can be decomposed as

(B27)\begin{equation} \boldsymbol{\nabla}\boldsymbol{u}- \boldsymbol{\mathsf{B}} = \boldsymbol{\mathsf{A}}+ \boldsymbol{\mathsf{C}}^{{-}1}\boldsymbol{\cdot} \boldsymbol{\mathsf{N}}, \end{equation}

with antisymmetric tensors $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{N}}$ (Fattal & Kupferman Reference Fattal and Kupferman2004). Tensor $\boldsymbol{\mathsf{N}}$ is proven to be irrelevant in the affine deformation of $\boldsymbol{\mathsf{C}}$ by inserting (B27) into the upper-convected time derivative of $\boldsymbol{\mathsf{C}}$. On the other hand, $\boldsymbol{\mathsf{A}}$ represents the rotation of the principal axes of $\boldsymbol{\mathsf{C}}$. From the affine deformation of $\boldsymbol{\mathsf{C}}$ in (2.4), the explicit expression of $\boldsymbol{\mathsf{A}}$ in the frame of the principal axes of $\boldsymbol{\mathsf{C}}$ is derived by Hulsen et al. (Reference Hulsen, Fattal and Kupferman2005) as

(B28)\begin{equation} {\mathsf{A}}_{ij}=\frac{\lambda_i{\mathsf{M}}_{ij}+\lambda_j{\mathsf{M}}_{ji}} {\lambda_i-\lambda_j},\quad i\neq j,\quad \lambda_i\neq\lambda_j \end{equation}

(the summation convention is not applied here). When $\lambda _i=\lambda _j$, ${\mathsf{A}}_{ij}$ is not uniquely determined in the decomposition of $\boldsymbol {\nabla }\boldsymbol {u}$ in (B27), but the affine deformation of $\boldsymbol{\mathsf{C}}$ and $\boldsymbol {\varPsi }$ by $\boldsymbol {\nabla }\boldsymbol {u}$ is still well defined, which case is explained next.

By using these tensors, the governing equation of $\boldsymbol {\varPsi }$ for the single-mode Oldroyd-B model is expressed as

(B29)\begin{equation} \frac{\partial \boldsymbol{\varPsi}}{\partial t} +\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varPsi} ={-} \boldsymbol{\mathsf{A}}\boldsymbol{\cdot}\boldsymbol{\varPsi} +\boldsymbol{\varPsi}\boldsymbol{\cdot} \boldsymbol{\mathsf{A}}+2 \boldsymbol{\mathsf{B}} + \boldsymbol{\mathsf{R}}\boldsymbol{\cdot} \left[\frac{1}{\lambda}(\boldsymbol{\varLambda}^{{-}1}- \boldsymbol{\mathsf{I}})\right] \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}^\textrm{T}. \end{equation}

When $\lambda _i=\lambda _j$, the corotational terms including ${\mathsf{A}}_{ij}$ are reduced as (in the frame of the principal axes of $\boldsymbol{\mathsf{C}}$)

(B30)\begin{align} -{\mathsf{A}}_{ij}\varPsi_{jj}+\varPsi_{ii}{\mathsf{A}}_{ij} &= (\ln\lambda_i-\ln\lambda_j ) \frac{\lambda_i{\mathsf{M}}_{ij}+\lambda_j{\mathsf{M}}_{ji}}{\lambda_i-\lambda_j}\nonumber\\ &\to {\mathsf{M}}_{ij}+{\mathsf{M}}_{ji}=2{\mathsf{D}}_{ij} \end{align}

(the summation convention is not applied here). With this treatment, the evolution equation (B29) of $\boldsymbol {\varPsi }$ works safely even when $\lambda _i=\lambda _j$ happens.

In the initial conditions, $\boldsymbol {u}=\boldsymbol {0}$ over the entire domain leads to $\boldsymbol{\mathsf{A}}=\boldsymbol {0}$ and $\boldsymbol{\mathsf{B}}=\boldsymbol{\mathsf{D}}$, and $\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{I}}$ results in $\boldsymbol{\mathsf{R}}=\boldsymbol{\mathsf{I}}$, $\boldsymbol {\varLambda }=\boldsymbol{\mathsf{I}}$ and $\boldsymbol {\varPsi }=\boldsymbol {0}$. The evolution of $\boldsymbol {\varPsi }$ according to (B29) is solved by numerical simulation; and $\boldsymbol{\mathsf{C}}$ is calculated from $\boldsymbol {\varPsi }$ via (B23) and (B24). The contravariant tensor expression corresponding to (B29) in the oblique coordinates is the following:

(B31)\begin{align} \frac{\partial\hat{\varPsi}^{\mu\nu}}{\partial \hat{t}}+\hat{\xi}^{\gamma}\hat{ \boldsymbol{\nabla}}_{\gamma}\hat{\varPsi}^{\mu\nu}&={-}\dot{\gamma} (\hat{\varPsi}^{\mu2}\delta^{1,\nu}+\hat{\varPsi}^{2\nu} \delta^{1,\mu} )-{\mathsf{G}}_{\gamma\zeta}\hat{A}^{\mu\gamma}\hat{\varPsi}^{\zeta\nu}+ {\mathsf{G}}_{\gamma\zeta}\hat{\varPsi}^{\mu\gamma}\hat{A}^{\zeta\nu}\nonumber\\ &\quad +2\hat{B}^{\mu\nu}+\hat{R}^{\mu\gamma}\left[\frac{1}{\lambda}( [\hat{\boldsymbol{\varLambda}}]^{{-}1}_{\gamma\zeta} -{\mathsf{G}}_{\gamma\zeta})\right]\hat{R}^{\nu\zeta}, \end{align}

where $[\hat {\boldsymbol {\varLambda }}]^{-1}_{\gamma \zeta }$ represents the covariant matrix component of $\hat {\boldsymbol {\varLambda }}^{-1}$, which is simply the matrix inverse of the contravariant matrix $\hat {\varLambda }^{\gamma \zeta }$; and ${\mathsf{G}}_{\gamma \zeta }$ is the covariant metric tensor, which is defined as

(B32)\begin{equation} {\mathsf{G}}_{\gamma\zeta}\equiv\hat{\boldsymbol{E}}_{\gamma} \boldsymbol{\cdot}\hat{\boldsymbol{E}}_{\zeta}=\left( \begin{array}{@{}ccc@{}} 1 & \gamma(t) & 0 \\ \gamma(t) & 1+\gamma^2(t) & 0 \\ 0 & 0 & 1 \end{array} \right). \end{equation}

Note that, in the oblique coordinates, there is an additional term originating from the moving coordinates (the first term on right-hand side of (B31)). Since (B31) does not explicitly depend on the coordinate variables, it can be discretized by a spectral method.

C.1. Polymer stress around a single particle

To test the validity of the developed method, a single-particle system is set up where a neutrally buoyant spherical particle is suspended in a sheared Oldroyd-B fluid (figure 1a). The cubic domain with a box length of $L$ is sufficiently large compared to the size of the particle used to represent the dilute particle system. Hereafter, for simplicity, the directions of the Cartesian coordinate basis vectors are denoted by $x$, $y$ and $z$ instead of the 1, 2 and $3$ notation used in appendices A and B, where $x$, $y$ and $z$ indicate the flow, velocity-gradient and vorticity directions, respectively. As shown in figure 1(a), because the particle is located at the centre of a simple shear flow, the net translational hydrodynamic force acting on the particle $\boldsymbol {F}^H$ vanishes, while the hydrodynamic torque $\boldsymbol {N}^H$ rotates the particle.

A flow condition is considered in the $\beta \to 1$, small-${Wi}$ and small-${Re}$ limit where an analytical solution is available. In this limit, the flow pattern is minimally affected by polymer stress, which is expressed analytically (Lin, Peery & Schowalter Reference Lin, Peery and Schowalter1970; Mikulencak & Morris Reference Mikulencak and Morris2004). Furthermore, when ${Wi}\ll 1$, the polymer stress distribution is approximated by the second-order fluid (SOF) theory: $\boldsymbol {\sigma }_p=2\eta _p\boldsymbol{\mathsf{D}}+4(\varPsi _1+\varPsi _2)\boldsymbol{\mathsf{D}} \boldsymbol {\cdot }\boldsymbol{\mathsf{D}} -\varPsi _1(\boldsymbol{\mathsf{D}}_{(2)}+4\boldsymbol{\mathsf{D}}\boldsymbol {\cdot } \boldsymbol{\mathsf{D}})$, where $\boldsymbol{\mathsf{D}}_{(2)}$ is the upper-convected derivative of $\boldsymbol{\mathsf{D}}$ (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987). Considering an Oldroyd-B fluid ($\varPsi _1=2\eta _p\lambda ,\ \varPsi _2=0$), the normalized polymer stress in the SOF limit is expressed as $\boldsymbol {\sigma }_p/(\eta _p\dot {\gamma })= 2\tilde {\boldsymbol{\mathsf{D}}}-2{Wi}\tilde {\boldsymbol{\mathsf{D}}}_{(2)}$. Here, $\beta =0.99$, ${Wi}=0.001$ and ${Re}=0.0142$. In this situation, the normalized polymer stress is approximated by $\boldsymbol {\sigma }_p/(\eta _p\dot {\gamma })\approx 2\tilde {\boldsymbol{\mathsf{D}}}$.

Different mesh resolutions of the particle interface are examined: $(a/\varDelta , \xi /\varDelta )= (8, 2)$, $(8, 1)$ and $(20, 2)$. At first, the overall trend of the polymer shear stress ($\sigma _{p,xy}$) distribution is similar for different resolutions (figure 20a–c) and the analytical solution (figure 20d). The only difference is that small $\sigma _{p,xy}$ oscillation is observed in the numerical solutions. When comparing the result for $a=8\varDelta$ and $\xi =2\varDelta$ (figure 20a) with that for $a=20\varDelta$ and $\xi =2\varDelta$ (figure 20c), it is clear that, as $a/\varDelta$ increases, the wavenumber of the small ripple in $\sigma _{p,xy}$ increases, but its amplitude decreases. This is due to the slow convergence of the Fourier series caused by the discontinuous change in $\sigma _{p,xy}$ at the solid–liquid interface. This artifact is a partly unavoidable intrinsic property of the spectral method. Regarding the interface thickness, when comparing the result for $a=8\varDelta$ and $\xi =2\varDelta$ (figure 20a) with that for $a=8\varDelta$ and $\xi =\varDelta$ (figure 20b), no significant difference in the overall trend is observed. However, when $\xi =\varDelta$, the distribution of ${\sigma }_{p,xy}$ near the interface is somewhat blurred, which is caused by the decrease in the number of mesh points that support the interface region.

Figure 20. Mesh resolution dependence of polymer shear stress distributions near the particle. In panels (a–d) the polymer shear stress normalized by $\eta _p\dot {\gamma }$ on the shear plane through the centre of the particle is drawn by colour contour for $(a/\varDelta ,\xi /\varDelta )=(8,2)$, $(8,1)$ and $(20,2)$ and for the analytical solution, respectively. The dotted lines around the particles for the DNS results in panels (a–c) denote the radial location of $a+\xi /2$.(It may be necessary to zoom-in to see these.) In panel (e), the normalized polymer shear stress along the line from the particle centre to the shear gradient direction $y$ is drawn: the exact solution (analytical) of $\boldsymbol {\sigma }_p/(\eta _p\dot {\gamma }) \approx 2\tilde {\boldsymbol{\mathsf{D}}}-2{Wi}\tilde {\boldsymbol{\mathsf{D}}}_{(2)}$ at ${Wi}\to 0$ (solid line), and the results for $(a/\varDelta ,\xi /\varDelta )=(20,2)$ (red), $(8,2)$ (green) and $(8,1)$ (blue), respectively. In panel (e), the interface region indicated by $0<\phi <1$ at around $r_y/a=1$ with thickness $\xi$ is coloured in the same manner as that for the symbols.

For a detailed evaluation, a one-dimensional profile of the polymer shear stress in the velocity-gradient direction from the particle centre is shown in figure 20(e). The steep increase in $\sigma _{p,xy}$ near the particle surface is reasonably reproduced as the mesh resolution increases, though the peak in $\sigma _{p,xy}$ is somewhat smeared due to the limited mesh points in the interface domain. Hereafter, considering a balance between the accuracy of the numerical solution and the required computational cost, the particle radius is set to $a=8\varDelta$ and the interface thickness to $\xi =2\varDelta$. As seen in appendix C.2 and § 3, this resolution is sufficiently valid for the problems investigated in this study.

C.2. Rotation of a particle under simple shear flow

Under simple shear in Stokes flow, as is well known, a suspended particle in a Newtonian medium rotates with an angular velocity that is half of the applied shear rate, $\omega _z/\dot {\gamma }=0.5$, where $\omega _z$ is the particle angular velocity in the vorticity direction. However, in viscoelastic fluids, this relative rotational speed decreases with increasing ${Wi}$ (Snijkers et al. Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2009, Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2011). D'Avino et al. (Reference D'Avino, Hulsen, Snijkers, Vermant, Greco and Maffettone2008) and Snijkers et al. (Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2009, Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2011) conducted numerical evaluation for this phenomenon using a finite element method (FEM) and surface-conforming mesh, reproducing the experimental rotational slowdown data with DNS. They observed that the distribution of local torque and pressure on the particle surface becomes asymmetrical with ${Wi}$. However, the physics of slowdown has not been elucidated. This result is often referred to as the benchmark problem for a newly developed numerical scheme of viscoelastic suspensions (Ji et al. Reference Ji, Jiang, Winkler and Gompper2011; Yang et al. Reference Yang, Krishnan and Shaqfeh2016; Vázquez-Quesada & Ellero Reference Vázquez-Quesada and Ellero2017; Fernandes et al. Reference Fernandes, Faroughi, Carneiro, Nóbrega and McKinley2019). To validate the method developed in this work, the angular velocity of a particle in an Oldroyd-B fluid is evaluated at the same numerical conditions as previously reported (Snijkers et al. Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2009); however, no walls are used in this study. The numerical set-up is the same as that in § 3.1.

Figures 21(a) and 21(b) show the $\beta$ and ${Wi}$ dependence of the normalized particle angular velocity $\omega _z/\dot {\gamma }$. To compare with the previous numerical result, figure 21(a) shows the result at $\beta =0.615$. As ${Wi}$ increases, $\omega _z/\dot {\gamma }$ decreases. At ${Wi}\lesssim 1$, the result converges with the theoretical prediction up to $O(Wi^4)$ made using asymptotic methods (Housiadas & Tanner Reference Housiadas and Tanner2011a,b, Reference Housiadas and TannerReference Housiadas and Tanner2018): $\omega _z/\dot {\gamma }=1/2-(1-\beta ){Wi}^2/[4(1-4{Wi}^2\tilde {\varOmega }_4)]$, where $\tilde {\varOmega }_4$ is the coefficient of the $(1-\beta ){Wi}^4$ term in the series solution. At ${Wi}\gtrsim 1$, the asymptotic prediction starts to overestimate $\omega _z/\dot {\gamma }$. In this region, the result agrees reasonably well with the previous FEM result (Snijkers et al. Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2009). Figure 21(b) shows the $\beta$ dependence of $\omega _z/\dot {\gamma }$. As $\beta$ decreases, which corresponds to an increase in the polymer stress contribution, the negative slope of $\omega _z/\dot {\gamma }$ increases. This trend is consistent with the asymptotic predictions shown in figure 21(b).

Figure 21. The ${Wi}$ dependence of the normalized particle angular velocity $\omega _z/\dot {\gamma }$. (a) The result of this work (red circles) compared with a previous numerical result (solid line; Snijkers et al. Reference Snijkers, D'Avino, Maffettone, Greco, Hulsen and Vermant2009) and the asymptotic solution (dashed line; Housiadas & Tanner Reference Housiadas and Tanner2011a,b, Reference Housiadas and Tanner2018) at $\beta =0.615$. (b) The $\beta$ dependence of $\omega _z/\dot {\gamma }$ for: $\beta =0.2$ (red circles), $0.5$ (orange squares), $0.615$ (green triangles) and $0.8$ (blue diamonds). The lines correspond to predictions of the asymptotic solution.

The slowdown of rotation with increasing ${Wi}$ and/or $(1-\beta )$ suggests that the energy partition from external work to elastic energy increases. Figure 22 shows the normalized energy dissipation rate around the particle on the shear plane through the particle centre at $\beta =0.5$ and ${Wi}=0.1$ and $1.0$. The dissipation rate is decomposed into viscous ($\varPhi _{s}$) and elastic ($\varPhi _{p}$) contributions (Vázquez-Quesada et al. Reference Vázquez-Quesada, Espa nol, Tanner and Ellero2019) as $\varPhi _{{t}}=\varPhi _{{s}}+\varPhi _{{p}}$:

(C1)\begin{gather} \varPhi_{s}=2\eta_s \boldsymbol{\mathsf{D}}\boldsymbol{:} \boldsymbol{\mathsf{D}}, \end{gather}

(C2)\begin{gather}\varPhi_{p}=\frac{\eta_p}{2\lambda^2} (\textrm{tr}\, \boldsymbol{\mathsf{C}}+\textrm{tr}\, \boldsymbol{\mathsf{C}}^{{-}1}-6 ). \end{gather}

At ${Wi}=0.1$, $\varPhi _{s}$ (figure 22a) and $\varPhi _{p}$ (figure 22b) present similar distributions. Since $\boldsymbol {\sigma }_p= (\eta _p/\lambda )(\boldsymbol{\mathsf{C}}-\boldsymbol{\mathsf{I}})\to 2\eta _p\boldsymbol{\mathsf{D}}$ at ${Wi}\to 0$, $\varPhi _{p}$ is reduced to $2\eta _p\boldsymbol{\mathsf{D}}\boldsymbol {:}\boldsymbol{\mathsf{D}}$, which is proportional to $\varPhi _{s}$. In contrast, at ${Wi}=1.0$, $\varPhi _{s}$ and $\varPhi _{p}$ develop differently; the high-$\varPhi _{p}$ region expands, whereas $\varPhi _{s}$ does not change much in comparison to the ${Wi}=0.1$ case. The distribution of $\varPhi _{p}$ expands towards the shear-flow direction and high-$\varPhi _{p}$ grows near the equator of the particle (figure 22d), thus clearly showing an increase in the fraction of elastic energy dissipation at high ${Wi}$. Particle rotation is caused by viscous stress. An increase in elastic energy leads to a decrease in the relative fraction of viscous dissipation. As a result, the angular velocity of the particle in the viscoelastic medium decreases in comparison with that in viscous media. This slowdown of the particle rotation is enhanced with ${Wi}$ and $(1-\beta )$.

Figure 22. The ${Wi}$ dependence of the normalized energy dissipation around a particle at $\beta =0.5$: (a,b) ${Wi}=0.1$ and (c,d) ${Wi}=1.0$. Panels (a,c) and (b,d) show the distributions of viscous dissipation $\varPhi _{s}/(\eta _0\dot {\gamma }^2)$ and polymer dissipation $\varPhi _{p}/(\eta _0\dot {\gamma }^2)$, respectively. The dotted lines around particles show the radial location of $a+\xi /2$. (Again, it may be necessary to zoom-in to see these.)

In this section, the agreement between the numerical results of this work and the previously reported numerical and theoretical results verifies that the presented numerical scheme can successfully capture the dynamic coupling between particles and a viscoelastic fluid.