Immersed boundary simulation of drop stability

The stability of the quasi-two-dimensional droplet flow is of great importance in microfluidic devices. We check the drop’s stability in the square box using the immersed boundary and lattice Boltzmann methods. We implement two-dimensional equations within the immersed boundary approach in the Palabos programming platform. We check the influence of the boundaries on the drop movement. We estimate fluctuations in the quantities while applying different initial conditions of the linear and angular velocities. We found that the level of fluctuations depends on the symmetrical displacement of drop at the initial state. The effect is connected with the hydrodynamic interaction of drop with the walls.


Introduction
The drop flows are subject to intensive studies both in the laboratory setups [1] and in vivo [2]. There are many intriguing peculiarities of the flow in a microfluidic water-in-oil droplet generator [3]. The experiment shows the generation of acoustic waves in a quasi-two-dimensional flow of the droplet chain and instability due to the long-range hydrodynamic interaction. The value of parameters corresponds to the low Reynolds number hydrodynamics. It is worth checking the lattice Boltzmann method's applicability to the chain's simulation and comparing results with the detailed experiment.
In the paper, we report the preliminary results. We build a computational model based on the immersed boundary method together with the lattice Boltzmann method. As the first step, we simulate the drop stability comparing the symmetrical and asymmetrical position of the drop in the channel. We found that the hydrodynamic interaction of the drop with the walls cannot be negligible for long enough simulations. We check how the fluctuations in the observables depend upon the drop's initial linear and angular velocities. Thus, we can estimate how long we can simulate the chain, avoiding the unnecessary effects of the walls' finite distance.
We consider different settings: eight combinations of zero/non-zero initial translational velocity along the x-axis, zero/non-zero initial angular velocity, and symmetrical/asymmetrical placement of the drop in the domain. As a part of a study, we provide the results for different lattice sizes.
The paper is organized as follows. In section 2, we describe the simulation model and simulation methods. In section 3, we present results of the drop stability varying distance to the walls and viscosity. Section 4 briefly summarize the results.

Model and Methods
We can consider the two-dimensional flow since in the quasi-two-dimensional experiment [3] the only effect of the friction of the drop with the flow and celling is the renormalization of the drag force. We consider the drop with the size R much smaller than the size of the box L, and it is in line with the experiment. We simulate drop-fluid interaction with the immersed boundary (IB) method [4].
We are using the lattice Boltzmann method (LBM) [5,6] with incompressible Bhatnagar-Gross-Krook (BGK) single relaxation time collision model with a D2Q9 velocity set. The simulations were done using IB-LBM method [7] utilizing the open-source program platform Palabos [8], in which we modify the immersed boundary method from 3D case to 2D case.

Immersed boundary method.
The immersed boundary (IB) method is used to simulate blood flow, bubble dynamics, and many other phenomena. It applies to the moving boundaries of the drop in the flow. Peskin [9] initially proposed it for blood flow simulation in 1972. Since that time method has been used for particle suspension simulations: coupled with difference scheme for the Navier-Stokes equation [9,10,11] and lattice Boltzmann equations [12,4,7,13,14,15].
The main idea of the IB method is the representation of the drop as a set of Lagrangian boundary points, which generally do not coincide with Eulerian lattice nodes. The combined IB-LBM approach is an iterative method. For each iteration, the value of body force for Lagrangian points is calculated. For that purpose, velocities from the Eulerian grid are interpolated into Lagrangian boundary points. Then body force, which acts from drop to fluid for the Eulerian point, is computed. The computed force is used for velocity computation in lattice points. Finally, velocities interpolated into the drop boundary point, and correction of the drop's body force is calculated. We set number of IB iterations to 4. The details of IB method usage are described in Refs. [4,7,16].

Drop motion.
The drop is represented by a set of evenly distributed points on the boundary of a rigid disk. The drop motion is defined by classical mechanics laws I is the total force on the drop, computed by IB method, I is the intertial tensor and for the rigid disk it is equal to , Ω is the angular velocity, T = − r × g is total torque acting on drop.
The boundary velocity u b calculated with where r b is coordinates of boundary point and X drop center of mass.
In two-dimensions, the torque has only one non-zero component and Ω = {0, 0, ω} T , so the Eq.(2) reduces to We integrate system of Eqs.(1-4) using the Euler method.

