Simulation of Fluid Particle Cutting - Validation and Case Study

In this paper we present the comparison of experiments and numerical simulations for bubble cutting by a wire. The air bubble is surrounded by water. In the experimental setup an air bubble is injected on the bottom of a water column. When the bubble rises and contacts the wire, it is separated into two daughter bubbles. The flow is modeled by the incompressible Navier-Stokes equations. A meshfree method is used to simulate the bubble cutting. We have observed that the experimental and numerical results are in very good agreement. Moreover, we have further presented simulation results for liquid with higher viscosity. In this case the numerical results are close to previously published results.

1. Introduction. Fluid particle cutting plays an important role in gas-liquid and liquid-gas contactors. In gas-liquid contactors, the bubble size distribution, determining the mass transfer area, is influenced by the local hydrodynamics, but also by measuring probes such as needle probes [1,2,3] and mesh based conductivity probes [4]. The shape of the probes are mainly cylindrical, while the probe may be in flow direction but also in a rectangular angle to it. The rising bubbles approach the immersed object and starts to change its shape. Depending on the position of the bubble to the wire, the bubble will pass the object or be cutted in two fragments (daughter bubbles). Beside the unwanted cutting at probes, a wire mesh can be used to generate smaller bubbles and homogenise the flow structure. Furthermore, in liquid-gas contactors, phase separation is often a problem. Demisters are then frequently used to prevent a phase slip (entrainment) of fine dispersed phase droplets in the continuous product phase. A loss of the total solvent inventory within one year is reported causing costs and environmental hazards. Entrainment can cause a significant reduction in separation efficiency. Demisters are based on wire meshes, where the small droplets should accumulate. Using an optimal design, the small droplet seperation efficiencies can be up to 99.9%. Nevertheless, bigger droplets tend to break up in the rows of wires.
Hence, particle cutting is a frequently observed phenomena in various separation processes ranging from low viscosity to high viscosity of the continuous fluid. Nevertheless, it can be hardly investigated under operation conditions due to the complex insertion of optical probes into the apparatus or the complex mesh structure e.g. of the demister, but also the operation conditions as high pressure, high dispersed phase hold ups make an experimental investigation challenging.
In this study, we focus on the simulation of particle cutting at a single wire strengthened by experimental investigations to generate the basis for further numerical studies at complex geometries and fluid flow conditions such as demister simulations. For the simulation of bubble cutting, a meshfree approach is applied. It overcomes several drawbacks of classical computational fluid dynamics (CFD) meth-ods such as Finite Element Method (FEM) , Finite Volume Method (FVM). The main drawback of the classical methods (FEM, FVM) is the relatively expensive geometrical mesh/grid required to carry out the numerical computations. The computational costs to generate and maintain the grid becomes particularly high for complex geometries and when the grid moves in time, as in the case of fluid particles with a dynamic interface or in case where the interface between fluids changes in time.
For such problems meshfree methods are appropriate. Here, we use a meshfree method, based on the generalized finite difference method, called Finite Pointset Method (FPM). The two phase flow is modeled by using the continuous surface force (CSF) model [8]. Each phase is indicated by the color of the respective particles. When particles move, they carry all the information about the flow with them such as their color, density, velocity, etc. The colors, densities and viscosity values of all particles remain constant during the time evolution. The fluid-fluid interface is easily determined with the help of the color function [9]. In [13] an implementation of the CSF model within the FPM was presented to simulate surface-tension driven flows. We have further extended the method to simulate wetting phenomena [11].
2. Experimental Setup. Bubble cutting is investigated in a Plexiglas column filled with reversed osmosis water up to a level of 10 cm. The column has a width and depth of 46 mm. A syringe pump (PSD/3, Hamilton) is used to inject air of known volume at the bottom of the column. The injection diameter is 8 mm. The schematic setup is given in Figure 2.1. In a distance of 60 mm from the bottom, a wire is mounted in the middle of the column. The wire has a diameter of 3 mm. Two cameras (Imaging Solutions NX8-S2 and Os 8-S2) are mounted in an angle of 90 • to track the bubble motion over an image sequence, respectively over time. Both cameras are triggered and allow a synchronous detection at 4000 fps and a resolution of 1600x1200 px 2 . By tracking the bubble motion from two sides, it is possible to analyse the side movement and to detect the exact position of bubble contact with the wire. Also, the bubble deformation can be analyzed in two direction and therefore leads to more precise results compared to single camera setups.

