CC BY
31
Вычислительные технологии
Том 1, т 1, 1996
NUMERICAL SIMULATION ON SIMPLE SALINE OSCILLATORY FLOWS WITH FREE SURFACES*
YUKIMASA TAKEMOTO Faculty of Economics, Yokkaichi University (Yokkaichi, Mie Prefecture, 512 Japan)
Satoshi Chiba Ocean Engineering Research, Inc.
(Shinjuku-ku, Tokyo, 169 Japan)
Mina Okamura, Kenichi Yoshikawa Graduate School of Human Informatics, Nagoya University (Nagoya, Aichi Prefecture, 464-01 Japan)
In recent years, as the problems related to environmental conservation of water area and coastal development have become more serious and more complicated, it has become important to clarify the characteristics of the various kinds of density currents. A density current is the flow of one fluid within another caused by the density difference between the two fluids.
In this paper, we would like to report a curious self-oscillatory phenomenon, so called saline oscillator, induced by gravitational instability. An oscillatory flow is generated here, when a cup having small orifice in its base filled with saline water is placed within an outer vessel containing fresh water (Fig. 8). This system exhibits various nonlinear characteristics such as bifurcation of the oscillatory flow, oscillation cycle time and entrainment among the oscillations according to various experimental conditions.
In order to explain the simple saline oscillatory phenomenon, numerical simulation in two dimensional system has been performed. We have developed a simulation program for incompressible fluids that is extended to calculate the two layer time-dependent density currents. In this simulation program, a low diffusion, stable and high-accurate algorithm called Cubic-Interpolated pseudo-Particle (CIP) method for solving hyperbolic equations was adopted to calculate the Navier-Stokes equation and the equation of mass conservation. The experimental trend on the behavior of oscillatory flow has been reproduced in this simulation in a qualitative manner. The comparison between the results of its computation and experiment is fairly good (Fig. 9).
1. Introduction
In recent years, interface between convective flow field and multi-component diffusion under non-equilibrium condition is one of the important problems in environmental hydraulics and
*© Yukimasa Takemoto, Satoshi Chiba, Mina Okamura, Kenichi Yoshikawa, 1996.
chemical engineering.
We would like to report a curious self-oscillatory phenomenon, saline oscillator, induced by gravitational instability. The experimental procedure is quite simple and highly reproducible. An oscillatory flow of water generated here, when a cup with orifice (ca. 1~ 5 mm diam.) in its base, filled with saline water, is located within an outer vessel containing pure water (Fig. 8). A rhythmic change is then observed between the up-flow of the pure water and the down-flow of the saline water (Fig. 8). The period of oscillation can be several seconds or several minutes, depending on the experimental condition, i. e., the density of NaCl solution, the orifice size at the bottom of a cup, the surface area covered by fluid within a cup, and so on. As a typical example, the period was 13.6 seconds with the density: p0 = 1-0 g/cm3 and ps = 1.13 g/cm3, the dimensions: a = 0.55 mm, b = 21.0 mm, d = 1.0 mm, h = 20.0 mm.
In this saline oscillator, as the salt water going down to the outer vessel forms a layer on the bottom of the vessel and the pure water coming up to the cup forms same layer at the top of the saline water, the density difference near the orifice is kept essentially constant. Therefore, the rhythmic changes continue for several hours.
In this paper, simulation programs in 2D case was developed as a basis for computations of density currents driven by salinity, turbidity and temperature inhomogeneities. These simulation programs developed here are suitable for the incompressible fluid that are extended to calculate the two dimensional, two layer time-dependent density currents.
In these simulation programs, a low diffusion, stable and high-accurate algorithm called Cubic-Interpolated pseudo-Particle (CIP) method for solving hyperbolic equations was adopted to calculate the Navier—Stokes equations and the equation of mass conservation.
As numerical examples, a density current going forward to the fresh water tank and a densimetric (or lock) exchange flow were calculated. The results of the comparison between another computation and experiments were satisfactory. An oscillatory flow of water generated by the density difference between saline water and pure water when a cup with a small orifice in its base is placed within an outer vessel was also computed.
The experimental trend on the behavior of oscillatory flow has been reproduced in this simulation in a qualitative manner. The comparison between the results of its computation and experiment was fairly good (Fig. 9).
2. Basic Equations
We consider a two dimensional rectangular flow field with free surfaces. The flow is described by the Reynolds-averaged equations written in the following conservation form :
du dv
dX + dy = °’ (1)
du + dM + d(vu) = _ d(P/P) . d_ f2( . ) du\ +
dt dx dx dx dx \ V V dx)
d ( (du dv\\
+dy + Vt){dy + ax)) + g"
dv d(uv) d(vv) d(p/p) d ( . sdv\
at +J8T + V = ^ + dy (2(v + ay) +
+dx ((,/+v,)( £ + d^)) +gv ■ (2)
In these equations, u is the horizontal velocity and v is the vertical one, p is the pressure, p is the density, g = (gu,gv) is the gravitational acceleration, v is the molecular viscosity, vt is the turbulent viscosity.
It is necessary to have a equation for the density changes in order to compute the stratified density current. Usually, numerical calculations for incompressible density flow and buoyant jet have employed the Boussinesq approximation, but we directly calculate all the effects of density changes applying the following equation of mass conservation :
dp 3(pu) d(pv) _
dt dx dy
while subject to the total derivative Dp/Dt = 0, equation (3) is rewritten as the incompressibility condition V • V = 0. And so, we can use equation (3) with satisfying exactly equation
(1). In equation (3), it has been assumed that the diffusion of the solute is insignificant in the present computations, otherwise diffusion terms are needed in the right hand side of equation
(3).
3. Finite Difference Algorithm
In this paper, a accurate and stable numerical modeling which is called Cubic Interpolated pseudo-Particle (CIP) method for solving general hyperbolic equations is used. This CIP scheme was proposed by Yabe (1988, 1991), and it has good properties for solving a shock wave propagation and the convection diffusion equation in compressible and incompressible fluids. In incompressible flow computations, usually the evaluation of a higher-order formula (for example, QUICK scheme) involves more points. But the CIP scheme is less points and more accurate because it has the first special derivative in each grid point.
For example, we apply the two dimensional CIP scheme to the equation of mass conservation. At first, equation (3) can be rewritten as
pt + upx + vpy + pux + pvy = 0, (4)
by expanding derivatives. In this equation, we have another concise expression for partial differential form that describes dp/dt as pt. In the following equations the same expressions are adopted.
The placement of the field variables relative to the grid is shown in Fig. 1. We split the above equation into two phases; the Eulerian non-advection phase and the Lagrangian advection phase.
3.1 First Eulerian Phase
In this phase, the half part of equation (4)
pt = —pux — pvy (5)
is solved using the finite differencing:
pi,j = pi,j + At(-pi,j (ui+1/2,j — ui-1/2,j )/Ax — pi,j (vij+l/2 — vi,j-l/2)/Ay).
(6)
Figure 1. Field variables layout of the computational grid.
In the above equation, p* means a density value after one time step in the nonadvection phase. After p* is computed, the first special derivative p is approximated as the following diffrencing in x-direction:
(p'i,j - pi,j)/At = ((pi+i,j - pi-i,j)/2Ax - p+i,j - pi-i,j)/2Ax))/A-t. (7)
And then we obtain an equation, by eliminating At in equation (7):
p'i,j \x — pi,j \x + (p*+1,j — p*i-l,j - pi+1,j + pi-1,j)/2Ax,
(8)
where \x denotes the x-direction and it means p \x equals dxp. In a similar way, p for y-direction
(dyp) is
p'i,j \y = pi,j \y +(p*i,j+i — pi,j-i — pi,j+i + pi,j-i)/2AV- (9)
Before proceeding further, we take p* and p* values into new p and p values:
pLj \x^ pij \
'j lx
*
hj \y
/* I / I
pi,j \y * pi,j \y ■
j lx’ jy
3.2Second Lagrangian Phase
After p and p are advanced in the previous non-advection phase, the CIP solver is applied to the advection phase. In this phase, next two equations as
pt + upx + vpy — 0,
p't + upX + vpy = 0
(10)
(11)
are solved by shifting the cubic polynomial.
At first, we consider the x-directional components of the above equations. Let us begin with a linear hyperbolic equation:
pt + upx — °-
(12)
The solution of this equation is p(x,t) = p(x — ut, 0), if u is constant. However, we might use the solution for the time evolution of u during a very short period At, that is estimated as
p(x, t + At) — p(x — uAt, t).
(13)
In the point of xi — uAt, the conventional cubic spline interpolation function ^(x) is calculated between xi_i and xi as follows:
^i(x) = ai(x — x'i-i )3 + bi(x — xi_i)2 + Ci(x — xi_ i) + pi-i,
where
ai = (si — Si_i)/Ax,
bi 3si_i j
Ci = (pi — pi_i)/Ax — Ax(si + 2si_i),
Si = ^"(x)/6. (14)
In this case, the spline function requires the continuity of p, dp/dx = p', and d2p/dx2 = p". In contrast, we eliminate the continuity of d2p/dx2, then ci in the above equation would be p' in the following way:
Ci = (pi — pi_i)/Ax = pi_i.
Therefore substitution of ci into equation (14) gives
^i(x) = ai (x — xi_i)3 + bi (x — xi_i)2 + pi_i(x — xi_i) + pi_i. (15)
How the ai and bi values must be determined in this equation? We suppose that p' value is
also transported by a hyperbolic equation:
pt + upx = 0. (16)
The solution during a short period At can be approximated by equations (12) and (16), then the predicted profile will be:
Pi+l |x = ^i(xi — uAt) = ai(3 + bi(2 + pi_iC + pi_ij (17)
pT+1 \x = ^i(xi — uAt) = 3ai(2 + 2biC + pi_ij (18)
where
Z = (xi — uAt) — xi_i = Ax — uAt.
And so, we determine ai and bi in terms of p and p' because of requiring the continuity of ^i and ^i
ai = (pi + pi_i)/Ax<2 — 2(pi — pi_i)/Ax^j (19)
bi = — (pi + 2pi_i)/Ax + 3(pi — pi_i)/Ax. (20)
Once we obtain ai and bi in the above equations from previous pi and pi the values of and p'ii +i at the next time step are simply given by shifting the polynomial interpolation function as in equations (17) and (18).
These expressions are derived for u > 0. For u < 0, similar expressions in equations (19) and (20) by replacing i ^ i + 1 and i — 1 ^ i are obtained.
Secondly, we also consider the y-directional components of the equations (10) and (11). The similar equations can be derived:
pnj+l \y = ^j(Vj — vAt) = ajV3 + bjn2 + pj_iV + pj_ij (21)
pjn+i \y = ^j(yj — vAt) = 3ajn2 + 2bjn + pj_ij (22)
where
n = (y — vAt) — yj_i = AV — vAtj
where aj and bj would be derived using the same equation in (19) and (20), but these expressions can be omitted here.
Finally, the total value of two dimensional form becomes
j = pTi \x +p;+i \y. (23)
And the total p' values in each direction are calculated from the following equation:
p't + p'ux + p'vy = 0. (24)
The above equation is divided into two equations:
(p' \x)t = —pT+l \x ux — pT+i \x vy =
= —pT+1 \x ux — p'ji+1 \y vxj (25)
(p' \y)t = —pT+1 \y ux — pT+i \y vy =
= —pT+l \x uy — pj+1 \y vy. (26)
Consequently, final pj1 \x and pj1 \y values can be obtained.
These two phases complete one time step of the whole calculation. Then after replacing
pi +1, p'i +1 \ x and p'i +1 \ y by p, p' \ x and p' \ y, we return to the non-advection phase again.
The CIP solver is applied to the x- and y-momentum equations in the same way. But the explanation is too long, so we don’t describe it in this paper.
3.3Modified CIP2D Scheme
In the previous sections, calculated values of p and p are like a one-dimensional form. We use here a new multi-dimensional CIP method in the modified CIP scheme.
The quantities p, dxp(= p \x) and dyp(= p \y) are advanced in Lagrangian phase instead
of equations (17), (21), (23) and (18), (22) according to
p7+ i = aiC3 + biC2 + dxpi,j C + aj n3 + bj n2 + dy pi,j n+
+Ci,jC n + di,jCn + ei,jCn + pi,j, (27)
dxpj = 3aiC2 + 2biC + dxpi,j + Cij Cn + dij n2 + eid n, (28)
dy pi,j 3ajn + 2bjn + dy pi,j + Ci,jC + di,jCn + ei,jC. (29)
The above equations are derived for u < 0 and v < 0, and therefore C = —uAt, n = —vAt.
And above ai, ... , ei>j parameters are as follows ;
ai (dxpi,j + dxpimi,j)/Ax 2(pi,j pimi,j)/Ax j
bi = —(dxpi,j + 2dxpimi,j)/ Ax + 3(pi,j — pimi,j)/ Ax2 j aj (dy pi,j + dy pi,jmi)/Ay 2(pi,j pi,jmi) / Ay j
bj = —(dy pi,j + 2dy pi,jmi) / AV + 3(pi,j — pi,jmi) / Ay2 j
P pi,j pimi,j pi,jmi + pimi,jmij
Q dy pij — dyp
imijj
R dxpi,j dxpi,jmij
Ci,j = P/Ax2Ay — R/AxAy,
dij = P/AxAy2 — Q/AxAy,
ei,j = Q/Ax — P Ax Ay + R/ Ay,
where im1 = i — sgn(u), jm1 = j — sgn(v). For u < 0 and v < 0, then they become im1 = i + 1 and jm1 = j + 1.
4. Density Flow Computation
We are now applying the present scheme to several incompressible density flow problems. At first, a density current coming forward to the fresh water tank is computed to test the various kinds of finite difference schemes. In the present case, 30 x 12 grid points are used in x- and y-directions and the grid size is Ax = Ay = 1.0. The density of the saline water is pi = 1.03 and that of the fresh water is p2 = 1.0. These include the components of the gravity acceleration, gu = 0, gv = —980.0 cm/s2. A inlet boundary is in 3 computing cells with the left lower boundary wall and the saline water comes in the region at the inflow velocity of 10 cm/s. The results are shown in Fig. 2 as the marker particles profile. Figures 2(A) and (B) show the discontinuity of particles behind the top point of the current head by using a first order upwind method and a ZIP type central differencing scheme respectively. On the contrary. Figure 2(C) shows good distributions of particles by the CIP scheme. At the same situation, another computation is carried out using the CIP scheme by the modified CIP 2DKAI method in Fig. 3. As this CIP2DKAI code has no particles, the density difference of the two fluid is clarified by 3 contour lines at the value of p = 1.01, 1.015, 1.02.
Secondly, in order to analyze the applicability of the present method, a lock exchange flow shown on Fig. 4 is chosen as an example. Some selected time-marching results are represented in Fig. 4 and Fig. 5 with the aim of describing the behaviour of a densimetric exchange flow with the 1.03 : 1 density ratio. The simulation result of this particle code in Fig. 4 shows the irregularity of a free surface because of the surface cell pressures being defined as zero at the center of each surface cell. Figure 5 contrasts the surface smoothness using the modified surface treatment of the CIP2DKAI code. Figure 6 shows plots of the distance traveled by the density front as a function of time for several density ratio. In the present method, the computation is carried out for the Ap/p = 0.03 density ratio problem. The agreement between predictions of both present cases and another calculations performed by Komoda and Takeuchi is farely good. In addition to the prediction of Takeuchi, the measured data in the experiments by Takeuchi (1985) agree well with the results of its computation in Fig. 6.
Figure 7 contains slight variations of water level at both fresh water side wall and saline one. The surface level of the heavy fluid drops on the left wall and that of the lighter fluid rises on the right wall initially. In a short time interval, the free surface orientation at each wall has completed about one-half cycle of oscillation in Fig. 7. These variations of water level might be caused by the propagation and reflection of the surface wave as the density current moves by the gravitational acceleration. In the present calculation, a slight difference of a periodic time of these waves has been observed. The dashed line is proposed by Takeuchi at the same
condition of our computation. However, it is much difficult to explain the difference, and it is not possible to prove it by experiment.
Finally, a saline-water oscillator that was discovered by Martin (1970) and its theoretical and experimental report proposed by Yoshikawa (1991) is simulated by using the present scheme in Fig. 9. In this simulation, a outer vessel is filled with pure water and saline water is poured into a upper cup with a small orifice in its base. Then a rhythmic change is observed and simulated between the up-flow of fresh water and the down-flow of saline water. The period of oscillation can be several seconds depending on the experimental devices and simulation conditions in Fig. 8.
The simulation conditions are as follows. The diameter of orifice with a cup is 5 mm (a = 2.5 mm) and the diameter of a cup is 90 mm (b = 45 mm). The thickness of the base of a cup is 5 mm (d = 5 mm). Density values of NaCI solution and pure water are ps = 1.1 g/cm3 and p0 = 1.0 g/cm3 respectively. Initial depth of water is 8 cm, and computational region of this axisymmetric flow is a half size of total axisymmetrical experimental device. The period of this saline oscillator in these conditions is 1.05 seconds for the simulation, and 1.20 seconds for the experiment. The coincident between simulation and experiment is fairly good.
References
[1] Takemoto Y. AND Nakamura Y. A Three-Dimensional Incompressible Flow Solver. Lecture Notes in Physics, 264, 1986, 594-599.
[2] Takeuchi T. A Study on the Hydraulic Characteristics of the Head of Bottom-Layer Density Currents. Bull. Nat. Res. Inst. Fisheries Eng., t 6, 1985, 27-87.
[3] Takemoto Y., Tanaka M. and Nakamura Y. Numerical Simulation for Free Surface Flows in Generalized Coordinates. In “Extended Abstracts 4th Int. Conf. on Computing in Civil and Building Eng.", Tokyo, 1991, 198.
[4] Yabe T. AND Aoki T. A Universal Solver for Hyperbolic Equations by Cubic-polynomial Interpolation. Computer Physics Commun, 66, 1991, 219-242.
[5] Yoshikawa K. et al. Use of a Saline Oscillator as a Simple Nonlinear Dynamical System: Rhythms, Bifurcation, and Entrainment. Am. J. Phys. Teach., 59, 2, 1991, 137-141.
[6] Takemoto Y. and Torii K. Numerical simulation on unsteady density flows with free surfaces. Computational Mechanics’95. Proc. Int. Conf. on Comput. Eng. Sci., Aug. 1995, Hawaii, USA, 1, 1041-1046.
[7] Okamura M., Yoshikawa K. and Takemoto Y.. Simple Saline Oscillator as a System of Nonequilibrium Multi-Component Diffusion Coupled with Convective Flow. In “Extended Abstracts 5th Int. Symp. on BMED’95 and 6th Int. Conf. on MEBC95". Okinawa, Nov. 1995, 119-120.
Received for publication August 24, 1996
90 Yukimasa Takemoto, Satoshi Chiba, Mina Okamura, Kenichi Yoshikawa
Figure 2. Results of the particle simulation for the advancing heads of the density current. Density current simulation by using several kinds of differencing methods to convective terms: A: Upwind Method,
B: Zip Type Differencing,
C: CIP Method.
Figure 3. Results of the same simulation which includes no particles for the advancing heads of the density current using IP2DKAI method.
Figure 4. Densimetric exchange flow simulation with particle CIP method.
12.0 0.0 4.0 8.0 12.0 0.0 4.0 8.0 12.0 0.0 4.0 8.0 12.0 0.0 4.0 8.0 12.0 0.0 4.0 8.0 12.0
Figure 5. Densimetric exchange flow simulation using CIP2DKA method.
Figure 6. Front position versus time of the densimetric exchange flow.
Ah (cm)
I Ho = 10cm A,0 = 0.03
------- Present Computation (SOLADEN)
-------Simulation (Takeuchi)
0.2 -■ Ahz Left wall level
Figure 7. Time variation of a free surface level of the densimetric exchange flow.
Figure 8. Experimental device for simple saline oscillator.
Figure 9. Some comparison between experiment and simulation for sapine oscillator.
