The evolution of a three-dimensional microbubble at a corner in a Maxwell fluid

Bubbles often appear in non-Newtonian liquids from nature, engineering to biomedical applications, but their study has been under research compared to their Newtonian counterpart. Here, we extend the axisymmetric modeling of Lind and Phillips to three-dimensional modeling. The approach is based on the boundary integral method coupled with the Maxwell constitutive equation. The flow is assumed to have moderate to high Reynolds numbers and, thus, is irrotational in the bulk domain. The viscoelastic effects are incorporated approximately in the normal stress balance at the bubble surface. The numerical model has excellent agreement with the corresponding Rayleigh – Plesset equation for spherical bubbles in a non-Newtonian liquid. Computations are carried out for a bubble near a corner at various angles. The numerical results agree very well with the experiments for bubbles in a Newtonian fluid in a corner. As the Deborah number increases, the amplitude and period of the bubble oscillation increase, the bubble migration to the corner enhances, and the bubble jet is broader, flatter, and inclined more to the further boundary. This implies an improvement to surface cleaning of all surrounding boundaries for ultrasonic cavitation cleaning and results in greater administration of noninvasive therapy and drug delivery.


I. INTRODUCTION
Bubbles in non-Newtonian fluids are widely encountered in various domains such as decompression sickness, 1 volcanic eruption, 2 glass manufacture, materials, petrochemicals, metallurgy, and different dispersed systems. 3Certain additions of polymers have been reported to affect cavitation damage [4][5][6] and noise. 72][13][14] The non-Newtonian properties of the fluid are important due to the rheological properties of the blood and other bodily fluids. 15,16he behavior of viscoelastic fluids is characterized by performing experiments in the laboratory using a rheometer.Examples of simple flows that are studied include simple shear flow, extensional flow, and oscillatory shear flow.The latter is used to determine the relaxation times of the fluid, a key measure of the viscoelastic nature of the fluid.There are two classes of fluids that often display Maxwell-like behavior in the normal measuring range (10 À2 À 10 2 Hz) with the expected shapes of the curves for the storage and loss moduli and a single relaxation time.(Note that this frequency range refers to the experiments conducted by the rheometer and not the oscillatory frequency of the bubble.)9][20] Fluids in the latter class are sometimes called "living polymers" because if they break under large stresses, they can reconstitute under conditions of rest or low stress.
Extensive research has been developed for bubbles in a Newtonian fluid. 21,22One of the computational models that has been used widely in bubble simulations is the boundary integral method (BIM) based on a potential flow theory for inertially dominant flows. 23,246][27][28][29] Weak viscous effects are approximated in this model using the viscous pressure correction, which is determined by the conservation of energy at the interface. 30,31ubble dynamics were also simulated based on the Navier-Stokes equations using the finite volume method (FVM) [32][33][34][35][36] or finite element method. 37,38These are mainly for axisymmetric profiles.
Contrary to Newtonian bubble dynamics, investigations on bubbles in non-Newtonian fluids have been sparse.A rising bubble in a non-Newtonian fluid was observed experimentally by Wagner et al. 39 This was numerically simulated by Lind and Phillips, 40 using the BIM along with Maxwell's constitutive equation.Lind and Phillips 41 also simulated bubbles near a free surface.Liu et al. 42 computed equally spaced bubbles rising in a viscoelastic fluid.
Recent developments have focused on a three-dimensional bubble rising in a viscoelastic fluid, 43,44 with particular emphasis on predictions of the velocity discontinuity at critical volume. 45,46However, the domain still remains radially symmetric, and three-dimensional effects are scarcely analyzed.
In Secs.II and III, we extend the axisymmetric viscoelastic BIM of Lind and Phillips 40,41,47 to three-dimensional modeling.The viscoelastic liquid is modeled using the Maxwell constitutive equation.The flow is assumed to have moderate to high Reynolds numbers and, thus, is irrotational in the bulk domain.The viscoelastic effects are approximately included in the normal stress balance at the bubble surface.In Sec.IV, the numerical model is validated by comparing with the corresponding Rayleigh-Plesset equation for spherical bubbles in a non-Newtonian liquid and experiments for bubbles in Newtonian fluid in a corner.In Sec.V, computations are performed for a bubble near a corner at various corner angles using Green's function for Laplace's equation in the corner domain.Viscoelastic effects are studied in terms of the Deborah number De.In Sec.VI, some conclusions are made with reference to the bubble radius amplitude, oscillation period, migration to the corner, Kelvin impulse, and bubble jet profile.

