Direct numerical simulation of a bubble suspension in small amplitude oscillatory shear flow

Bubble suspensions can be found in many different fields and studying their rheology is crucial in order to improve manufacturing processes. When bubbles are added to a liquid, the magnitude of the viscosity changes and the behavior of the material is modified, giving it viscoelastic properties. For the purpose of this work, the suspended bubbles are considered to be monodisperse. It is assumed that Brownian motion and inertia can be neglected and that the fluid of the matrix is Newtonian and incompressible. The suspension is subject to an oscillatory strain while remaining in the linear regime. The resulting equations are solved in 3D with direct numerical simulation using a finite element discretization. Results of an ordered and random distribution of bubbles of volume fractions up to 40% are presented. The presence of bubbles has an opposite effect on the rheology of the suspension depending on the applied frequency. When the frequency is low, bubbles act as rigid fillers giving a rise to viscosity. On the contrary, when the frequency is high, the strain rate is being accommodated by the gaseous phase. Hence, bubbles deform, leading to a decrease of the viscosity.


Introduction
Dispersions of bubbles in a liquid are very common in many different fields and understanding the rheology of such materials is essential in order to improve manufacturing processes. Bubble suspensions are encountered in nature as magmas (Manga and Loewenberg 2001) and in food industry where the gas phase can affect the texture, spreadability and provide a more uniform distribution of taste (Thakur et al. 2003). In polymeric foams, the change of bubble volume fraction due to bubble growth greatly affects material properties and has been studied extensively by many authors (Aloku and Yuan 2010;Everitt et al. 2003). Usually, they consist of a continuous liquid phase with bubble volume fractions up to 50% and in foams even more.
Adding bubbles to a liquid not only changes the magnitude of the viscosity but also modifies the behavior of the material giving it viscoelastic properties (Llewellin et al. 2002). Both Rust and Manga (2002) and Stein and Spera (2002) identify two different regimes of such a material in flow. When the applied strain rate is small, the bubbles act as rigid fillers increasing viscosity, while at large strain rates, where the flow line distortion is smaller due to bubble deformation, bubbles provide free slip surfaces and thus decreasing viscosity. Llewellin et al. (2002) highlighted the importance of distinguishing between steady and unsteady flow. If shear conditions remain constant for a time much longer than the relaxation time of the bubble, the flow is assumed to be steady.
There is a lot of work done in the past century on the shear viscosity of suspensions of particles and emulsions. The relative shear viscosity η r is the viscosity of the suspension normalized by the viscosity of the matrix μ and is a function of the volume fraction of the dispersed phase, φ. Taylor (1932) generalized Einstein's equation for dilute suspensions of rigid particles (Einstein 1906) to drops where η d is the viscosity of the dispersed phase. If we assume that the viscosity of the dispersed phase is small, i.e., η d μ, for a dilute suspension of bubbles we get η r = 1 + φ.
(2) Llewellin et al. (2002) provided an extensive review of previous models developed so far together with their own model for suspensions of bubbles under steady and oscillatory shear conditions.

Problem definition
For the purpose of this work, the suspended bubbles are considered to be monodisperse of radius R and they are not subject to Brownian motion. The suspension is subject to an oscillatory strain γ (t) = γ 0 sin(ωt) while keeping γ 0 small, so that we remain in the linear regime. The fluid of the matrix will be Newtonian. The initial geometry of suspended spherical bubbles inside a cubic domain is shown in Fig. 1. The bubble interfaces are described by the surfaces S i and periodic boundaries of the box ∂ j for j = 1, ..., 6. The periodic domain will be considered as a representative volume element of the bulk and the effect of the volume fraction on the relative viscosity will be investigated.

Governing equations
To describe the flow dynamics, we assume that inertia can be neglected and that the fluid is incompressible. As a result the momentum and the mass balance reduce to where u is the velocity vector and σ the stress tensor. For incompressible Newtonian fluids, the stress tensor σ is written as: where p is the pressure, I is the identity tensor, μ is the viscosity of the matrix fluid, and D = (∇u + (∇u) T )/2 is the symmetric part of the velocity gradient tensor. The gas inside the bubbles is assumed to be incompressible. For that reason, we impose a constant volume with S i the ith bubble boundary on which we define n as the outwardly directed unit normal vector. The Nth bubble constraint is fulfilled automatically because of the incompressibility of the matrix fluid. Thus inside bubble N, we have zero pressure and for the remaining bubbles the pressure is such that the volume remains constant.

