CC BY
11
DOI: 10.17516/1997-1397-2020-13-6-670-677 УДК 532.542
Three-Dimensional Simulation of a Tank Filling With a Viscous Fluid Using the VOF Method
Evgeny I. Borzenko* Efim I. Hegaj1"
Tomsk State University Tomsk, Russian Federation
Received 13.06.2020, received in revised form 04.09.2020, accepted 04.10.2020 Abstract. This paper presents the results of 3D modeling of a Newtonian fluid flow with a free surface. The PLIC-VOF algorithm, which is developed to solve the problems of two-dimensional fluid flows with a free surface, is generalized to the case of three-dimensional flows. Efficiency of the developed algorithm and reliability of the obtained results are justified by comparing with available data in literature and by testing approximation convergence.
Parametric calculations of a rectangular channel filling show that the free surface assumes a steady convex shape over time and then moves along the channel at a constant velocity. As a result of parametric studies, the dependences of geometric characteristics of the free surface shape on problem parameters have been plotted.
Keywords: Newtonian fluid flow, filling of a rectangular channel, free surface, 3D modeling, numerical simulation, VOF method, flow structure.
Citation: E.I. Borzenko, E.I.Hegaj, Three-Dimensional Simulation of a Tank Filling With a Viscous Fluid Using the VOF Method, J. Sib. Fed. Univ. Math. Phys., 2020, 13(6), 670-677. DOI: 10.17516/1997-1397-2020-13-6-670-677.
Introduction
Technological processes associated with casting of items, filling of tanks or draining of polymer compounds are characterized by the presence of a free surface in a fluid flow. Adequate engineering for such processes requires a detailed study of the flow features and free surface behavior [1].
Tracking of free surface evolution in a fluid flow is known as a complex problem in hydromechanics [2]. One of the first successful attempts to determine free surface dynamics in the two-dimensional approximation has been made in [3]. In this paper, a numerical method based on the VOF (Volume of Fluid) approach is proposed, which allows one to determine a free surface position at any time instant using a scalar function defined in the cells of a regular grid. Moreover, at a discrete level, the free boundary in a control volume is represented as a segment, which is parallel to one of the volume's faces. The method was further developed in [4]. A modification of the method, namely PLIC-VOF (Piecewise-Linear Interface Calculation), is proposed, according to which a free surface in a control volume is represented as a set of arbitrary oriented segments. The PLIC-VOF method has been successfully applied for two-dimensional flows in [5-10].
Many fluid flow features cannot be taken into account nor adequately assessed when modeling the process in the plane or axisymmetric approximation. Two-dimensional problem formulations
* [email protected] https://orcid.org/0000-0001-6264-1776 [email protected] https://orcid.org/0000-0002-7973-9435 © Siberian Federal University. All rights reserved
allow one to study flow mechanics only for some channel designs. However, most of the real pipelines have complex three-dimensional geometry. Therefore, there is a need to develop algorithms for calculating and modeling fluid flows in a full three-dimensional formulation.
In this work, testing of a modified VOF method on the three-dimensional fluid flow modeling is implemented.
1. Formulation of the problem
A three-dimensional flow, which occurs when a vertical rectangular channel is being filled in a gravity field, is considered. A solution domain is schematically shown in Fig. 1 a. In this case, the fluid is supplied from the bottom through the inlet section at given constant flow rate.
At the initial time instant, the channel is partially filled with a fluid, and the free surface is represented as a plane 2 = const confined by the walls.
When filling the channel, the free surface becomes curved and assumes a convex shape. The maximum height of the free surface, x (Fig. 1 b), is taken as convexity characteristics, which is determined by the values of the parameters Re and W.
Fig. 1. Solution domain (a) at the initial time instant and (b) during the filling process
Mathematical formulation of the problem includes the Navier-Stokes and continuity equations written in a dimensionless form as
Re Re Re
dUx ~dt~ + Ux dUx dx + Uy dUx dy + Uz dUx dz
dUy ~dt~ + Ux dUy dx + Uy dUy dy + Uz dUy dz
U dt + Ux U dx + Uy U dy + Uz dUzj dz
= dP
dx
= dP dy
+
+
dz2
id2ux , d2ux , d2a
dx2
( d2 uy
dx2
+
dy2
d2Uy d2Uy
dy2
= dP (d2Uz , d2uz dz
+
dUx dUy
duz
H--- +
dx dy dz
dx2
0
+
dy2
+
dz2
d2Uz dz2
(1)
- W
(2)
No-slip conditions are assigned on the solid walls. On the free surface, the continuity conditions for normal and shear stresses are used. The fluid is supplied through the inlet section at a velocity equal to unity.
The following quantities are used as length, velocity, time, and pressure scales: L (a characteristic size of the channel), U0 (an average velocity at the inlet section), and the complexes of L/U0 and ¡iU0/L, respectively. The problem formulation includes dimensionless criteria: the
Reynolds number Re = pU0L/^ and the parameter W = pL2g/^U0 = Re/Fr, which is equal to the ratio of the Reynolds and Froude numbers.
2. Method of solving
The formulated problem is solved numerically using the finite volume method implemented on a staggered grid. Kinematic and dynamic characteristics of the flow are determined using the SIMPLE algorithm [11]. In this case, to approximate of convective and non-stationary term, an exponential scheme was used. Tracking of the free surface evolution is carried out by a modified VOF method, which is generalized to the three-dimensional case with account for an arbitrary inclination of the free surface in a control volume.
The VOF method implies introducing of a scalar function F, whose value is equal to unity at all points occupied by the fluid and equal to zero at the rest of the points. At a discrete level, when averaging over the control volume, the value of F is equal to a volume fraction of the control volume occupied by the fluid. In particular, when the control volume is entirely filled with a fluid, F = 1, and when the control volume does not contain any fluid, F = 0. If there is a free surface in the control volume, 0 < F < 1. The values of this function can be determined from the equality to zero of the total derivative of F with respect to time, which reflects the law of conservation of mass
dF dF U dF U dF
dt x dx y dy z dz
When integrating this equation over the control volume, it is necessary to determine the fluxes through the faces, which are calculated using the obtained velocity values on these faces and the orientation of the free surface in the control volume.
It is assumed that at a discrete level, the free surface in the control volume represents a cutting plane, whose position is determined by its normal and by the fraction of the volume filled with the fluid. In the three-dimensional case, there are eight options for the free surface to cross the control volume (Fig. 2), which are characterized by positive components of the normal vector. The normal to the free surface is assumed to be directed along the gradient of the function F. When the value of F and the direction of the normal to the free surface are known, the plane approximating the free surface can be drawn in the boundary control volume. Therefore, in addition to the boundary control volume tracking, the function F is used to detect the fluid location inside the volume.
— + U^— + Uy— + Uz— = 0 (3)
Fig. 2. Possible locations for the free surface inside the control volume
Other possible free surface locations are reduced to those shown in Fig. 2 by rotating the control volume and / or by reflecting in coordinate planes.
Variation of the function F with time is determined by equation (3), which can be solved numerically after calculating its fluxes through the control volume faces. An illustration of the method in use is shown in Fig. 3. Two adjacent control volumes are considered, where the fluid
flows through a common side. The velocity U on the adjacent face determines which one is a donor and which is an acceptor. Afterwards, a plane, which is parallel to the common side, is plotted at a distance of UAt. The fluid fraction in the donor cell AF, which is enclosed between the plotted plane and the adjacent face, is transferred to the acceptor. The value of this fluid fraction (AF) is calculated by analytical formulas for polyhedra volume.
Donor Donor Acceptor
U At
Fig. 3. Illustration of the VOF method operation The fluxes through the other five faces of the control volume are determined similarly.
3. Verification of numerical method
To verify operational capability of the developed algorithm and reliability of the obtained results, the approximation convergence is tested on a sequence of grids.
The calculations show that at a time instant of t = 2 in a cross section of z = 2, the absolute value of the maximum transverse velocity is of the order of 10~6. It is supposed that the inlet boundary and the free surface do not affect the flow in this section, where a steady-state flow is observed. Thus, to verify the obtained velocity distributions, a well-known solution is used [12]:
Uz = 3.665 Y^
(-1)n
(2n +1)3
ch
l
(2n +1) n (2y - 1) 2
ch
(2n + 1) n 2
cos
(2n + 1) n (2x - 1)
(4)
The velocity error can be calculated as
AUz = max
Uz- Uz
UUz
100%.
(5)
Since the fluid velocity is equal to unity in the inlet section, at a time instant of t = 2, the volume of the inflowed fluid should equal V = 2. The error in the calculation of the volume is determined by formula
A^
2- Vt
t=2
2
100%.
(6)
To show the approximation convergence, a velocity profile along a straight line of y = 0.5 at z = 2 (Fig. 4 a) and a free surface shape in a cross section of y = 0.5 (Fig. 4 b) are calculated at different grid steps.
According to Fig. 4 b, a maximal difference in the free surface shapes is observed on the solid walls. Therefore, to control the approximation convergence, a parameter H (the free surface height at the point of x = 0, y = 0.5) is introduced.
Thus, the errors in the calculations of the fluid velocity, fluid volume, and free surface height on the solid wall H are selected as controlled characteristics in the computational method verification.
2
max
Fig. 4. (a) Velocity distribution along a straight line of y = 0.5 at 2 = 2 and (b) the free surface shape in a cross section of y = 0.5 at a time instant of t = 2 for Re = 0.1 and W = 32: h = 0.1 (the dashed line), h = 0.05 (the dotted and dashed line), and h = 0.025 (the solid line)
The obtained results presented in Tab. 1 demonstrate the approximation convergence for the selected characteristics.
Table 1. Dependence of the values of the controlled characteristics on the grid step at t = 2, Re = 0.1, W = 32
h AUZ, % A^, % H
0.1 2.53 0.42 3.944
0.05 2.15 0.41 3.915
0.025 2.01 0.40 3.900
4. Results of calculations
The initially flat free surface becomes curved over time and assumes a convex shape, which then moves upward along the channel and remains unchanged. Fig. 5 shows the evolution of the free surface shape.
Fig. 5. Free surface shapes for Re = 0.1, W = 32 at various time instants: (a) t = 0, (b) t = 0.2, (c) t = 0.5, (d) t = 1, (e) t = 2, and (f) t = 3
Fig. 6 demonstrates distributions of the velocities and pressure in a longitudinal section of the channel y = 0.5 for Re = 0.1 and W = 32 at a time instant of t = 2. A fountain flow is observed in the section under consideration. Three zones can be distinguished in the flow: a hydrodynamic flow stabilization zone near the inlet section; a fountain flow zone near the free surface; and a one-dimensional flow zone. The calculations showed that the lengths of the stabilization and fountain flow regions are less than unity for the selected values of the dimensionless criteria.
0 1 0 1 0 1
a b c
Fig. 6. Distribution of the kinematic characteristics along the channel for Re = 0.1, W = 32 at a time instant of t = 2 in a cross section of y = 0.5: (a) velocity Ux, (b) velocity Uz, and (c) pressure
In the other longitudinal sections, distributions of the kinematic characteristics qualitatively coincide with those presented in Fig. 6. In the section of x = y, velocity distributions (Ux and Uy) coincide with each other, which also confirms the efficiency of the computational algorithm.
The obtained results are compared with those of other authors. In particular, the results of studying of a channel filling in the creeping flow approximation are presented in [12]. Fig. 7 a shows the calculated values of the free surface height as a function of the parameter W presented in [13] (the solid line) and the results obtained by the VOF method (the dots).
Fig. 7 b illustrates a difference in the free surface convexity depending on the parameter W.
X X (x, 0.5)
Fig. 7. (a) Comparison of the dependence of free surface convexity on the parameter W in the current work (the dots) and in [13] (the solid line); (b) free surface shapes in a cross section of y = 0.5 for Re = 0.1 at the same time instant: 1 — W = 0, 2 — W = 10, and 3 - W = 100
The maximal difference in the calculated results does not exceed 3%. As a result of comparison, qualitative agreement and little quantitative deviations are observed.
Conclusion
In this work, testing of a modified VOF method on the three-dimensional fluid flow modeling is implemented. In particular, the PLIC-VOF algorithm, which is developed to solve the problems of two-dimensional fluid flows with a free surface, is generalized to the case of three-dimensional flows. Comparing with available data in literature and by testing approximation convergence are justified by efficiency of the developed algorithm. As a demonstration of the operation of the calculation program, the results of a study on filling a rectangular channel are given. Parametric calculations show that the free surface assumes a steady convex shape over time and then moves along the channel at a constant velocity.
The research is implemented at the expanses of the Russian Science Foundation (project no. 18-19-00021).
References
[1] I.A.Glushkov et al., Simulation of the molding of articles from free-molding compositions, Moscow, Arkhitektura-S, 2007 (in Russian).
[2] N.D.Katopodes, Free-Surface Flow: Computational Methods, 2019.
[3] B.D.Nichols, C.W.Hirt, R.S.Hotchkiss, SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries, Los Alamos Scientific Laboratory Report, LA-8355, 1980.
[4] W.Jang, J.Jilesen, F.S.Lien, H.Ji, A study on the extension of a VOF/PLIC based method to a curvilinear co-ordinate system, International Journal of Computational Fluid Dynamics, 22(2008), no. 4, 241-257.
[5] A.Issakhov, Y.Zhandaulet, A.Nogaeva, Numerical simulation of dam break flow for various forms of the obstacle by VOF method, International Journal of Multiphase Flow, 109(2018), 191-206. DOI: 10.1016/J.IJMULTIPHASEFLOW.2018.08.003
[6] D.M.Hargreaves, H.P.Morvan, N.G.Wright, Validation of the Volume of Fluid Method for Free Surface Calculation: The Broad-Crested Weir, Engineering Applications of Computational Fluid Mechanics, 1(2007), no. 2, 136-146. DOI: 10.1080/19942060.2007.11015188
[7] S.Hansch, D.Lucas, T.Hohne, E.Krepper, Gu.Montoya, Comparative Simulations of Free Surface Flows Using VOF-Methods and a New Approach for Multi-Scale Interfacial Structures, Proceedings of the ASME 2013 Fluids Engineering Summer Meeting (FEDSM2013-16104), 2013. DOI: 10.1115/FEDSM2013-16104
[8] S.Saincher, J.Banerjee, A Redistribution-Based Volume Preserving PLIC-VOF Technique, Numerical Heat Transfer, 67(2015), 338-362. DOI: 10.1080/10407790.2014.950078
[9] X.Yin, I.Zaricos, N.K.Karadimitriou, A.Raoof, S.M.Hassanizadeh, Direct simulations of two-phase flow experiments of different geometry complexities using Volume-of-Fluid (VOF) method, Chemical Engineering Science, 195(2018) 820-827.
DOI: 10.1016/j.ces.2018.10.029.
[10] M.Karimi, H.Droghetti, D.L.Marchisio, Multiscale Modeling of Expanding Polyurethane Foams via Computational Fluid Dynamics and Population Balance Equation, Macromol. Symp., 360(2016), 108-122. DOI: 10.1002/masy.201500108
[11] L.G.Loytsyansky, Fluid and Gas Mechanics, Moscow-Leningrad, Gostekhizdat, 1950.
[12] S.Patankar, Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, 1980.
[13] G.R.Shrager, A.N.Kozlobrodov, V.A.Yakutenok, Modeling of hydrodynamic processes in polymer processing technology, Tomsk, TGU Publ., 1999.
Моделирование пространственного заполнения емкости вязкой жидкостью с использованием VOF-метода
Евгений И. Борзенко Ефим И. Хегай
Томский государственный университет Томск, Российская Федерация
Аннотация. В настоящей работе представлены результаты моделирования пространственного течения ньютоновской жидкости со свободной поверхностью. Алгоритм РЫС УОР, предназначенный для решения задач о течении жидкостей со свободной поверхностью в двумерной постановке, обобщен на случай пространственных потоков. Работоспособность разработанного алгоритма и достоверность получаемых результатов продемонстрированы путем сравнения с литературными данными и проверкой аппроксимационной сходимости.
Параметрические расчеты заполнения канала с прямоугольным сечением показали, что с течением времени свободная поверхность принимает установившуюся выпуклую форму, которая перемещается вдоль канала с постоянной скоростью. В результате параметрического исследования построены зависимости геометрических характеристик формы свободной поверхности от параметров задачи.
Ключевые слова: течение ньютоновской жидкости, заполнение прямоугольного канала, свободная поверхность, трехмерное моделирование, численное моделирование, УОР-метод, структура потока.