II. PHYSICAL AND MATHEMATICAL MODEL
We consider a bubble in a viscoelastic fluid located near a corner formed by two rigid flat boundaries.We assume that the associated liquid flow has moderate to high Reynolds numbers, and thus, viscoelastic effects are negligible in the bulk.The viscoelastic effects are not negligible in a thin viscous layer at the bubble surface, which is approximated only through the normal stress balance.Viscous fluid dynamics can be described approximately by potential flows when the vorticity is small or is confined to a narrow layer near the boundary. 48It is particularly useful for a gas-liquid two-phase flow with an interface.
Assuming that the flow is irrotational and incompressible in the bulk, the velocity, u Ã ¼ r Ã / Ã satisfies Laplace's equation in the fluid given by r Ã2 / Ã ¼ 0: ( The kinematic boundary conditions on the two rigid boundaries, S W , the bubble surface S B , and the far field are given, respectively, by where r Ã p is a material point at the bubble surface and r Ã is a field point.
At any instant in time, the fluid pressure at the bubble surface is related to the bubble pressure and viscoelastic effects, through the normal stress balance 40,47 where P Ã L is the liquid pressure on the bubble surface, r is the surface tension coefficient, j Ã is the curvature of the bubble surface, P Ã B is the internal pressure of bubble gases, and T Ã nn represents the component of deviatoric stress on the bubble surface.Here, P Ã vc is proportional to the normal stress T Ã nn induced by the irrotational velocity P Ã vc ¼ ÀCT Ã nn where the constant C, given by Manmi and Q. Wang, 30 is Maxwell's constitutive equation is used to model the non-Newtonian properties of the fluid.The equation yields a general irrotational equation of motion and provides no contribution to stress in the bulk 49 and so the viscoelastic effects only appear in the normal stress balance condition.However, the model comes with the limitation of only being applicable for small deformations. 50Thus, a "material" Maxwell model is used.The material Maxwell model means that the material time derivative collapses to a time derivative in the particle reference frame and can, therefore, be easily calculated. 40he Maxwell model can be applicable based on the following reasons: the bubble is approximately spherical during most of its lifetime due to surface tension.It may become non-spherical during a very short period at the end of collapse when the inertial effects are dominant and the viscoelastic effects are negligible.In addition, as for the moderate Reynolds number Newtonian case, jet formation is completely suppressed in the case of a bubble in a viscoelastic fluid near a rigid boundary. 51Jet suppression has also been seen experimentally in viscoelastic fluids. 52,53axwell's constitutive equation is given by 54 where k 1 is the relaxation time of the fluid.
For general constitutive equations, we cannot find a function w that satisfies r Á T ¼ rw for general irrotational flows. 49In these situations, we have r Â r Á T 6 ¼ 0 even though r Â u ¼ 0, for velocity u.We cannot always find a function w satisfying the above equation so that not all fluids will satisfy a Bernoulli equation.However, under the assumptions described above, the admissibility condition can be satisfied in an approximate sense for moderate to large Reynolds numbers since r Á T becomes small compared to inertial terms in the momentum equation in the bulk of the flow. 47Hence, a Bernoulli equation is admissible for inviscid Newtonian, viscous Newtonian, and linear Maxwell fluids with appropriate Reynolds numbers, where w ¼ 0. Therefore, under such a condition, the potential flow theory can accurately provide a three-dimensional description for a bubble in a Maxwell fluid.
Assuming that the gas-bubble is adiabatic, the pressure inside the bubble is given by where P g0 is the initial partial pressure of the gas, V 0 is the initial volume of the bubble, and f is the heat capacity ratio.Using Bernoulli's equation for the velocity potential and equation ( 5), we obtain where q is the density of the surrounding fluid, g is the gravitational constant, and P 0 is the ambient pressure of the fluid surrounding the bubble.
From now on the non-spherical bubble model is converted into a nondimensional problem with the length scale being the initial bubble radius R 0 and the pressure scale being P 0 À P v .The stand-off distances, c N and c F , are the normalized distances between the bubble center and the near and far wall, respectively.Hence, we obtain the dimensionless model for the non-spherical bubble where d ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gR m q=ðP 1 À P v Þ p is the buoyancy parameter, e ¼ P g0 = ðP 0 À P v Þ is a measure of the initial bubble gas pressure, We ¼ R 0 ðP 0 À P v Þ=r is the Weber number measuring the effect of surface tension, 56 and Re ¼ R 0 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðP 0 À P v Þq p =l L is the Reynolds number. 30The Deborah number, De, of the fluid is given by 57