Boundary conditions
To close the balance equations, the interface of the bubbles is assumed to be a material interface and triperiodic boundary conditions are imposed for the velocity and traction. The interfacial boundary conditions for a sharp interface can be expressed as where the subscript m denotes the matrix fluid, is the surface tension coefficient and κ is the curvature. The periodicity of the domain can be formulated as whereγ is the imposed shear rate and is time dependent, L is the length of the cubic computational domain, and additionally we have imposed where t = σ · m is the traction vector with m the outwardly directed unit normal vector on the outer boundaries of the periodic domain. It is required to indicate a reference value for the velocity in a single point in the periodic domain to remove the rigid body motion of the domain. The periodic boundary conditions and the bubble volume constraints are imposed using Lagrange multipliers.

Interface tracking
The boundary of the sharp interface between two phases can be described explicitly by a surface mesh. We can track the interface by following the motion of the surface mesh. An interface in 3D can be described by a moving curvilinear coordinate system given by whereξ = (ξ 1 , ξ 2 ) are the curvilinear coordinates andx is the function that maps the coordinatesξ onto the spatial coordinates x. For a material interface, the velocity of the interfaceẋ must be such that: where u is the material velocity at the interface and n is the normal vector of the surface. The interface velocity is given bẏ Note, thatẋ is equal the velocity of the material in the normal direction only. The following equation for the motion of an interface has been implemented: where u c is the velocity of the center of volume of the bubble. The main purpose for introducing the field u c is to make the movement of the nodes on the interface Galilean frame independent. For implementation, Eq. 13 is rewritten as (Villone et al. 2014) x + (u − u c ) · (I − nn) = u.
By substituting the expression for the surface identity tensor I − nn = g α g α , where g α , α = 1, 2 the dual base vectors and g α = ∂x/∂ξ α , α = 1, 2 are the covariant base vectors. For the term g α g α , we use the Einstein notational convention for simplifying the expression of summation. We can rewrite Eq. 14 tȯ where ∇ s = g α ∂ ∂ξ α is the surface gradient operator.

Techniques to obtain rheology
The stress is one of the most important rheological properties and we would like to know the stress generated by imposing a certain strain history to the bulk. However, the stress response will vary with position of the bulk and in time. Thus, the "bulk stress" only exists in an average sense (Batchelor 1970). To calculate the average stress, we assume that the triperiodic domain is a representative volume element which is small compared to macroscopic variation of the suspension but still large enough compared to the distances between bubbles. We describe the macroscopic problem by a representative fully periodic domain.
The box over which we calculate the stresses is the periodic domain. By integrating the stress over the box and dividing by its volume, we can obtain the volume-average stress. The volume-average stress is given by Batchelor's formula (Batchelor 1970) One common way to probe the rheological properties of a fluid is the use of oscillatory measurements (Macosko 1994). In oscillatory measurements, a sinusoidal strain deformation of magnitude γ 0 is applied at a frequency ω.
For a linear response, the resulting stress response remains sinusoidal with a phase shift and amplitude that depends on the material and the frequency. Thus, the stress response can be split in two parts, the in-phase response and the out-ofphase response, as follows: where G (ω) and G (ω) are the elastic and viscous modulus respectively. To ensure that the response is linear for a given strain amplitude γ 0 , we check if by doubling the strain amplitude, the stress response doubles as well If we write the applied oscillatory strain deformation in the form of a complex function γ (t) = γ 0 e iωt , then the two responses would be determined respectively by the real and imaginary part of the stress response, written as where G * (ω) is the complex shear modulus From the complex shear modulus, the complex viscosity can be obtained from

Numerical description Weak form
The weak form can be obtained by multiplying Eqs. 3, 4, and 15 with test functions v, q and w: Find u, p, x such that where the inner products are defined as Using partial integration and Gauss' theorem, we obtain the following weak form: Find u, q, and x such that using appropriate spaces for u, p, x, v, q, and w. Furthermore, similar to Glowinski et al. (1999), we use the discrete inner product on the bubble interface · , · S which will be defined in "Discretization" section. In the weak form of the momentum balance, Eq. 29, D v is defined as

