Liquid bridge simulations with OpenFOAM

Liquid bridge simulations with OpenFOAM Guillermo Guirao Aguilar SUPERVISED BY Hendrik Kuhlmann Universitat Polit` ecnica de Catalunya Master in Aerospace Science & Technology May 2011 Liquid bridge simulations with OpenFOAM BY Guillermo Guirao Aguilar DIPLOMA THESIS FOR DEGREE Master in Aerospace Science and Technology AT Universitat Polit` ecnica de Catalunya SUPERVISED BY: Hendrik Kuhlmann Technische Universit¨ at Wien I would like to thank Hendrik Kuhlmann and Laureano Ram´ ırez de la Piscina for their kindness and support, Frank Muldoon for sharing his knowledge on the topic with me, and specially Ernst Hofmann without whom this work would not have been possible. ABSTRACT Any liquid-gas interface being under a temperature gradient is subject to a thermocapillary flow generated by the differences in the surface tension. In the particular case of liquid bridges, this flow can evolve into a oscillatory or travelling flow known as hydrothermal wave. This works aims to achieve a better understanding of the formation of this kind of flow by means of computational fluid dynamics techniques using a specific software package. For that, a simplified model of the liquid bridge is proposed, reducing the problem its most elementary parameters. The numerics used to solve the system make use of the NavierStokes equations for fluid mechanics. In order to validate the simulations developed, a series of analysis will be performed to test the accuracy and reliability of the results obtained, as well as to test the capabilities of OpenFOAM running this kind of problems. The software employed to run the simulations is called OpenFOAM (Open Field Operations And Modifications), which is a free and open-source CFD package with a wide range of usability in many engineering and science fields. Keywords: Computational fluid dynamics, liquid bridge, OpenFOAM, thermocapillary flow, hydrothermal wave, Marangoni effect CONTENTS ABSTRACT ...................................... 7 Preamble ...................................... 1 1. Introduction .................................. 3 1.1. Liquid bridges .................................. 3 1.2. Thermocapillary flow .............................. 4 1.3. The hydrothermal wave ............................ 5 2. Numerics .................................... 7 2.1. The Navier-Stokes equation .......................... 7 2.2. PISO algorithm ................................. 8 2.3. Rhie-Chow interpolation ............................ 10 2.4. Solver ...................................... 11 3. 2D case ..................................... 13 3.1. Set up ...................................... 13 3.2. Analysis ..................................... 15 4. 3D case ..................................... 19 4.1. Initial conditions ................................ 19 4.2. Grid ....................................... 19 4.3. Simulations ................................... 21 5. Analysis ..................................... 25 5.1. Grid convergence ................................ 25 5.2. Analysis for different Re ............................ 28 5.3. Parallel processing ............................... 30 4 Liquid bridge simulations with OpenFOAM Figure 1.1: Half-zone (a) and full-zone (b) models. Taken from Leypold [5] Difficulties on developing tall liquid bridges (large aspect ratio) and large liquid bridges (large radius) in terrestrial experiments avoid further research on such configurations. However, low aspect ratio liquid bridges with Γ1are easier to produce, and a lot of research has been done from which we can obtain the data to compare with the numerical results. Therefore, an aspect ratio of Γ=0.66 will be used for the present work. 1.2. Thermocapillary flow Any spatial gradient of the surface tension due to nonuniform temperature distribution along an interface inevitably generates a fluid motion known as thermocapillary flow. Under these conditions, the fluid is forced to move from hotter toward the colder end wall along the free surface, resulting in a return flow through the inner part of the liquid bridge. This is known as the Marangoni effect. Thermocapillary flows can be characterized by the Marangoni number, which at the same time can be defined as the product of the thermocapillary Reynolds number and the Prandtl number, such that Ma =RePr =γd ρ0ν2 ν κ∆T=−∂σ ∂T d µκ∆T(1.2) The Reynolds number establishes the ratio between inertial forces and viscous forces and consequently quantifies the relative importance of these two types of forces for given flow conditions Re =γd ρ0ν2∆T(1.3) Introduction 5 while the Prandtl number states the ratio of kinematic viscosity and the thermal diffusivity Pr =ν κ(1.4) In the equations above we can find the surface tension coefficient (γ) factor and the kinematic viscosity (ν), respectively defined as γ=−∂σ ∂T,ν=µ ρ0 (1.5) where σis the surface tension, Tis the absolute temperature, µis the dynamic viscosity and ρ0is the fluid density. With decreasing volume, the surface to volume ratio increases. Hence, the capillarity and the thermoand soluto-capillary convections play more important roles in the micro-fluid hydrodynamics. For the liquid free surface, we assume Bo 1and Ca 1, being the former the Bond number, also known as the Etv¨ os number, which relates the importance of the buoyancy forces versus the surface tension, while the Capillary number relates the relative effect of viscous forces versus surface tension forces. These number become vanishingly small as the gravitational force, which drives the buoyancy flows, can be neglected in small scale geometries common in liquid bridge, or as in this case, in a zero-gravity environment. Therefore a non-deformed surface is a good assumption that greatly simplifies the problem while keeping the simulation close to the real model. These number are defined as Bo =ρ0gd2 σ0 ,Ca =γ∆T σ0 (1.6) 1.3. The hydrothermal wave Existing researches have revealed that a steady two-dimensional flow can undergo a transition to a three-dimensional (azimuthally non-uniform) either stationary or oscillatory when the temperature difference ∆Tachieves a critical value ∆Tcr (Schwabe et al. [1]). A necessary condition for that to occur is to have a relatively small applied temperature difference compared to the absolute temperature of the system. By changing the Reynolds number, the ∆Tcis affected proportionally. Likewise, the critical Marangoni number Maccan be related to these two as Mac=RecPr ∝ δTc(1.7) The hydrothermal instability is oscillatory and starts as a result of a supercritical Hopf bifurcation as either travelling or standing wave, often simply referred to as hydrothermal wave. Its angular velocity approaches asymptotically a constant value. 6 Liquid bridge simulations with OpenFOAM Figure 1.2: Different azimuthal wave number flows in terms of the aspect ratio. Taken from Kawamura [2] The three-dimensional flow shows a modal structure described by an azimuthal wave number m. The latter has been found to depend both on the aspect ratio Γof the liquid bridge and the temperature difference (Schwabe et al. [3]). Several experiments has been made in order to find the conditions under which the liquid bridge evolves into one of these structured flows. Low Prandtl numbers are preferred for numerical simulations, but these are mostly characteristic of metals and semiconductors, and not commonly employed in experiments as the flow cannot be seen through them, so less data is available. However, numerical modeling for Pr <7predicted m≈2.0/Γ(see Fig. 1.2), as explained in Melnikov et al. [4]. Numerics 7 CHAPTER 2. NUMERICS 2.1. The Navier-Stokes equation Since we are dealing with a liquid fluid, an incompressibility assumption is feasible. Therefore the continuity and momentum equations are given by ∇·u=0(2.1) ∂u ∂t+u·∇u=−∇p+ν∇2u(2.2) The left hand side of the equation 2.2 represents the inertia per volume of the fluid, where u·∇uis the convective acceleration, while the right hand side is the divergence of stress divided into the pressure gradient ∇pand the viscosity ν∇2u This is written in OpenFOAM as: fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) -fvm::laplacian(nu, U) ); On it, the three main operators can be easily identified, being these the derivative with respect to time, the divergence (ddt, and div respectively) and the laplacian. However, there is no right hand side, and there is a field named phi (φ). This term is (for incompressible flows) the volume velocity flux defined on the faces of each cell, and it is used because OpenFOAM can utilize the Gauss theorem, which is frequently used in applied mathematics, and defines the transform of a volume integral into a surface integral. ZV ∇·(uϒ)dV =ZS (uϒ)f·ˆ ndS (2.3) =∑ i uf,iϒf,i·S fi=∑ i uf,iφi(2.4) where φ=ϒf·S f (2.5) ϒis the velocity that will be held constant when the equation for pressure is solved, while uis the vector velocity that will be solved for. It is important to note the difference in the 8 Liquid bridge simulations with OpenFOAM subscript when the surface integral is introduced; subscript findicates that the term should be evaluated on the face. φis defined as the scalar product of the cell face velocity and the cell face normal (Eq. 2.5). The magnitude of the cell face normal is the cell face area. We can now recognize OpenFOAM’s divergence in the left hand side of the equation 2.4. 2.2. PISO algorithm The solution of these equations is complicated by the lack of an independent equation of the pressure, whose gradient contributes to each of the three momentum equations. Furthermore, the continuity equation does not have a dominant variable in incompressible flows. Mass conservation is a kinematic constraint on the velocity field rather than a dynamic equation. One way out of this difficulty is to construct the pressure field so as to guarantee satisfaction of the continuity equation. It must be noted that the absolute pressure is of no significance in an incompressible flow; only the gradient of the pressure (pressure difference) affects the flow. To do that, several implicit iterative methods can be used. Many of these methods for steady problems can be regarded as solving an unsteady problem until a steady state is reached. The principal difference is that, when solving an unsteady problem, the time step is chosen so that an accurate history is obtained while, when a steady solution is sought, large time steps are used to try to reach the steady state quickly. Implicit methods are preferred for steady and slow transient flows, because they have less stringent time step restrictions than explicit schemes. The method employed for this solver is a derivative of the SIMPLE algorithm called PISO (Pressure-Implicit with Splitting of Operators). It uses a pressure (or pressure-correction) equation to enforce mass conservation at each time step. Rather than solve all of the coupled equations in a coupled or iterative sequential fashion, the PISO algorithm splits the operators into an implicit predictor and multiple explicit corrector steps. Very few corrector steps are necessary to obtain desired accuracy. So firstly, to obtain the pressure equation for an incompressible flow, it has to be derived from the momentum and continuity equations. To start, the momentum equation is discretized: Aui Pun+1 i,P+∑ l Aui lun+1 i,l=Qn+1 ui−∂pn+1 ∂xiP (2.6) where P is the index of an arbitrary velocity node, Ais a coefficient, mis the number of the current iteration, and the index ldenotes the neighbor points that appear in the discretized momentum equation. The source term Qcontains all of the terms that may be explicitly computed in terms of un ias well as un+1 i. Due to the non-linearity and coupling of the underlying differential equations, the equation 2.6 cannot be solved directly as the Numerics 9 coefficients Aand, possibly the source term, depend on the unknown solution un i+1. An iterative approach is the only choice. The iterations executed within one time step, in which the coefficient and source matrices are updated, are called outer iterations to distinguish them from the inner iterations performed on linear system with fixed coefficients. On each outer iteration, the equations solved are: Aui Pum∗ i,P+∑ l Aui lum∗ i,l=Qm−1 ui−∂pm−1 ∂xiP (2.7) As the pressure used in these iterations was obtained from the previous outer iteration or time step, the velocities computed from equation 2.7 do not normally satisfy the discretized continuity equation. To enforce the continuity condition, the velocities need to be corrected; this requires modification of the pressure field; the manner of doing this is described next. First off, the velocity at node P, obtained by solving the linearized momentum equations 2.7, can be formally expressed as: um∗ i,P=Qm−1 ui−∑lAui lum∗ i−l Aui P − 1 Aui P∂pm−1 ∂xiP (2.8) For convenience let’s assume: ˜um∗ i,P=Qm−1 ui−∑lAui lum∗ i,l Aui P (2.9) As previously stated, these velocities do not satisfy the continuity equation and must be corrected. A pressure-correction can be used instead of the actual pressure. The velocities computed from the linearized momentum equations and the pressure pm−1are taken as provisional values to which a small correction must be added: um i=um∗ i+u0(2.10) and pm=pm−1+p0(2.11) If these are substituted into the momentum equation, we obtain the relation between the velocity and pressure corrections by means of the SIMPLE method: u0 i,P=˜u0 i,P− 1 Aui P∂p0 ∂xiP (2.12) 10 Liquid bridge simulations with OpenFOAM In the PISO algorithm, a small time-step is assumed, hence the pressure-velocity coupling is much stronger than the non-linear coupling, and therefore is possible to repeat a number of pressure correctors without updating the discretization of the momentum equation. In such a setup, the first pressure corrector will create a conservative velocity field, while the second and following will establish the pressure distribution. Continuity is enforced by inserting this expression for um iinto the continuity equation to yield a discrete Poisson equation for the pressure: ∂ ∂xiρ Aui P∂p0 ∂xiP =∂(ρ˜u0 i) ∂xiP (2.13) After solving this equation for the pressure, the final velocity field at the new iteration um i is calculated from equation 2.12. At this point, we have a velocity field which satisfies the continuity condition, but the velocity and pressure fields do not satisfy the momentum equations. A number of iterations given as an input parameter in OpenFOAM is performed until a velocity field which satisfies both the momentum and continuity equations is obtained. Since multiple pressure correctors are used with a single momentum equation, it is not necessary to under-relax neither the pressure nor the velocity. On the negative side, the derivation of PISO is based in the assumption that the momentum discretization may be safely frozen through a series of pressure correctors, which is true only at small time-steps. Experience also shows that the PISO algorithm is more sensitive to mesh quality than the SIMPLE algorithm. 2.3. Rhie-Chow interpolation The Rhie-Chow interpolation is a necessary step when using a colocated finite volume method formulation, since it removes oscillations in the solutions. These oscillations occur if the pressure gradient does not depend on the pressure in adjacent cells, and thus allowing a jigsaw pattern. It is defined as a correction proportional to the difference between the pressure gradient at the face and the interpolated pressure gradient at the face. This would be as follows for a velocity correction: uj=uj−∆ 1 Auj p! ∂p ∂xj −∂p ∂xj!(2.14) where the overbar indicates interpolation, and ∆is related to the mesh size. But such term is not explicitly found in OpenFOAM. A work around is done instead by replacing the corrected velocity from equation 2.12 by the velocity flux φ, since the face velocities will be used to evaluate the term. The resulting equation is then written as follows in OpenFOAM: Numerics 11 volScalarField rUA = 1.0/UEqn.A(); U = rUA*UEqn.H(); phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); fvScalarMatrix pEqn ( fvm::laplacian(rUA, p) == fvc::div(phi) ); As before, OpenFOAM makes use of the Gauss theorem, and therefore it is not necessary to calculate a second derivative of p, but only a first derivative. 2.4. Solver For the sake of clarity, the scheme of the calculations made by the solver for each time step is presented here: The conservative fluxes derived from the previous time step are used to discretize the momentum equation. The momentum equation is solved using the pressure from the previous time step, obtaining the momentum predictor step. The PISO loop starts: •The velocity field is computed without the pressure gradient and the interpolated face fluxes are calculated from the approximate velocity field (corrected to be globally conservative so that there is a solution to the pressure equation). •An inner loop for the non-orthogonal corrector is ran for a fixed number of times. In the final iteration, φis finally corrected for the next pressure-corrector step. •The continuity error is calculated. •The approximate velocity field is corrected using the new pressure gradient. •The calculation of the pressure-corrector is repeated until the continuity equation is satisfied. Finally the equation for the temperature is solved. At this point, all the parameters for the current time step are known and a new iteration can be started. 12 Liquid bridge simulations with OpenFOAM 2D case 13 CHAPTER 3. 2D CASE 3.1. Set up To test the behavior of the solver with the numeric scheme explained in the previous chapter, a simple bidimensional case was ran. Considering that OpenFOAM can’t handle pure bidimensional cases, the grid had to be made with one cell in depth, being the third dimension perpendicular to the plane of the two dimensions of the problem. The boundary conditions will be set afterwards in such a way that the data is written only at the cell centers and cell faces at the top, bottom and free surface boundaries, as in a real 2D case. To approximate the cylindrical geometry of the original case as close as possible, one of the side faces of the grid were collapsed into one edge, resulting in a wedge shaped grid. This edge corresponds to the axis of the cylinder and the outer face is the free surface. This grid is made of one layer of structured hexahedra along the plane, with a cell compression factor in the radial and axial directions. The compression factor is made to refine the mesh gradually in a given direction, and it is defined as the ratio between the start δs and end δecells lenghts in that direction: Compression ratio =δs δe Compression ratios of 5 were applied in the radial direction toward the free surface, and 2 from the middle cross section toward the top and bottom walls. A refinement of these corners of the cylinder is necessary to obtain a better resolution of those areas where the gradients and velocities are much higher. The inner area around the axis of the cylinder has little movement, hence the use of bigger cells helps reducing computational effort without affecting the behaviour of the flow. An example of the geometry of the grid is shown in figure 3.1 for a relatively low number of cells. To the end of making easier comparing results among simulations, the main parameters have been adimensionalized. This way, the range of values will remain the same for all the cases, whatever are the real values corresponding to them. For that, scaling constants are needed to convert from the dimensional values to the dimensionless ones. In the case of temperature and pressure, these are P0=ρ0U2 0=ρ0ν2 d2,∆T=Ttop −Tbot (3.1) Therefore, their respectives non-dimensional values are p=p∗ P0 ,T=T∗−T0 ∆T(3.2) 20 Liquid bridge simulations with OpenFOAM Figure 4.1: Top view of the perturbed temperature field at z=0 in such cases. To ensure the grid independence of these results, the axisymmetric grid arrangement shown in Fig. 4.3 was selected as a first choice. Nevertheless, a problem arises for this kind of grid due to the narrowing of the cells near the center and outer surface of the cylinder. As the faces of the cells must match one to one with its neighboring cells, a reduction of the number of cells around the axis (or otherwise increase of cells near the free surface) can’t be performed with such distribution. The size and shape of the cells is conditioned by the Courant-Friedrichs-Lewy condition. This condition states that the time a travelling particle takes to cross a cell should be less than the time step used in the calculation in order to assure that the scheme can access the information required to form the solution. Otherwise the simulation could produce wildly incorrect results, and eventually diverging in the iterative process. Therefore, the size of the cells and velocity of the flow can be related with the time step as U·∆t ∆x=ν≤C(4.2) being νcalled the Courant number, which as a rule of thumb should remain below 1. According to this expression, if a cell is stretched in the direction perpendicular to the flow, the time step needed to perform the calculation becomes smaller while the resolution of the resulting field of parameters is as bad as the distance of the long side of the cell. For this reason the slender cells should be avoided in non-laminar flows. 3D case 21 Figure 4.2: Side view of the perturbed temperature field Grid B (Fig. 4.4) minimizes this problem by changing the distribution of the inner area into better shaped cells. But non-orthogonal grids demand another loop within the PISO algorithm as explained before, meaning some extra computational time per time step. Nonetheless, the time step needed decreases by two orders of magnitude with respect to the one used for grid A. So the overall savings in time for the whole computation fully justifies this change. After several runs only a few cases of those using grid B arose an even number of lobes, differing from the results obtained with grid one, concluding that the computational effort saved by using grid B was more significant than the proportion of affected results. Both grids were built with the same cell compression factor employed in the bidimensional grid. 4.3. Simulations Using the same parameters as for the bidimensional case, the critical Reynolds number is known to be Rec=1080 (Hofmann [6]). The grid used has a total of 737280 cells, and two different Reynolds numbers were used, one for a subcritical case (Re =1000) and one for a supercritical (Re =1800). In both cases, the time needed to reach a stable 2D flow or travelling wave respectively was of two dimensionless units, i.e. half of the momentum diffusion time. The temperature at a probe cell has been plotted versus time. From it, the angular velocity of the hydrothermal wave can be calculated from the periodic changes of temperature once the flow is fully developed. For the present case, it was found to be Ω=10.235 in a counter-clockwise direction. 22 Liquid bridge simulations with OpenFOAM Figure 4.3: Top view of grid A Figure 4.4: Top view of grid B 3D case 23 Figure 4.5: Temperature field at the middle section of a liquid bridge showing m=4for Γ=0.66 0 0.5 1 1.5 2 2.5 3 3.5 4 −0.2 −0.1 0 0.1 0.2 Dimensionless time Temperature Figure 4.6: Evolution of the temperature at r=1,φ=0and z=0 24 Liquid bridge simulations with OpenFOAM Analysis 25 CHAPTER 5. ANALYSIS 5.1. Grid convergence In order to obtain reliable results, one must ensure that the cells size is small enough to avoid the accumulation of errors. If the distance between cells is too large, the gradients lying in between will not be present in the calculation, leading to some error in the parameters around it. To avoid that, the cell must be refined until the gradient is small enough to be neglected and to not produce any appreciable difference in the results of the neighboring cells. To determine the size of the cells from which the results become independent to further refinements of the grid, an analysis of grid convergence has been made. In this analysis, the same case from previous chapter was ran, increasing gradually the refinement of the grid, keeping constant the ratio between number of cells along the three axis, such that N(r,φ,z)0 N(r,φ,z)1 =constant This ratio was doubled each time, increasing eightfold the number of total cells in the grid. A total of five different grids were made ranging from 1,440 to 5,898,240 cells, being the former too small to even develop the hydrothermal wave, and the latter big enough to take a considerable amount of time to develop a short period of the simulation. Merely the Reynolds number was changed, using three different values to compare the grid convergence for subcritical and supercritical values. In the first place, the subcritical case was evaluated using a Re =1000. Knowing that the 2D flow is the only possible solution, we can plot the field parameters over a line, as it was done for the 2D case, and compare the differences between grids giving a basic idea of how rough are the results for a coarse number of cells. The samples taken lie again in R=rfor z=0and h=zfor R=1.515. Results are shown in Figs. 5.1-5.4. The line corresponding to the first grid was omitted for being too inaccurate, and useless for the purpose of this test. It can be clearly seen how the difference between lines decreases as the number of cells becomes larger, approaching to the real values. In the case of the axial velocity along the free surface, an upward peak appears near the bottom wall for the coarse grids, differing from the line shown for fine grids. This is a deviation of the real values caused by the relative coarseness of the cells in an area where the local velocities are very high. This high velocities form a jet stream originated by the sudden change in direction of the flow going downwards along the free surface and hitting the bottom wall of the liquid bridge. 26 Liquid bridge simulations with OpenFOAM 0 0.5 1 1.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 r T 11520 cells 92160 cells 737280 cells 5898240 cells Figure 5.1: Temperature distribution along z=0 0 0.5 1 1.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 r T 11520 cells 92160 cells 737280 cells 5898240 cells Figure 5.2: Urdistribution along z=0 Analysis 27 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 z T 11520 cells 92160 cells 737280 cells 5898240 cells Figure 5.3: Temperature distribution along r=1.515 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 −60 −50 −40 −30 −20 −10 0 z Uz 11520 cells 92160 cells 737280 cells 5898240 cells Figure 5.4: Uzdistribution along r=1.515 28 Liquid bridge simulations with OpenFOAM 5.2. Analysis for different Re An analysis has been made to determine how an increasing Reynolds number affects the formation of the hydrothermal wave. Using the same parameters as before, four different Reynolds numbers were employed, all of them supercritical: 1800, 2400, 4000, 8000. From the pictures shown in Figs. 5.5–5.7, an increase in the temperature gradients can be appreciated. The angular velocity has been observed to change as well, obtaining Ω1800 =10.235 Ω2400 =10.277 and Ω4000 =10.310. Beyond a given Reynolds number, the flow exhibits a chaotic behavior. The surface temperature variation mostly loses its periodic nature, and its power spectrum broadens. This was the case for Re =8000. Figure 5.5: Middle cross section for Re =1800 Analysis 29 Figure 5.6: Middle cross section for Re =2400 Figure 5.7: Middle cross section for Re =4000