Drag of an elliptic intruder
in a two-dimensional granular environment

As shown in Fig. 2, the movement parallel to the major axis is unstable. Even if the initial direction of the major axis of the intruder $\theta$ is parallel to the drag direction, this direction tends to be perpendicular to the drag direction ($\theta=\pi/2$ or $-\pi/2$ depending on the initial configuration). This can be understood as following: If there exists a small perturbation $\Delta\theta$, the torque works on the intruder with the same sign of $\Delta\theta$, which means that the small perturbation enhances the rotation of the intruder for $\theta\simeq 0$. Therefore, from now on, we present results when we pull the intruder with the initial angle $\theta=\pi/2$.

Figure 3 shows the relationship between the drag force $F_{\rm drag}$ and the steady speed $V$ When we compare the results of the elliptic case and the cylindrical case with $D=15d$, the resistance becomes larger for the elliptic case. This can be understood as follows: If the surrounding particle collides with the intruder at a certain $y_{0}$, the corresponding $x$ coordinates become $R_{2}\sqrt{1-y_{0}^{2}/R_{1}^{2}}$ and $\sqrt{R^{2}-y_{0}^{2}}$ for the elliptic and cylindrical cases, respectively, where $R=D/2$ is the radius of the cylinder. If we put the angle between the $x$ direction and the normal vector at the colliding point as $\beta$ the detailed derivation is given in the next subsection, and the work from the particle to the intruder is proportional to $\cos^{2}\beta$ [15, 33]. Then, this angle becomes smaller for the elliptic case, which also means that $\cos\beta$ becomes larger. Summing up all the contributions around the surface of the intruder, the drag resistance becomes larger for the elliptic case.

Let us derive the relationship between the dynamic part of the drag force $F_{\rm dyn}$ and the steady speed $V$. We consider the coordinate that the intruder is set at the origin with the angle $\theta=\pi/2$. Moreover, we consider a situation where the surrounding particles move towards the negative $x$ direction with the steady speed $V$ and collide with the intruder. First, we consider a point $(R_{2}^{\prime}\cos\vartheta,R_{1}^{\prime}\sin\vartheta)$ ($-\pi/2\leq\vartheta\leq\pi/2$), where we have introduced $R_{1}^{\prime}\equiv R_{1}+d/2$ and $R_{2}^{\prime}\equiv R_{2}+d/2$. This point indicates the center of the colliding particle whose diameter is $d$. Here, this system can be regarded as the collision between the mass point and the elliptic intruder with $R_{1}^{\prime}$ and $R_{2}^{\prime}$. Then, the unit normal vector at this point becomes

with

It should be noted that the angle of the vector $\bm{n}$ from the $x$ axis is not $\vartheta$ but $\beta$. The momentum change in the $x$ direction due to this collision is given by $\Delta p_{x}\simeq\{[(2/3)+e_{n}-e_{t}/3]\cos^{2}\beta+(1+e_{t})/3\}m_{\rm r}V$ with the reduced mass $m_{\rm r}\equiv Mm/(M+m)$ [33]. In addition, the volume of the collision cylinder corresponding to the parameter $\vartheta$ becomes $\mathcal{V}_{\rm coll}=V(R_{1}^{\prime 2}\cos^{2}\vartheta+R_{2}^{\prime 2}\sin^{2}\vartheta)^{1/2}\cos\beta$. Integrating the momentum change in the $x$ direction over the surface of the (modified) intruder, we get

with

where we have introduced the eccentricity $\epsilon\equiv[1-(2R_{2}+d)^{2}/(2R_{1}+d)^{2}]^{1/2}$. In the above derivation, we regard all particles around the intruder as monodisperse with the diameter $d$ and the mass $m$. This simplification does not affect the results. We note that we assume that the momentum arm is fixed as a constant independent of the collision point. Our collision model (7) recovers that for a cylindrical intruder [7] if we take the limit $R_{2}\to R_{1}$. This estimation works well as shown in Fig. 3 at least when the dimensionless speed $V^{*}$ is sufficiently smaller than unity.

In this subsection, we compare streamlines obtained from the simulation with those of perfect fluid around an elliptic intruder. Now, let us derive the latter from Ref. [34]. In this case, it is well known that a complex velocity function $W(z)=\phi(z)+i\psi(z)$ is a powerful tool to describe a two-dimensional flow [35, 36], where $\phi(z)$ and $\psi(z)$ are the velocity potential and the stream function, respectively, with $z=x+iy$, and $i$ is an imaginary number. First, we consider a flow around a cylinder. Assuming a uniform flow with the speed $-V$ at infinity in the $x$-direction, the stream function can be described by [35]