Discretization
The Eqs. 29-31 are discretized using the finite element method employing a mesh of quadratic tetrahedra. The interface mesh consists of quadratic triangles which are conforming with the volume mesh. Quadratic interpolation for the velocity and linear interpolation for the pressure (Taylor-Hood) is used. We solve the resulting equations coupled so that our timestep is not limited by the characteristic mesh size time limit. In a decoupled approach, the timestep is limited by the capillary time which is t = xμ/ where x is the element size, μ is the viscosity of the matrix and the surface tension (Courant et al. 1967). A second-order time integration scheme (BDF) is used for the time derivative of the position of the bubble interface. Equations 29-31 form a non-linear system of equations which needs to be linearized so that it can be solved with a Newton-like iterative scheme.
It is important to note here that since the position x as an unknown exists only on the bubble interface the method is not fully Newton-Raphson because the domain is not perturbed. To linearize velocity and position, u and x can be written as follows where subscript n indicates the current iteration and δu, δx are the differences of velocity and position after one iteration respectively. For implementation the surface tension term after manipulation is rewritten taking into account the change in position and area as The discrete inner product · , · S is analogous to imposing a constraint using a collocation method and is given by where n coll is the number of collocation points. The collocation points coincide with the nodes on the bubble interface. Furthermore, since the pressure term is linear, we choose to solve for p. Substituting in Eqs. 29-31 and dropping the subscript n leads to To take the change of area into account, Eq. 31 is multiplied by ∇ s · δx and added to Eq. 38, while keeping only linear terms. The converged solution is obtained when the maximum difference of r x,k = x k − x k−1 , r u,k = u k − u k−1 and r p,k = p k − p k−1 of two consecutive iterations is less than r max < b 2 where = 10 −8 and b 2 the 2-norm of the right hand side of the linear system of equations. By numerical experimentation, we found that multiplying the matrix terms containing the unknown δx in Eq. 34 by −1 greatly increases convergence. Since the resulting system of equations can be large, the necessary physical memory for using direct solvers can grow fast. For that reason in most of the simulations presented in this paper, we use iterative solvers, such as GMRES, from the SPARSKIT library (Saad 2001) in combination with an ILUT preconditioner. Since the connectivity of the nodes does not change, the structure of the system matrix remains the same and thus it is possible to re-use the ILUT preconditioner for many timesteps or even for the whole simulation. To increase the convergence rate by reducing the fill-in of the preconditioner, the mesh nodes are renumber using MeTiS (Karypis and Kumar 1998). Nevertheless, selecting the proper parameters for the iterative solver can be challenging. In those cases, we used a direct solver from the HSL library (HSL 2013).

Mesh movement
Since our problem contains moving boundaries, i.e., bubble interfaces and boundaries of the domain, it is essential to define the motion of the mesh by using a smooth displacement field. We do this by solving a Laplace equation in a periodic domain, similar to Hu et al. (2001). Firstly, we have to solve the Laplace problem to obtain the position of the periodic domain in the next time step where z 1 , are the nodal mesh displacements and k e is a constant diffusion coefficient. The periodicity of the domain can be formulated as where γ is the strain difference between two consecutive time steps The mesh displacement can be then used to find the new mesh coordinates as Then, the bubble interfaces are moved during the Newton iterations according to the solution δx of the interface motion. The rest of the periodical boundaries remain constant during the iterations. This leads to the following set of equations Once more, we obtain the new mesh coordinates using the mesh displacements