Simulation setup.
The grid size is N ∆x in both box directions, with N chosen from the set of pairs 125, 126, 250, 251, and 500, 501 -to investigate the symmetrical and asymmetrical drop displacement against the walls. The drop radius R = 10∆x, with ∆x = 1 is the lattice spacing of LBM. The drop density ρ is equal to the fluid density ρ = ρ f luid = 1. The drop boundary is defined as the set of 75 nodes, evenly distributed over the circle of radius R.
We use in the paper the LBM lattice units [5,6]. We set the values of LBM parametersthe relaxation parameter τ = 23/10 and time step ∆t = 1. Therefore, the viscosity ν = 0.6, the speed sound c s = 1/ √ 3, and Reynolds number is equal to unity. To investigate the influence of symmetry on drop motion in the box, we consider eight possible cases shown in Table 1. Those cases are combinations of zero/non-zero initial translational velocity along the x-axis, zero/non-zero initial angular velocity, and symmetrical/asymmetrical placement of the drop in the domain. The computational domain is a square box with a bounceback boundary condition on walls. v 0 is the initial translational velocity along the x-axis in lattice units, ω 0 is the initial angular velocity in lattice units.
Surrounding liquid is in the rest at the initial moment of simulations, and translational and rotational movement of the drop generates the waves in the liquid [17,18]. The velocities decay because of the finite viscosity of the liquid.

Simulation results
We check how much the simulation depends on the drop placement's asymmetry in the box and on the distance between the drop and the wall. Asymmetry is realized using the odd and even number of lattice spacings with the center of drop always placed at the LBM lattice vertex. This way, we can estimate how long we can simulate the drop chain's movement and investigate the chain instability avoiding the finiteness of the box influence on each drop. It is necessary to exclude the chain's instability triggered by the setup asymmetry while simulating the drop chain movement using a combined IB-LBM approach. In Figures 1-3, the solid lines correspond to the asymmetrical case and dashed lines to the symmetrical case. Colors mark the following sets of lattice sizes. The blue color marks set N = (125, 126), the salmon color marks set N = (250, 251), and the magenta color marks set N = (500, 501).

Asymmetry influence.
The initial values of the translational and angular velocities are set to 0 or 0.06, see Table 1.
The drop center is always placed at integer values (N/2, N/2). We distinguish two setups. For even box size N , the drop is placed precisely at the computational domain center, and it is the symmetrical setup. For the odd box size N , the drop is shifted by (−0.5∆x, −0.5∆x), and it is the asymmetrical setup. The distance from the two sides of the drop to the walls is different in the asymmetrical setup. Parameters for computational experiments are presented in Table 1.
The details of the asymmetry influence on the velocity decay are: • Cases 0 and 1. Both initial velocities are zero, and there is no motion at all, hence no difference for symmetrical and asymmetrical drop placement. The drop remains unmoving during the simulation. • Cases 2 and 3. The angular velocity displays no difference for symmetrical and asymmetrical cases as can be seen from Figure 1b. The translational velocity shown in Figure 1a behaves differently for symmetrical and asymmetrical cases, decaying slowly for the asymmetrical case. • Cases 4 and 5. Symmetry of the drop displacement in the box does not reflect in the translational velocity behavior as shown in Figure 2a, while the angular velocity behavior is quite different. In the symmetrical case, it is zero with the accuracy of simulation, and it is small but finite in the asymmetrical case, see Figure 2b. The level of angular velocity excitation is small and may be attributed to the IB-LBM method's accuracy.  for the sizes 250 and 500 behaves quite slowly -see Figure 1b. The translational velocity shown in Figure 1a is sensitive to both the asymmetry and the box size. • Cases 4 and 5. Velocity behavior is shown in Figure 2a and 2b. The box size does not substantially influence the angular velocity behavior. The translational velocity behaves differently with the box size. There are visible oscillations of velocity with the frequency connected with the wave propagation from the drop to the wall and back. The period of oscillation is proportional to the ratio of the wave's traveled distance to the sound speed. • Cases 6 and 7. Velocity behavior is shown in Figures 3a and 3b. In this case, the translational velocity behavior is similar to those for cases 4 and 5. Moreover, the angular velocity behavior is similar to cases 2 and 3.

Conclusions
We perform simulations of the drop in the square box, analyzing the boundary conditions' influence. We found the asymmetry can be essential for the nonzero linear and angular velocity. We choose the same initial values for angular velocity and translational velocity, which corresponds to the different Reynolds numbers. The circular velocity associated with the angular velocity is R = 10 times larger than the translational one, which leads to the large Reynolds number. Indeed the effects associated with the initial angular velocities are more pronounced. We present the preliminary analysis defining the safe time length of simulations for the drop chain instability analysis.