De
We note that as k 1 approaches zero, that is, the relaxation time of the fluid becomes instantaneous, De also approaches zero, and the Maxwell model reverts back to the Newtonian model.

III. NUMERICAL MODEL
The initial mesh is first approximated using an icosahedron with 20 equal-sized triangles and 12 nodes situated on a spherical bubble surface.The mesh is then improved by dividing each triangle into four smaller triangles with new nodes added at the midpoint of each line segment and projected onto the surface of a sphere.On each triangular element, the field point, r, velocity potential, /, and normal velocity, @/=@n, are linearly interpolated. 25 high-quality surface mesh of the bubble surface is maintained by implementing a hybrid of the Lagrangian method and elastic mesh technique. 58The bubble surface and potential distribution were interpolated using a polynomial scheme coupled with the moving least squares method for calculating the surface curvature and tangential velocity on the surface. 59,60he boundary integral method (BIM) can be implemented with Green's second identity given by cðr; tÞ/ðr; tÞ ¼ ð @X @/ðq; tÞ @n Gðr; qÞ À /ðq; tÞ @Gðr; qÞ @n dS; (14)   where cðr; tÞ is the solid angle, n is the unit outward facing normal at the surface, and @X represents every boundary in the domain, X.For a corner angle of p=k, for some natural number, k, Green's function is given by where q 0 is the source point and q j are the image points of q 0 for j ¼ 1; 2; …; 2k À 1.The images were provided by Kucera and Blake 61 and Tagawa and Peters. 62A proof is provided by Wang et al. 63 To calculate the normal stress, T nn , we need to calculate @ 2 /=@n 2 .This is calculated using the following equation: Since / satisfies Laplace's equation ( 1) and Green's second identity (14), so do / x ; / y and / z can be calculated using the BIM.Once @ 2 /=@n 2 is found, the normal stress, T nn , can be obtained using a backward Euler approximation 57 for Maxwell's constitutive equation (12), given by De T nn ðxðtÞ; tÞ À T nn ðxðt À DtÞ; where the variable time step is chosen as for some constant D/ found experimentally.
The BIM is grid-free in the flow domain making it computationally efficient; thus, the approach has been used extensively in many bubble dynamics simulations.At each time step, the bubble surface and potential distribution on the surface are known.The tangential velocity on the bubble surface is computed from the gradient of the potential using a polynomial interpolation combined with a least squares method.An advanced Linear Algebra Package (LAPACK) is used to solve the linear system ( 14) with a seven-point Gaussian quadrature being used to integrate the off diagonal elements of the matrix of coefficients.The weak singularities that occur in Eq. ( 15) are dealt with via a change to polar coordinates. 64,65The remaining singularities are eliminated by adopting a 4p rule. 60Solving the boundary integral equation ( 14) yields the normal velocity on the surface.The bubble surface and potential distribution are updated on the bubble surface using Eqs.(10) and (18), with a second-order Runge-Kutta scheme.The stress tensor component, T nn , is then calculated using a backward Euler method.A variable time step (18) is chosen for accuracy and efficiency, with the maximum change of potential at each time step being dependent on constant D/. 23,66A hybrid approach for the Lagrangian method and elastic mesh technique is implemented to maintain a high quality surface mesh throughout the simulations. 58

