# Numerical Simulation of Micro-Bubbles Dispersion by Surface Waves
## 3 Numerical Method
First, an intermediate velocity, \(U_{i}^{*}\), is computed using the velocity fields at the preceding time steps [26]:  
\[U_{i}^{*}=U_{i}(t_{k})+\bigg{(}\frac{3}{2}F_{i}(t_{k})-\frac{1}{2}F_{i}(t_{k-1 })\bigg{)}\Delta t, \tag{16}\]  
where the flux, \(F_{i}\), is evaluated as,  
\[F_{i}=-\frac{\partial(U_{i}U_{j})}{\partial x_{j}}+\frac{1}{\text{Re}}\frac{ \partial^{2}U_{i}}{\partial x_{j}\partial x_{j}}+\sum_{n=1}^{N_{b}}S_{i}^{n}. \tag{17}\]  
Further, the new pressure, \(P(t_{k+1})\), is computed by solving its respective Poisson equation in the form,  
\[\frac{\partial^{2}P(t_{k+1})}{\partial x_{j}\partial x_{j}}=\frac{1}{\Delta t }\frac{\partial U_{j}^{*}}{\partial x_{j}}. \tag{18}\]  
Equation (18) is solved by iterations by performing, at each iteration step (_j_) the FFT in the horizontal directions and Gaussian elimination in the vertical direction. The iteration procedure stops when the condition \(\big{|}P_{j+1}-P_{j}\big{|}/P_{j}<0.1\%\) is satisfied. Usually this condition is met after \(j\) = 3-5 iterations. The new velocity at \(k+1\) time step satisfying the incompressibility condition (2) is then computed as:  
\[U_{i}(t_{k+1})=U_{i}^{*}-\frac{\partial P(t_{k})}{\partial x_{j}}\Delta t. \tag{19}\]  
At the upper boundary, \(\widetilde{\eta}=\eta=0\), the no-slip (Dirichlet) condition for the velocity is prescribed:  
\[ \begin{array}{l}U_{x}(\xi,y,0)=U_{d}-cka\ \cos k(x(\xi,0)-ct),\\ U_{y}(\xi,y,0)=0,\\ U_{z}(\xi,y,0)=-cka\ \sin k(x(\xi,0)-ct),\end{array}  \tag{20}\]  
where the surface drift velocity, \(U_{d}\), is expressed as [27]:  
\[U_{d}=c+cka\ \cos k(x(\xi,0)-ct)-\big{[}(c+cka\ \cos k(x(\xi,0)-ct)\big{)}^{2}-q(2 c-q)\big{]}^{1/2}. \tag{21}\]  
Two different cases are considered. In one case, parameter \(q\) is put to zero, \(U_{d}=0\), and the wind stress effects are not taken into account. In this case, the water velocity at the boundary coincides with orbital velocities of the fluid particles in the surface wave. In another case, \(q=0.05c\) is prescribed, and thus \(U_{d}\) is finite so that the wind-stress effects upon the water surface are accounted for. Periodic conditions for all fields are prescribed at the side boundaries, and the no-slip (Dirichlet) condition for the water velocity is prescribed at the bottom boundary (\(\widetilde{\eta}=-1\)). The bubble equation of motion is solved by employing the Adams-Bashforth method:  
\[V_{i}^{n}(t_{k+1})=V_{i}^{n}(t_{k})+\bigg{(}\frac{3}{2}F_{i}^{n}(t_{k})-\frac {1}{2}F_{i}^{n}(t_{k-1})\bigg{)}\Delta t \tag{22}\]  
\[F_{i}^{n}=3\frac{DU_{i}^{n}}{Dt}+\frac{f(\text{Re}_{b}^{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{23}\]  
In Equation (23), the surrounding fluid velocity, its acceleration, and the ambient-flow vorticity at the location of each bubble are obtained by a Hermitian 4 th-order accuracy interpolation procedure. The dimensionless bubble response time, \(\tau_{n}\), (or the Stokes number) is expressed as:  
\[\tau_{n}=\frac{d_{n}^{2}U_{0}}{36\
u L_{0}}. \tag{24}\]  
The bubble coordinate equation is advanced in time by employing the Adams method as:  
\[r_{i}^{n}(t_{k+1})=r_{i}^{n}(t_{k})+0.5(V_{i}^{n}(t_{k+1})+V_{i}^{n}(t_{k})) \Delta t \tag{25}\]  
The water velocity is initialized as a random, zero-mean field with an amplitude of 0.1% (normalized by the velocity scale, \(U_{0}\)). After initiation, a transient occurs during which the velocity field adjusts to boundary conditions, and a statistically stationary flow state is reached (at dimensionless time \(t=100\)). The results (not shown) indicate that the the mean, wave velocity field is established comparatively early (\(\sim\)O(10) wave periods) and the transient to the stationary flow state is mainly related to the development of a near-surface, wave-induced turbulent layer (of \(\sim\)O(100) duration, cf., e.g., [17]).