Mesh generation
For a given bubble distribution, a triperiodic mesh has to be generated. For numerical reasons, we want to do this such that there are no intersections between the bubbles and the outer surfaces. However, due to the bubbles being randomly distributed, it is difficult (or even impossible) to find straight planes that do not intersect any of the bubbles. We have therefore designed an algorithm that generates surfaces that curve around the bubbles. The algorithm is based on an approach used in Jaensson et al. (2015), where a convectiondiffusion problem is solved for the nodes on the outer curves of a 2D mesh to obtain smooth boundaries that do not intersect particles. We start the mesh generation by finding the corner points of the mesh. Note, that by choosing one corner point, the remaining seven corner points are fixed due to the periodicity of the problem. The corner points are typically chosen far away from any bubble boundary, to avoid problems with meshing the volume. Next, coarse surface meshes are generated of three mutually perpendicular straight surfaces (the remaining three surfaces are periodic copies). A number of nodes on these surfaces are likely to be located inside the bubbles, therefore we solve the following convection-diffusion equation on these surfaces where z are the nodal coordinates, ∇ s is the surface gradient operator, α is a diffusion coefficient, u a is an artificial velocity field, and m is the normal to the surface. The diffusion coefficient α is chosen such that the resulting surface is smooth. The velocity field u s points away from the center of the bubbles, with a magnitude that quickly decays to zero outside the bubbles. This ensures that nodes are "pushedout" of the bubbles, but not into adjacent bubbles. Similar to Jaensson et al. (2015), a tanh-function is employed to determine the magnitude of u a . Equation 48 is integrated in time using an explicit Euler scheme. During the time integration, the planes are iteratively refined in regions where bubbles are close, ensuring that the element size on the surfaces matches the required element size of the volume mesh. An example of the three surfaces after running the algorithm is shown in Fig. 2. Finally, the periodic surfaces are used to generate a triperiodic volume mesh using Gmsh (Geuzaine and Remacle 2009).

Convergence test
Convergence is tested for a problem of a single bubble in a periodic domain with length L = 1, which is equivalent to regularly stacked distribution of bubbles. For the constant parameters, we choose φ = 0.3, μ = 1, = 0.1, γ 0 = 0.04. Three different meshes are used with varying number of elements on the equator of the bubble, n eq = 10, 20 and 40. In Fig. 3, the error e r of the viscosity is plotted for the cases of n eq = 10, 20, which is relative to the case of n eq = 40. We can see that the error quickly converges for n eq = 20 giving a maximum relative error of less than 0.5%. Choosing now the number of elements on the equator to be n eq = 20, we perform simulations with the same constant parameters. This time, we use different number of steps n t in each period. In Fig. 4, we plot the error for n t = 20, 40 relative to the case of n t = 80. This clearly indicates that Fig. 3 Relative error e r of the viscosity using different meshes with the number of elements on the equator of the bubble being n eq = 10, 20. As a reference value, we use the result obtained from the case of n eq = 40 Fig. 4 Relative error e r of the viscosity using different number of steps in each period of the oscillations. The error is plotted for the cases of n t = 20, 40 relative to the one obtained with n t = 80 our numerical method is converging and from now on we will use n eq = 20 and n t = 40, which is a good trade off between accuracy and numerical efficiency.
In the case of multiple bubbles and for higher volume fractions, one should take into account the thin layers that are present between bubble interfaces. In order to accurately describe the flow dynamics between bubbles, it is essential to perform mesh refinement. Figure 5 shows amplitude of the viscous stress tensor in the thin layer with use of 1, 3, and 6 elements between bubbles surfaces. It can be seen that three elements are enough to accurately describe the solution in the thin layer between the bubble interface. Hence, in the remainder of this paper, we will use at least three elements between bubbles.

Results
To study the effect of suspending bubbles in a Newtonian fluid, many simulations were performed using our model with regularly distributed, as well as randomly distributed bubbles. The problem is described by the following dimensionless parameters, the relative viscosity as the ratio of the magnitude of the complex viscosity η * (ω) over the fluid viscosity μ and the frequency scaled with the relaxation time λω. The relaxation time is defined as which means that λω represents the ratio of viscous forces versus surface tension forces. The results are presented for Starting from the dilute regime, it is possible to validate our results, by comparing with analytical results that already exist in the literature (Llewellin et al. 2002). We then present results from simulations that were performed with a volume fraction well above the dilute limit up to φ = 0.4. This can be used to study the effect of the volume fraction on the rheology of the suspension.