Bubble Motion Analyses.
For bubble motion analyses, the tool box ImageJ (https://imagej.nih.gov/ij/) is used. The raw images (Fig. 3.1 a/b) are therefore binarized, followed by a watershed segmentation. The tracks of the single bubble and cutted particles are analysed using the Plugin Mtrack2 (http://imagej.net/MTrack2). Two particle tracks are tracked, one by each camera and reconstructed using Matlab software toolbox ( ( Fig. 3.1 c) These are the basis for three dimensional reconstruction of the bubble motion. The conversion of pixels to metric length is done by a afore performed calibration.
Matlab is used to reconstruct the bubble in a three dimensional domain. The images are converted to greyscale and further to binary images. Possible holes (white spots in a surrounded black bubble structure) are filled to get a better identification of the bubbles. To detect the bubble position, a distance transform is performed, followed by a watershed segmentation to separate the bubble from the pipe structure. Finally, the bubble size and shape is determined from each image. In a next step, the basic grey scale images are again converted to binary images, followed by a watershed algorithm [5]. Finally, the resulting structures are transformed into 3D space. Therefore, the detected structures are extruded into the third dimension resulting in overlapping structures. The overlapping structures represent the bubble and the wire and are visualized in Fig. 3.1d. Applying the assumption of an ellipsoidic structure for the bubble, results finally in Fig. 3.1e.

Mathematical Model .
We consider a one-fluid model for two immiscible fluids which are liquid and gas. We model the equation of motion of these fluids by the incompressible Navier-Stokes equations, which are given in the Lagrangian form where v is the fluid velocity, ρ is the fluid density, p is the pressure, τ is the stress tensor given by τ = 1 2 (∇v + (∇v) T ), g is the external force and F S is the surface tension force. The quantity F S is force density, which acts on the interface and its neighbor of the interface between gas and liquid. We compute the surface tension force CSF model of Brackbill et al ( [8]) and is given by where, σ is the surface tension coefficient, n is the unit normal vector at the interface and its neighbor, κ is the curvature and δ S is the surface delta function. We note that δ S is quite strong in the interface and its surroundings. We solve the equations (4.3) with initial and boundary conditions. 5. Numerical methods. We solve the equations (4.3) by a meshfree Lagrangian particle method. In this method, we first approximate a computational domain by discrete grid points. The grids points are divided into two parts as interior and boundary particles. The boundary particles approximate boundaries and we prescribe boundary conditions on them. The interior particles move with the fluid velocities. (e) Visualization based on the assumption of ellipsoidic structure.

Fig. 3.1: 3D bubble reconstruction
Particles may come very close to each other or can go far away from each other leading to very fine or very coarse approximations. This problem has to be tackled carefully due to stability reasons. To obtain a uniform distribution of particles in each time step one has to add or remove particles, if necessary. We refer to [10] for details of such a particle management.

5.1.
Computation of the quantities in surface tension force. For meshfree particle methods the interfaces between fluids are easily tracked by using flags on the particles. Initially, we assign different flags or color function c of particles representing the corresponding fluids. We define the color function c = 1 for fluid type 1 and c = 2 for fluid type 2. On the interface and its vicinity, the Shepard interpolation is applied for smoothing of the color functions using where m is the number of neighbor of arbitrary particle having position x, c i are the color values at neighboring particle i and w i is the weight as a function of distance from x to x i given by where α is a positive constant After smoothing the color function, we compute the unit normal vector Finally, we compute the curvature by The quantity δ s is approximated as Here, δ s is non-zero in the vicinity of the interface and vanishes far from it.

Numerical scheme.
We solve the Navier-Stokes equations (4.3) with the help of Chorin's projection method [7]. Here, the projection method is adopted in the Lagrangian meshfree particle method. Consider the discrete time levels t n = n dt, n = 0, 1, 2, . . . with time step dt. Let x n be the position of a particle at time level n.
In the Lagrangian particle scheme we compute the new particle positions at the time level (n + 1) by and then use Chorin's pressure projection scheme in new positions of particles. The pressure projection scheme is divided into two steps. The first step consists of computing the intermediate velocity v * with neglecting the pressure term Since we use the Lagrangian formulation, we do not need to handle the nonlinear convective term. The second step consists of computation of pressure and the velocity at time level (n + 1) by solving the equation where v n+1 should obey the continuity equation We observe from the equation (5.8) that the new pressure p n+1 is necessary in order to compute the new velocity. v n+1 . Now, we take the divergence of equation (5.8) on both sides and use of the continuity constraint (5.9), we obtain the pressure Poisson equation In order to derive the boundary condition for p we project the equation (5.8) on the outward unit normal vector n at the boundary Γ and then we obtain the Neumann boundary condition where v Γ is the value of v on Γ. Assuming v · n = 0 on Γ, we obtain We note that we have to approximate the spatial derivatives at each particle position as well as solve the second order elliptic problems for the velocities and the pressure. The spatial derivatives at each particle position are approximated from its neighboring clouds of particles based on the weighted least squares method. The weight is a function of a distance of a particle position to its neighbors. We observe that in Eq. 5.7 there is a discontinues coefficient µ inside the divergence operator since the viscosities of two liquid may have the ratio of up to 1 to 100. Similarly, the density ratio also has 1 to 1000, which can be seen also in Eq. 5.10. This discontinous coefficients have to be smoothed for stable computation. This is done using a similar procedure as for smoothing the color function. We denote the smoothed viscosity and density byμ andρ, respectively. We note that we smooth the density and viscosity while solving Eqs. 5.7 and 5.10, but keep them constant on each phase of particles during the entire computational time. If the density and viscosity has larger ratios, we may have to iterate the smoothing 2 or 3 times. Finally, Eq. 5.7 and Eq. 5.10 can be re-expressed as Note that, for constant density, the first term of Eq. 5.15 vanishes and we get the pressure Poisson equation. Far from the interface we haveμ = µ andρ = ρ. The momentum and pressure equations have the following general form 5.3. A meshfree particle method for general elliptic boundary value problems. In this subsection we describe a meshfree method for solving second order elliptic boundary value problems of type Eqs. 5.16 -5.17. The method will be described in a two-dimensional space. The extension of the method to three-dimensional space is straightforward. Let Ω ∈ R 2 be the computational domain. The domain Ω is approximated by particles of positions x i , i = 1, . . . , N , which are socalled numerical grid points. Consider a scaler function ψ(x) and let ψ i = ψ(x i ) be its discrete values at particle indices i = 1, . . . , N . We approximate the spatial derivatives of ψ(x) at an arbitrary position x ∈ {x i , i = 1, . . . , N }, from the values of its neighboring points. We introduce a weight function w = w(x i − x, h) with a compact support h. The value of h can be 2.5 to 3 times the initial spacing of particles such that the minimum number of neighbor is guaranteed in order to approximate the spatial derivatives. But it is user defined quantity. This weight function has two properties, first, it avoids the influence of the far particles and the second it reduce the unnecessary neighbors in the computational part. One can consider different weight function, in this paper we consider the Gaussian weight function defined in (5.2), where α is equal to 6.25. Let P (x, h) = {x j : j = 1, 2, . . . , m} be the set of m neighboring particles of x in a circle of radius h. We note that the point x is itself one of x j .
We consider Taylor expansions of ψ( ∂ 2 ψ ∂y 2 (y j − y) 2 + e j (5.18) for j = 1, . . . , m, where e j is the residual error. Let the coefficients of the Taylor expansion be denoted by a 1 = ψ(x, y), a 2 = ∂ψ ∂x , a 3 = ∂ψ ∂y , a 4 = ∂ 2 ψ ∂x 2 , a 5 = ∂ 2 ψ ∂x∂y , a 6 = ∂ 2 ψ ∂y 2 . We add the constraint that at particle position (x, y) the partial differential equation (5.16) should be satisfied. If the point (x, y) lies on the boundary, also the boundary condition (5.17) needs to be satisfied. Therefore, we add Eqs. 5.16 and 5.17 to the m equations (5.18). Equations 5.16 and 5.17 are re-expressed as where B = (B 1 , B 2 ) and n x , n y are the x, y components of the unit normal vector n on the boundary Γ. The coefficients a i , i = 1, . . . , 6 are the unknowns. We have six unknowns and m + 1 equations for the interior points and m + 2 unknowns for the Neumann boundary points. This means, we always need a minimum of six neighbors. In general, we have more than six neighbors, so the system is overdetermined and can be written in matrix form as  . , e m , e m+1 , e m+2 ) T and dx j = x j − x, dy j = y j − y. For the numerical implementation, we set n x = n y = 0 and ψ ΓN = 0 for the interior particles. For the Dirichlet boundary particles, we directly prescribe the boundary conditions, and for the Neumann boundary particles the matrix coefficients are given by Eq. 5.22. The unknowns a i are computed by minimizing a weighted error over the neighboring points. Thus, we have to minimize the following quadratic form The minimization of J with respect to a formally yields ( if M T W M is nonsingular) In Eq. 5.24 the vector (M T W )b is explicitly given by Equating the first components on both sides of Eq. 5.24, we get where Q 1 , Q 2 , . . . , Q 6 are the components of the first row of the matrix (M T W M ) −1 .
Rearranging the terms, we have We obtain the following sparse linear system of equations for the unknowns In matrix form we have where R is the right-hand side vector, Ψ is the unknown vector and L is the sparse matrix having non-zero entries only for neighboring particles. We solve the sparse system (5.29) by the Gauss-Seidel method. In each time iteration the initial values of ψ for time step n + 1 are taken as the values from previous time step n. While solving the equations for intermediate velocities and the pressure will require more iterations in the first few time steps. After a certain number of time steps, the velocities values and the pressure values at the old time step are close to those of new time step, so the number of iterations is dramatically reduced.
We stop the iteration process if for τ = 0, 1, 2, . . . and is a small positive constant and can be defined by the user.
6.1. Validation: low viscosity. In a first step, we validate the simulations with the experimental results of the single bubble cutting in reversed osmosis purified water. We apply the same bubble diameter from the experiments (6.5mm) and the wire diameter of 3mm in the simulation. The viscosity of the fluid (water) is µ l = 0.001P a.s and the interfacial tension between water and air is σ = 0.072N/m.. The density of water is ρ l = 998.2kg/m 3 and the density of air is approximated by ρ g = 1kg/m 3 and the dynamic viscosity of air is µ g = 2e −5 . . For the numerical simulations we consider a two-dimensional geometry of size 36mm × 63mm. The initial bubble position has the center at x = 18mm and y = 10mm and the wire has the center at x = 19.5mm and y = 45mm as shown in Fig. 6.1. The bubble is approximated by red particles, the liquid is approximated by blue particles. The white circular spot is the position of the wire. We have considered the total number of boundary particles, (including 4 walls and the wire) equal to 527 and the initial number of interior particles equal to 18310. The constant time step t = 5e −6 is considered. Here the horizontal distance between the initial center of bubble and the wire is d x = 1.5mm. In all four walls and the wire we have considered no-slip boundary conditions. Initially, the velocity and pressure are equal to zero. The gravitational force is g = (0, −9.81)m/s 2 .
The comparison between simulation and experiment is depicted in Fig. 6.2 and 6.3. We extracted a time sequence from the experiments and the corresponding simulations, starting at 0.19 seconds simulated time to 0.27 seconds. The temporal distance between each image is 0.02 s. The rising bubble approaches the wire and starts to deform. There is no direct contact during this phase between the bubble and the wire. Due to the non central approach to the wire, the bubble is cut in a smaller daughter bubble (right) and a larger bubble (left). The larger bubble has three times the diameter of the smaller bubble. The comparison of the cutting process gives a qualitatively good agreement between the experiment and the simulation. Also the shape and size of the mother and daughter bubbles are qualitatively very good agreement. A detailed comparison of the bubble path from experiment and simulation is shown in Fig. 6.4. The bubble position during first contact is important, which agrees well between experiment and simulation. Nevertheless, the path of the larger bubble in the simulation shows after the cutting a slightly different behaviour than in the experiment. In the experiment, the larger bubble moves inwards again, while the bubble in the simulation moves horizontally away from the wire, which may arise from a slight horizontally movement of the bubble in the experiment. The cutting of the bubble also depends on the bubble velocity, plotted in Fig. 6.5. The bubble accelerates in the simulation and finally reaches the same end velocity that is observed in the experiment. After the splitting into two daughter bubbles, the larger daughter bubble raises faster than the smaller one. The experimental results are governed by higher fluctuations especially for the smaller bubble, which results from very short temporal distances between the images and fluctuations by detecting the bubble interface. Neverthless, the average velocity for the smaller bubble fits with 0.12 m/s quiet well to the simulated result.
6.2. Case study: high viscosity. The computational domain is the same as in the previous section. However, the position of the wire is changed. The data has been taken from chapter 6 of [6]. The liquid has density ρ l = 1250kg/m 3 , dynamical viscosity µ l = 0.219P a.s. Similarly the gas density ρ g = 1kg/m 3 and the viscosity µ g = 2e −5 P a.s. The surface tension coefficient σ = 0.0658N/m. We consider a bubble of diameter 9.14mm with its initial center at (18mm, 9mm). We con- sider a wire (in 2D a circle) of diameter 3.1mm with different centers at y = 45mm and x = 18mm, 18.5mm, 19mm and 19.5mm. This means, we consider the initial distance d x between the center of the bubble and the center of the wire equal to d x = 0mm, 0.5mm, 1mm, 1.5mm. Fig. 6.6 shows the initial geometry with d x = 0mm. The initial number of particles and the time step are the same as in the low viscosity case. The initial and boundary conditions and the rest of other parameters are same as in the previous case.
In Figs. 6.7 -6.10 we have plotted the positions of the bubble and the wire for d x = 0mm, 0.5mm, 1mm and 1.5mm at time t = 0.288, 0.320, 0.352 and t = 0.384 seconds, respectively. For d x = 0 we observe the wire located in the middle of the bubble as expected. When we increased the distance d x from 0.5mm to 1.5mm, we observed that the left part of the bubble is increasing and the right part becomes smaller. We clearly observe that the daughter bubbles are symmetric for d x = 0mm in contrast to the other cases. We further observe a small layer between the wire and the bubble. After t = 0.352 seconds we observe the cutting of the bubble, see Figs. 6.9 and 6.10. Two daughter bubbles arise, a larger one on the left and a smaller one on the right side of the wire. The overall numerical results are comparable with the results presented in [6].
In Fig. 6.11 we have plotted the trajectories of the mother and bubble droplets. We observed that the mother droplet is cutted into two daughter bubble slightly below the wire, compare with Fig. 6.9. The trajectories are plotted up to time t = 0.5s. We see that when the size of the daughter bubble is increasing, it travels longer than the smaller bubbles. The reason is that the rising velocity of the larger bubble is larger than the smaller ones, see Fig. 6.12 for the velocities of mother and daughter bubbles. hydrodynamics of bubble cutting. For a first validation, we compared the solver to experimental data using the system air and water. A suffiently good agreement could be found in regard to bubble shape, bubble movement and cutting process itself. To study the effect of higher viscosity and bubble position, a case study was done. One observes that the inital position of the bubble to the wire has a high impact on the final daughter bubble size ratio. A centric approach of the bubble to the wire leads to a cutting of the bubble in two equally sized daughter bubbles. By increasing the inital distance to the wire, the daughter-bubble size ratio increases and the deviation between the velocities of daughter bubbles increases. Also, the movement of the bubbles directly behind the wire changes. While the bubbles split behind the wire at the centric approach, with increasing unsymmetry, the bubbles start to move inwards after the initial separation. In future, further studies with overlapping wires are planned.