When the $\zeta$-plane is isometrically mapped to the $z$-plate by

the flow around a cylinder (9) is transformed into the flow around an elliptical cylinder whose minor axis coincides with the flow ($x$-)direction. Actually, in the $\zeta$-plane, the circle $\zeta=ce^{i\theta}$ of radius $c$ located at the origin is projected to

in the $z$-plane (see Fig. 4). We note that $c\geq a$ should be satisfied. Now, both $c$ and $a$ should be

respectively. From the complex velocity function (9) with the mapping (10), we can obtain the complex velocity function which describes the flow around an elliptic intruder as

Then, the streamlines are characterized by $\psi(z)=\Im W(z)={\rm const}$ [35].

Figure 5 exhibits the comparison of the streamlines from the simulation with those of the perfect fluid given by Eq. (13), where we also plot those from the viscous fluid which is explained in the next paragraph. The theoretical streamlines well reproduce the simulation results in front of the intruder. Now, we note that the streamlines from the simulations are inconsistent with the assumption which is used when we derive the collision model, where the trajectories of a particle colliding with the intruder are straightforward. Actually, this assumption is not right because particles collide with each other before they collide with the intruder. However, this assumption works well as a zeroth-order approximation to derive the drag law because the collision model works well as compared to the simulation results.

We also compare the results with streamlines of the viscous fluid. It is well known that the stream function $\Psi$ of the viscous fluid satisfies the biharmonic equation $\Delta\Delta\Psi=0$. Following the procedure reported in Ref. [37], we can write the solution of the biharmonic equation as

where the expressions of the coefficients $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$, and $\mathcal{D}$ are listed in Table 1 (see Appendix B for details). As shown in Fig. 5, these streamlines cannot capture the behavior near the surface of the intruder, where the streamlines are bent further from those of the perfect fluid.

In this subsection, we focus on the stress fields around the intruder when the drag force $F_{\rm drag}$ is smaller than the threshold $F_{\rm Y}$, i.e., the intruder stops due to the balance between the drag force and those from the surrounding particles.

To characterize the stress, we adopt the similar method used in our previous paper [15]. Once we adopt the theory of elasticity, the stress component is calculated from the simulation as

with $\alpha,\beta=x,y,z$, where we have introduced the coarse graining method with the width $2d$ [38]. Now, it is convenient to consider the stress components in the elliptic coordinates $(\varrho,\vartheta)$ rather than the Cartesian coordinates $(x,y)$, either of which is converted from the other as

with $c=\sqrt{R_{1}^{2}-R_{2}^{2}}$.

Figure 6 shows typical profiles of the components of the stress on the surface of the intruder.

To characterize the stress profiles, let us introduce the Fourier components of the stress on the surface of the intruder as

for $n=0,1,2,\cdots$, where we have introduced the normalized stress $\widetilde{\sigma}_{\varrho\vartheta}=\sigma_{\varrho\vartheta}/(F_{\rm drag}/\pi d)$. We also note that $b_{0}^{(\varrho\vartheta)}=0$ is always satisfied by its definition. Figure 7 shows the behaviors of the Fourier components obtained from the simulation, where we have only shown $a_{n}^{(\varrho\varrho)}$, $b_{n}^{(\varrho\vartheta)}$, and $a_{n}^{(\vartheta\vartheta)}$ due to their symmetric properties. The coefficients for $n\geq 4$ become smaller than those for $n\leq 3$, which suggests that it is sufficient to consider terms up to the third-order even when we consider the stress around the intruder.

Once we adopt this assumption, the expressions of the stress components in the elliptic coordinates are listed in Table 2, where we have introduced $\vartheta^{\prime}\equiv\vartheta+\pi/2$ (see Appendix C for details).

In Table 2, we have ten unknown coefficients; $E$, $F$, $G$, $H$, $A_{2}$, $B_{2}$, $C_{1}$, $C_{3}$, $D_{1}$, and $D_{3}$. Now, let determine these coefficients from the information of them on the surface of the intruder [15]. From Eq. (17) and Table 2, the following equation is satisfied on the surface of the intruder:

where the vectors $\bm{X}$ and $\widetilde{\bm{\Sigma}}_{\rm surf}$ are defined by

respectively, and the each component of the matrix $\mathcal{M}$ (the size of $\mathcal{M}$ is $10\times 10$) is listed in Table 3.

Now, we can numerically solve Eq. (18) as

Using this solution, we can evaluate the stress fields. Figure 8 shows the comparison of $\sigma_{xy}$ obtained from the simulation and that from the expressions listed in Table 2 with Eqs. (45).