Validating against an analytical model
As mentioned earlier in "Techniques to obtain rheology" section, one convenient way to obtain rheological properties is to apply an oscillating strain γ and measure the resultant stress response. In Llewellin et al. (2002), the linearized Frankel and Acrivos model was used to describe a system which is subjected to forced oscillations. They give the following expression for the dynamic viscosity η and loss viscosity η where ω is the angular frequency and α 1 , β 1 , β 2 are related to physical properties of the bubble suspension We are going to use this model, which is the exact asymptotic result under the assumptions of the dilute limit and small bubble deformation, to verify our numerical results. For the constant parameters, we choose φ = 0.005, μ = 1, = 0.1, γ 0 = 0.04 and we perform simulations for a range of frequencies ω scaled by relaxation time λ from 10 −2 up to 10 2 . To validate our method, the result of our simulation for the dynamic η and elastic η components of the complex viscosity are compared with the expression given by the analytical model. The results are presented in Fig. 6 and both the components show a good match for the whole range of the scaled frequency λω. The small difference between simulations and analytical results can be explained by the fact that the analytical model has been derived for the case of a single bubble in an infinitely large domain, while in our numerical model there is a finite distance between the bubble and the periodic boundaries. In Fig. 7, the relative viscosity η r = η 2 + η 2 /μ over the scaled frequency λω is compared with the analytical results from the linearized Frankel and Acrivos model. As we can see, there is an excellent agreement between the two under the same assumptions of volume fraction in the dilute limit and small bubble deformations.

Regularly stacked bubbles
Having verified our numerical model, we continue with more simulations and investigate cases where the assumption of the dilute limit does not apply any more. We will begin by using only one sphere in a periodic domain as shown in Fig. 8, which represents a regularly stacked distribution of bubbles. All of the following simulations for Fig. 6 In the left figure, the dynamic component η is compared for the case of one bubble of φ = 0.005 between analytical results and simulations. While on the right the elastic component η is plotted the case of regularly stacked bubbles were performed with constant parameters μ = 1, = 0.1, γ 0 = 0.04 for a wide range of volume fractions and frequencies. In Fig. 9, we can see the effect of the increasing volume fraction to the relative viscosity η r . When the volume fraction increases the relative viscosity diverges more and more from what would be a constant line at η r = 1 for the case of a Newtonian fluid without any suspended bubbles. As mentioned in Llewellin et al. (2002), that behavior can be explained as follows. For low λω, the change in the bubble deformation is much smaller in comparison with the applied strain rate and the bubbles can be seen as rigid fillers, thus the suspension medium takes most of the strain. When φ increases, flow lines are increasingly distorted, which causes viscosity to increase with φ. On the other hand, for very high λω, most of the strain is taken by the bubbles and since their viscosity is negligible an increase of the volume fraction would lead to an decrease of the relative viscosity. To quantify that, it Fig. 7 Comparison between simulations and analytical results of the relative viscosity over a range of scaled frequency of regularly stacked bubbles of φ = 0.005 is possible to calculate the strain rate accommodated by the fluid similarly to the bulk stress calculation using Eq. 56 whereγ f is the average of the maximum strain rates for all oscillations. The ratio ofγ f /γ , as depicted in Fig. 10, is in agreement with the claim previously stated.
Since the matrix viscosity μ is known as a direct input of the model it is interesting to compare the results after subtracting it, as shown in Fig. 11. Considering that the analytical model contains only the linear terms of φ, it is possible to check what is the effect of the higher order terms (Fig. 12). Using a function for η that contains the higher order terms φ 2 and φ 3 and fitting the parameters a, b and c which are all functions of frequency λω, it is possible to create the master curves as shown in Fig. 13. These master curves indicate the importance of the higher order terms. The fit was performed for 0.005 ≤ φ ≤ 0.4 and the square of the correlation between the actual values and the predicted ones, r 2 is greater than 0.99. As expected the parameter a approximates the values obtained from the analytical expressions.