IV. VALIDATION A. Comparison with the Rayleigh-Plesset equation
We first compare the predictions of our numerical model with the solution of the Rayleigh-Plesset equation (RPE) for a spherical bubble with radius R(t) in a non-Newtonian fluid given by The microbubble has an initial radius of R Ã 0 ¼ 4:5 lm, and a total number of 5494 nodes of the surface mesh are used.
As shown in Fig. 1, the BIM predictions agree excellently with the RPE for the bubble radius history over two oscillation cycles for a Deborah number of De ¼ 2.0.A similar agreement has been found for De ¼ 0.0, 0.2, and 20.0.
Figure 2 displays that an increase in the Deborah number results in a significant increase in the maximum radius of the bubble over the first oscillation from 0.8203 to 0.8616 and an even greater increase over the second oscillation from 0.7395 to 0.8615.This is due to more elastic energy being stored and released in the viscoelastic fluid for a larger Deborah number and less energy being lost due to viscosity.This is apparent from Eq. ( 17), because the model approaches the inviscid case as De increases.Similarly, an increase in De also significantly increases the oscillation period of the bubble.

B. Comparison with experiments
To evaluate the BIM model, we compare the computational results with experimental data, 55 for a laser beam generated bubble with maximum radius R m ¼ 0:88 mm at a right angled corner of a ¼ p=2.The stand-off distances from the near and far walls are c N ¼ 0:88 and c F ¼ 1:08, respectively.The dimensionless parameters of the fluid are given by e ¼ 100 and j ¼ 1:4, where Re ¼ 8:4 Â 10 3 ; We ¼ 8:7 Â 10 À3 , and d ¼ 0:009.As shown in Fig. 3, the bubble expands relatively spherical except for the flattening of the surface nearest the closer wall.Examining frames four and five, as the bubble collapses with the near wall surfaces remaining in contact with the walls, the distal surface collapses inwards and a jet begins to form.Finally, at collapse, a wide jet forms and penetrates the bubble surface pointing toward the corner.The BIM computations are in very good agreement with the experimental images.

C. Convergence test
To analyze the convergence of the BIM model, convergence tests were conducted for a bubble near a corner of angles a ¼ p=2 and p=4.As shown in Fig. 4, the jet profiles for numbers of elements M ¼ 5412 and M ¼ 5724 are excellent in both cases.As such, all the remaining calculations in this study are performed using M ¼ 5412.

