# Numerical Simulation of Micro-Bubbles Dispersion by Surface Waves
## 2 Governing Equations  
Figure 1 outlines the schematic of the problem under consideration. A domain with sizes \(L_{x}=6\lambda\), \(L_{y}=4\lambda\), \(L_{z}=\lambda\), with periodic side boundaries, a solid bottom and a waved upper boundary, \(z_{s}(x,t)\), is considered, and the Cartesian coordinates are employed, \(0\leq x\leq L_{x}\); \(-L_{y}/2\leq y\leq L_{y}/2\); \(-L_{z}\leq z\leq z_{s}\). The motion of the fluid (water) is driven by the upper boundary where a progressive, two-dimensional (2D) stationary wave of amplitude \(a\), celerity \(c\) and length \(\lambda\) propagating in the positive \(x\)-direction is prescribed. Figure 1: Schematic of the problem: (\(x,y,z\)) are the Cartesian coordinates; \(a\), \(\lambda\), \(c\) are the surface wave amplitude, length, and celerity; \(g\) is the acceleration due to gravity; \(L_{x},\ L_{y},\ L_{z}\) are the domain sizes. Symbols (open circles) show injected bubbles (not to scale). A Eulerian-Lagrangian approach is adopted where the Navier-Stokes equations of the water motion are solved in a Eulerian frame, and the bubbles are tracked simultaneously by solving their respective equations of motion in a Lagrangian frame. The Navier-Stokes equations for the carrier fluid are written in the dimensionless form [24]:  
\[\frac{\partial U_{i}}{\partial t}+\frac{\partial(U_{i}U_{j})}{\partial x_{j}}=- \frac{\partial P}{\partial x_{j}}+\frac{1}{\mathrm{Re}}\frac{\partial^{2}U_{i }}{\partial x_{j}\partial x_{j}}+\sum_{n=1}^{N_{b}}S_{i}^{n}, \tag{1}\]  
where \(x_{i}\equiv(x,y,z)\), \(U_{j}=(U_{x},U_{y},U_{z})\) are water velocity components, \(P\) is the pressure, and \(S_{i}^{n}\) is a momentum source term (cf. Equation (8) below) contributed by the \(n\)-th bubble, and \(n=1,\,\ldots,\,N_{b}\), the latter being the total, constant, number of tracked bubbles. Variables in Equation (1) are normalized with velocity scale, \(U_{0}\), set equal the surface wave phase speed, \(c\), and length scale, \(L_{0}\), equal to the wave length, \(\lambda\). The pressure is normalized with \(\rho U_{0}^{2}\) where \(\rho\) is the water density (\(\approx\)1 g/cm\({}^{3}\)). The carrier flow Reynolds number, Re, measures the ratio of the inertial vs. viscous effects and is defined as:  
\[\mathrm{Re}=\frac{U_{0}L_{0}}{\upsilon}, \tag{2}\]  
where \(\upsilon\) is the water kinematic viscosity (\(\approx\)0.01 cm\({}^{2}\)/s). Equation (1) is supplemented by the incompressibility condition:  
\[\frac{\partial U_{j}}{\partial x_{j}}=0, \tag{3}\]  
which implicitly defines the pressure field, \(P\) (cf. Equation (18) below). The bubbles equations of motion are written as [15]:  
\[\frac{dr_{i}^{n}}{dt}=V_{i}^{n}, \tag{4}\]  
\[\frac{dV_{i}^{n}}{dt}=3\frac{DU_{i}^{n}}{Dt}+\frac{f(\mathrm{Re}_{i}^{n})}{ \tau_{n}}(U_{i}^{n}-V_{i}^{n})+\frac{1}{2}\;\varepsilon_{ijk}\Big{(}U_{j}^{n} -V_{j}^{n}\Big{)}\omega_{k}^{n}+2g\delta_{iz}. \tag{5}\]  
In Equations (4) and (5), \(r_{i}^{n}\), \(V_{i}^{n}(i=x,y,z)\) are the \(n\)-th bubble coordinate and velocity components, \(U_{i}^{n}\) is the fluid velocity at the bubble location, \(d/dt=\partial/\partial t+V_{j}^{n}\partial_{j}\) and \(D/Dt=\partial/\partial t+U_{j}^{n}\partial_{j}\) are the material derivatives along the bubble trajectories and the surrounding fluid Lagrangian paths, respectively; \(g\) is the dimensionless gravitational acceleration, and \(\omega_{i}^{n}=\varepsilon_{ijk}\partial_{i}U_{k}^{n}\) is the surrounding fluid vorticity; \(\delta_{ij}\) and \(\varepsilon_{ijk}\) are the Kronecker and Levi-Civita tensors. The forces acting on the bubble (in the order as in the right hand side of Equation (5)) are the fluid acceleration, viscous drag, lift, and buoyancy.