Random distribution of bubbles
Ordered structures as the one used in the previous section are rarely encountered in nature. For that reason, studying the change in rheology due the increase of volume fraction Fig. 10 The ratio of the average of the maximum strain rates for all oscillation over the applied strain rateγ f /γ , versus the scaled frequency λω Fig. 11 The dynamic viscosity η minus the matrix viscosity μ as a function of volume fraction φ, for different scaled frequencies λω of a random distribution of bubbles is essential. The representative volume element (Fig. 14) consists of a periodic mesh where the periodic boundaries curve around the bubbles as explained in "Mesh generation" section. To avoid extreme cases where bubbles nearly overlap, the random distribution is constraint such that the minimum distance between bubble interfaces can not be less than 0.1R. To complete our simulations in a feasible time, we had to constrain the number of bubbles to 6. For all the simulations of the multiple bubble problem, the following constant parameters where chosen. The viscosity of the matrix and surface tension are the same as in the case of regularly stacked bubbles (μ = 1, = 0.1), where a lower strain amplitude is preferred (γ 0 = 0.001) to ensure numerical Fig. 12 The difference of the dynamic viscosity η between numerical and analytical results as a function of volume fraction φ, for different scaled frequencies λω Fig. 13 Parameters a, b, c of Eq. 57 are plotted over the whole range of scaled frequencies stability. That is especially more imperative for the cases of high volume fractions. The range of volume fractions is again 0.1 ≤ φ ≤ 0.4 and we perform simulations for a wide range of frequencies. To obtain the correct rheological response when using representative volume elements, it is important to average over multiple realizations. Since our simulations are rather time consuming, it would be impossible to verify that for the whole range of volume fractions and frequencies. We instead calculated the standard deviation for volume fraction φ = 0.3 and ω = 50 which was 0.0112. Hence, we are confident that our results represent accurately the true rheological response. In Fig. 15,   Fig. 14 Geometry of a random distribution of bubbles with a boundary fitted mesh Fig. 15 The difference of the relative viscosity η r between regularly and randomly distributed bubbles as a function of scaled frequency λω, for different volume fractions φ the two different cases, of regularly and randomly distributed bubbles are plotted together. As one might expect, the relative viscosity of the randomly distributed bubbles show an increase of viscosity in the low frequency plateau compared with the case of an ordered distribution. This can be explained if we keep in mind the explanation given by Llewellin et al. (2002) and illustrated in "Regularly stacked bubbles" section. Since the strain is taken by the viscous phase, a random distribution of bubbles would cause bigger disturbance in the flow lines than an ordered structure and thus would lead to an increase of viscosity. On the other hand, in the high-frequency regime, this increase is smaller since most of the strain is accommodated by the inviscid phase. Note that the vertical axis in Fig. 15 has a logarithmic scale.

Conclusion
We have presented 3D direct numerical simulations of bubbles suspended in a Newtonian fluid subject to oscillatory strain. The influence on the rheology of the suspension due to the increase of the volume fraction was investigated for the cases of regularly stacked bubbles and a random distribution of bubbles in a fully periodic domain. The resulting stress, as an output of our simulation, was used to calculate the viscosity of the bubbly fluid. Using the linearized Frankel and Acrivos model which is derived for the dilute limit and small strain amplitudes, we were able to validate our results. Comparing the analytical results, where only linear terms of volume fraction are included, with our simulations we were able to show the significance of higher order terms. The increase in volume fraction has an opposite effect on viscosity in the high and low frequency regimes. When the frequency is small, the bubbles act as rigid fillers, thus increasing the viscosity. In contrast, in the high-frequency regime the bubble deformation rate compared to the bulk strain rate is much larger, hence, most of the strain is taken by the inviscid phase leading to a decrease in viscosity.
We ran preliminary tests to investigate the effect of the number of bubbles in the domain. The first results show that the effect is minor, especially for low volume fractions. Further optimization of the numerical technique will allow us to study more bubbles and higher volume fractions. For this work, we only consider bubbles as the suspended phase. However, the model can be easily extended to simulate droplets, and the results could be compared to the theory of Oldroyd (1953). Studying large amplitude oscillatory shear in bubble suspensions is part of ongoing research. To achieve this, it is essential to apply a remeshing technique to avoid highly deformed elements. Furthermore, a multilevel adaptive refinement method is necessary to describe flow dynamics between two bubble interfaces that approach each other. In this paper, we studied only the case of a Newtonian matrix. It should be noted that, it would be relatively easy to extend our numerical model to the case of viscoelastic matrices.