V. NUMERICAL RESULTS
First, we look at the evolution of a microbubble near corners of angle a ¼ p=2 and a ¼ p=4 in a Maxwell fluid with Deborah number De ¼ 2.0.The stand-off distances are given by c N ¼ 1:0 for the nearer horizontal wall below the bubble and c F ¼ 1:5 for the further vertical wall to the left of the bubble.As shown in Fig. 5, the bubble expands approximately spherically until it reaches its maximum volume at t ¼ 0.87.Due to the Bjerknes forces from the two rigid walls, the majority of the non-spherical collapse occurs on the side of the bubble furthest away from the corner.So, as the jet starts to form at t ¼ 1.61, the sides of the bubble nearest to the walls remain relatively spherical while the side furthest from the corner becomes flatter.Note that there is more flattening on the side of the bubble closest to the far wall as opposed to the side of the bubble closest to the near wall.The collapse continues in this fashion until the bubble becomes toroidal at t ¼ 1.76 with the jet rotated toward the far wall.As can be seen, the bubble has migrated toward the near wall with the bottom of the bubble remaining nearly spherical.
Figure 6 shows the evolution of a bubble in a Maxwell fluid near a corner of angle a ¼ p=4.The evolution is similar to the bubble near a corner of angle p=2 with the collapse occurring mostly on the bubble surface furthest away from the corner, and the jet rotated toward the far wall.Note the increase in the oscillation period is from 1.76 to 1.96.
We consider the shape of the bubble profile in a right-angled corner at collapse.The jet shapes are given in Fig. 7 for De ¼ 0.0, 0.2, 2.0, and 20.0, with the black lines below and to the left side of the bubbles representing the location of the rigid boundaries of the corner.The jet shape is broader, and the bubble is much larger as the jet penetrates through the surface as the Deborah number increases.This is a result of the higher maximum radius as De increases, as there is a greater Bjerknes attraction between the bubble and the nearest wall below the bubble.This is apparent from the fact that the bottom of the bubble is much closer to the near wall as the Deborah number increases.Furthermore, because of this Bjerknes forces of attraction, the jetting angle is more toward the far wall to the left as the Deborah increases, with a more asymmetric profile overall, as the jet becomes more angled.
The bubble profile at collapse for a corner of angle a ¼ p=4 is provided in Fig. 8.As shown, the collapse follows a similar trend for a ¼ p=2.As the Deborah number increases, the jet is broader and the bubble collapses closer to the nearer boundary.Once again, the jet is rotated upwards toward the far wall, an effect that is more apparent here given the restricted domain size.
The Kelvin impulse for a bubble is defined by The Kelvin impulse corresponds to apparent inertia of the bubble, and its direction indicates the direction of bubble migration and jetting. 67,68he Kelvin impulse toward the near and far boundaries over the first oscillation and subsequent bubble collapse is given, respectively, in Figs. 9 and 10 for De ¼ 0.0, 0.2, 2.0, and 20.0.As shown, an increase in the Deborah number results in a large increase in the absolute value of the Kelvin impulse and, therefore, corresponds to a higher apparent inertia toward both boundaries.We now consider bubble migration toward both boundaries in the corner.The bubble centroid migration toward the near boundary is shown in Fig. 11.As shown, an increase in the Deborah number results in a greater migration of the bubble toward the near boundary.This point is apparent from Fig. 7 where it is clear that the bubble is closer to the near boundary at collapse as De increases.There is also the bubble migration toward the far wall.This is shown in Fig. 12.As can be seen, there appears to be not much of a difference between the bubble migration in this direction and what was seen in the direction toward the near wall.This is because the bubble centroid motion toward the far wall is comparable in all cases due to the bubble being larger as the Deborah number increases.Note from Fig. 7 that indeed the left hand side of the bubble is closer to the vertical wall, but so also the right hand side of the bubble is further away as the Deborah number increases.So the bubble centroid in this direction is roughly similar for all Deborah numbers.
Figure 13 compares the equivalent bubble radius histories for a ¼ p=2 and a ¼ p=4 for De ¼ 0.0 and 20.0.The maximum bubble radius and oscillation period increase with De.As a increases, the amplitude does not change significantly but the period decreases.
Figure 14 shows the jet shape for bubbles at collapse in viscoelastic fluids of Deborah number De ¼ 20.0 and 200.0.As shown, an increase from De ¼ 2.0 to 20.0 results in a slightly larger bubble and a jet angle pointed more toward the far wall.However, increasing the Deborah number beyond De ¼ 20.0 results in minimal change to the bubble profile at collapse.This can be seen by observing that the bubble profile for De ¼ 200.0 has no significant difference from the bubble profile for De ¼ 20.0.From Maxwell's constitutive equation (12), we can see that DT nn =Dt is inversely proportional to De.As De is large, DT nn =Dt and T nn are small and negligible.To see this, let us observe Eq. ( 12).As De increases, the Maxwell model converges to the non-viscous model where T nn % 0 for all time.
Furthermore, Fig. 15 shows the corresponding radius history for the cases considered in Fig. 14.As shown, there is only a slight increase in the radius and oscillation period between De ¼ 20.0 and De ¼ 200.0.

