Governing equations

Translational and angular momentum balances

In addition to the Navier-Stokes equations

\[\begin{split}\der{u_i}{x_i} &= 0, \\ \der{u_i}{t} + u_j \der{u_i}{x_j} &= - \der{p}{x_i} + \frac{1}{Re} \der{}{x_j} \der{u_i}{x_j} + a_i,\end{split}\]

which govern the behaviour of the liquid, we consider the Newton-Euler equations to describe the motions of rigid particles (density: \(\rho_p\), volume: \(V_p\), surface: \(S_p\) or \(\partial V_p\)):

\[\rho_p V_p \oder{U_i}{t} = - \int_{\partial V_p} a_i d S_p + \oder{}{t} \int_{V_p} u_i d V_p + \left( 1 - \rho_p \right) \frac{1}{Fr^2} + F_i,\]

\[\oder{}{t} \left( I_m \Omega_i \right) = \int_{\partial V_p} \epsilon_{ijk} r_j a_k d S_p + \oder{}{t} \int_{V_p} \epsilon_{ijk} r_j u_k d V_p + T_i.\]

Note that they are already normalised by the reference scales, and \(Fr\) is the Froude number defined as

\[Fr \equiv \frac{U}{\sqrt{gL}}.\]

The first equation describes the translational velocity \(U_i\) of the gravity center \(X_i\), while the second equation tells how the angular velocity \(\Omega_i\) is evolved in time. The first terms in the right-hand-sides (including \(a_i\)) denote the force and torque caused by the fluid-structure interactions, while the second terms are responsible for the motion of fictitious fluid inside objects (see Breugem, J. Comput. Phys. (231), 2012). The last terms \(F_i\) and \(T_i\) are the external force and torque, which take into account the gravitational acceleration and collisions between the other objects in this project.

Since I am interested in ellipses, it is worthwhile to introduce the following basic parameters and their relations:

\[\begin{split}\begin{cases} \text{major axis} & a, \\ \text{minor axis} & b, \\ \text{volume} \, V_p & \pi a b, \\ \text{mass} \, m_p & \rho_p \pi a b, \\ \text{moment of inertia} \, I_m & \frac{1}{4} \pi \left( a^2 + b^2 \right) \end{cases}\end{split}\]

Note

Since I limit my focus to two-dimensional objects for now, the volume and surface integrals reduce to the surface and line integrals, respectively. Also the moment of inertia inside the temporal derivative is a constant for each object and thus can be taken out, i.e., the left-hand-side of the equation of angular velocity leads

\[\oder{}{t} \left( I_m \Omega_i \right) = I_m \oder{\Omega_i}{t}.\]

Collision

The lubrication force coming from the hydrodynamic interaction between particles (or a particle and a wall) tends to be underestimated especially when the resolution between the objects is insufficient, and as a result particles can penetrate (too much) to each other. In order to avoid this, this library employs a collision model between objects, which is based on a simple spring model:

\[m \oder{U_i}{t} + k X_i = 0,\]

in which a spring whose spring constant is \(k\) is assumed between objects and a corresponding force whose direction is parallel to the normal and magnitude is proportional to the penetration depth pushes away the objects to each other. This force is only activated when a penetration is observed and try to resolve the penetration gently (in \(T\)), which is the so-called soft-sphere collision model (Tsuji, KONA Powder and Particle J. (11), 1993).

A general solution of the above equation with boundary conditions

\[X_i \left( t = 0 \right) = X_i \left( t = T \right) = 0\]

leads

\[X_i = A \sin \sqrt{\frac{k}{m}} t,\]

where

\[k = m \frac{\pi^2}{T^2}\]

and thus the repellent force as a function of the penetration depth in the normal direction \(x_i\) is given by

\[f_i = k x_i.\]

Since the particles are rigid, \(T = 0\) in theory, which is impractical within the current numerical scope (so-called event-driven collision model deals with \(T = 0\), which is not easy to integrate in time simultaneously as the fluid behaviour).

As a remedy, in this project, I consider a timescale

\[T = \sqrt{\frac{\rho r^3}{\sigma}} \propto We^{\frac{1}{2}} r^{\frac{3}{2}},\]

which is analogous to the capillary timescale for deformable objects whose \(\sigma\) is the surface tension coefficient (and \(We\) is a pseudo Weber number). Note that \(r\) is a particle-based length scale (e.g., harmonic average of the colliding particle radii).

In theory, \(\sigma \rightarrow \infty\) (or equivalently \(T \rightarrow 0\)) since the objects are rigid. Giving larger \(T\) stabilises the integration process because the repellent force becomes smaller, which is especially useful when the suspension volume fraction is relatively large (e.g., \(\gt 40 \%\)). However, the outcome might be unphysical because of the large penetration depths, and thus a proper timescale should be determined by the user eventually.

The above spring model does not include a dumper, indicating that the collision is perfectly elastic (restitution coefficient is \(1\)) and no dissipation is involved. Obviously this is not true in most cases and also elastic, dumping, and sliding motions in the tangential direction should be included (see Costa et al., Phys. Rev. E (92), 2015 for extensive analyses).

In this project, however, I neglect all these aspects just for simplicity, and include only springs in the normal directions to obtain plausible results.

An extension of the above concept to elliptic objects is not straightforward, which will be discussed later.