VI. CONCLUSION
Bubble dynamics in non-Newtonian liquids has implications for a wide range of applications in engineering and biomedical ultrasonics.A numerical model is described based on the boundary integral method coupled with the Maxwell constitutive equation.Here, the flow is assumed to have high Reynolds number and thus is irrotational in the bulk domain.The viscoelastic effects are not negligible in a thin viscous layer at the bubble surface and are approximated only through the normal stress balance at the bubble surface.An Euler method is used to iteratively solve the Maxwell constitutive equation for the normal component of stress.Computations are carried out for a bubble near a corner for various corner angles using Green's function for the Laplace equation in the corner domain.
The numerical model has been validated and excellent agreement with predictions obtained using the corresponding Rayleigh-Plesset equation for spherical bubbles in a non-Newtonian liquid is demonstrated.The numerical results also agree very well with experiments for bubbles in a Newtonian fluid near a corner.Viscoelastic effects are studied in terms of Deborah number De, which is as the ratio of the time it takes for a material to adjust to applied stresses to the characteristic bubble oscillation time.Simulations were performed for De ¼ 0.0 (Newtonian fluid), De ¼ 0.02, 2.0, and 20.0 (viscoelastic fluid).The following new phenomena/features were observed: 1.The maximum radius and oscillation period increase with the Deborah number.This is due to more elastic energy being stored and released in the viscoelastic fluid for a larger Deborah number and less energy being lost due to viscosity.2. The wall effects due to the corner are enhanced with a larger De as a result of an increase in the maximum radius.Bubble migration and oscillation amplitude increase with the Deborah number.3. The bubble jet shape is broader and flatter for a bubble near a rigid boundary or in a corner as De increases.The bubble jet is inclined to the further boundary when it is in a corner as De increases.4. Results converge as the Deborah number increases as the elastic effects become dominant compared with the viscous damping effects of the fluid.
In conclusion, an increased Deborah number greatly enhances the effects of the surrounding boundaries on bubble collapse and provides a broader, flatter jet.This implies that a higher Deborah number in the fluid will enhance the interactions between a bubble and boundaries or other bubbles.This also suggest that, with certain parameter ranges, elasticity and non-Newtonian fluids may benefit ultrasonic cleaning, therapy, and drug delivery.

FIG. 5 .
FIG. 5.The evolution of a three-dimensional bubble near a corner of angle a ¼ p=2 in a Maxwell fluid for De ¼ 2.0 with the remaining parameters being the same as in Fig. 1.

FIG. 7 .
FIG. 7. The influence of De on the jet profiles and location of a bubble initiated in a corner of angle a ¼ p=2, with the remaining parameters being the same as in Fig. 1.

FIG. 4 .
FIG.4.The jet profiles and location of a bubble initiated with 5412 nodes and 5724 nodes for angles p=2 and p=4 with the remaining parameters being the same as in Fig.1.

FIG. 9 .FIG. 10 .FIG. 12 .
FIG. 9.The influence of De on the Kelvin impulse toward the near wall, for a bubble near a corner of angle a ¼ p=2 with c N ¼ 1:0 and c F ¼ 1:5.The remaining parameters are the same as in Fig. 1.

FIG. 13 .
FIG. 13.The evolution of the bubble radius in an corner of angle a ¼ p=2 (black) and a ¼ p=4 (red) with De ¼ 0.0 and 20.0 with the remaining parameters being the same as in Fig. 1.

FIG. 15 .
FIG. 15.The influence of large De on the radius history, for a bubble near a corner of angle a ¼ p=2 with c N ¼ 1:0 and c F ¼ 1:5.The remaining parameters are the same as in Fig. 1.