- Email: [email protected]
Finite element analysis of bubbly flow around an oscillating square-section cylinder
Finite Elements in Analysis and Design 38 (2002) 803 – 821 www.elsevier.com/locate/ nel
Finite element analysis of bubbly #ow around an oscillating square-section cylinder Tomomi Uchiyama ∗ Center for Information Media Studies, Nagoya University Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
Abstract The bubbly #ow around a square-section cylinder, which is forced to oscillate normal to a uniform #ow, is simulated by an incompressible two-#uid model. The model is solved by the Arbitrary Lagrangian–Eulerian (ALE) nite element method. The width of the cylinder is 30 mm. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase upstream of the cylinder is 0.05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. With an increment in the cylinder oscillation frequency Sc, the vortex shedding frequency St approaches the Sc value. When the Sc is increased beyond a critical value, the lock-in occurs, in which the vortex shedding is synchronized with the cylinder oscillation. The time-averaged drag coe2cient and the root-mean-square of the lift coe2cient for the cylinder take their maximum values at the lock-in state. The volumetric fraction of the gas-phase g is higher in the wake of the cylinder, irrespective of the Sc value. The g value is lower below the lock-in state and higher in the lock-in state than the value for the stationary cylinder. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Multiphase #ow; ALE method; Forced oscillation; Lock-in; Phase distribution
1. Introduction In shell-and-tube type heat exchangers and steam generators, vibrations of the tube bundle are frequently caused by a gas–liquid two-phase #ow perpendicular to the tube axis. The vibrations result in tube failure or fretting-wear, depending on the circumstances. Therefore, experimental investigations on the #uid forces acting on the tube have been thus far carried out to make clear the mechanism of the vibrations [1,2], as reviewed by Nakamura et al. [3]. With the intention of clarifying such ∗
Corresponding author. Tel.:=fax: +81-52-789-5187. E-mail address: [email protected] (T. Uchiyama).
0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 1 ) 0 0 1 0 7 - X
804
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Nomenclature a B CD CL Cp f fc fv g n p Re Sc St t t∗ u uˆ x
cylinder oscillation amplitude width of square-section cylinder drag coe2cient of cylinder lift coe2cient of cylinder pressure coe2cient =(p − p0 )=( l ul20 =2) interfacial forces cylinder oscillation frequency vortex shedding frequency acceleration of gravity unit vector normal to boundary pressure Reynolds number = l l0 ul0 B=l dimensionless oscillation frequency =fc B=l0 ul0 Strouhal number =fv B=l0 ul0 time dimensionless time =l0 ul0 t=B velocity mesh velocity orthogonal coordinates volumetric fraction viscosity density viscous stress
Subscripts o g i; j l
upstream boundary or stationary cylinder gas-phase component in direction of xi or xj liquid-phase
#ow-induced vibration in simpler systems, Hara [4] observed the #ow around a single vibrating cylinder subjected to a two-phase cross-#ow, and made clear the vibration characteristics. Hara et al. [5] also investigated the vibration phenomena of circular cylinders in tandem, subjected to a bubbly cross-#ow, to examine the eKect of the pitch-to-diameter ratio. Joo and Dhir [6] measured the drag coe2cient of a vibrating circular cylinder under a two-phase cross-#ow, and they proposed a correlation formula for the drag coe2cient. To make clear the #ow-induced vibration of a bluK body under single-phase #ow conditions, the #ow around the body and the unsteady #uid forces on the body are investigated, in which the body is forced to oscillate either normal or parallel to the uniform #ow with speci ed amplitude and frequency [7]. It is reported that the vortex shedding is synchronized with the oscillation of the body when the body oscillates normal to the #ow at a frequency close to the vortex shedding frequency.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
805
This phenomenon is called the lock-in, and the mechanism is gradually being clari ed. But, such forced-oscillation method has not been applied to study the two-phase #ow-induced vibrations, and therefore the lock-in under the two-phase #ow condition has scarcely been investigated. Since a numerical method to solve moving boundary problems related with gas–liquid two-phase #ow has not presented, the author proposed a nite element method for two-phase #ow around a moving body in prior papers [8,9]. In the method, an incompressible two-#uid model is expressed in the ALE (Arbitrary Lagrangian–Eulerian) form, and the model is solved by the Galerkin nite element method. The method has been applied to calculate the bubbly #ow around a single hydrofoil under the heaving motion [8]. The #ow around a circular cylinder, which is forced to oscillate sinuously in a quiescent two-phase mixture, has also been simulated [9]. In this study, the bubbly #ow around a square-section cylinder, which is forced to oscillate normal to a uniform #ow, is simulated by the above-mentioned ALE nite element method. The lock-in phenomena under two-phase #ow condition are discussed in relation to the #ow and the #uid forces acting on the cylinder.
2. Governing equations 2.1. Moving boundary problem Fig. 1 illustrates a moving boundary problem. A solid body is moving with velocity in a #ow eld S. It is supposed that the boundary of S consists of a moving boundary C3 , which is the body surface, and xed boundaries C0 ; C1 and C2 . The domain S transforms in accordance with the body motion. In the ALE nite element method [8–10], the mesh is rearranged at every time step to discretize ˆ The ALE method calls uˆ the mesh the domain S and accordingly it moves with a velocity u. velocity. In the case of dividing S into quadrilateral elements as shown in Fig. 1, the mesh velocity uˆ and the #uid velocity u are de ned on the nodes. 2.2. Incompressible two-3uid model It is postulated that the gas- and liquid-phases are incompressible and that the mixture #ows isothermally without phase change. A one-pressure two-#uid model [11] gives the conservation equations of the mass and momentum for each phase: @k @ (k ukj ) = 0; + @t @xj
(k = g; l);
@uki 1 @p 1 @ @uki + ukj =− + (k kij ) + fki + gi ; @t @xj
k @xi k k @xj
(1) (k = g; l);
(2)
where g + l = 1:
(3)
806
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 1. Moving boundary problem.
Subscript k denotes the gas (k = g) or liquid (k = l) phase, and subscripts k, g, l disobey the summation convention in this paper. In Eq. (2), k is the viscous stress and fk is the interfacial forces. As the forces fk , an interfacial drag force fkiD , a virtual mass force fkiV , and a lift force fkiL are taken into account [12,13]. fkiD = −k (ugi − uli ); (k = g; l); @ugi @ugi @uli @uli + ug j fkiV = −k ; − − ul j @t @xj @t @xj fkiL = − k ijm (ugj − ulj )!m ;
(k = g; l);
(4) (k = g; l);
(5) (6)
where ijm and !m in Eq. (6) are the permutation tensor and the vorticity of the liquid-phase, respectively. The coe2cients k ; k and k are expressed as 1 a l CD 1 aCD |ug − ul |; l = − |ug − ul |; (7) g = 8 g g 8 l g = ( l = g )CV ;
l = −(g =l )CV ;
(8)
g = ( l = g )CL ;
l = −(g =l )CL :
(9)
Here, a is the interfacial area concentration, and CD ; CV and CL are the drag coe2cient, the virtual mass coe2cient and the lift coe2cient, respectively. The viscous stress kij is determined by the following equation [14,15]. @ukj @uki 2 @ukm ; (k = g; l): (10) + − ij kij = k @xj @xi 3 @xm
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
807
The time diKerentials with the spatial coordinate x held constant are involved in Eqs. (1), (2) and (5). When they are rewritten by considering the mesh velocity u, ˆ the ALE description of the incompressible two-#uid model is obtained [8,9]. @k @k @ (k ukj ) − uˆ j = 0; + @t @xj @xj
(k = g; l);
@uki @uki + (ukj − uˆ j ) = ! k "g i + # k " l i ; @t @xj
(11)
(k = g; l);
(12)
where 1 @p 1 @k kij + + fkiD + fkiL + gi ;
k @xi k k @xj −1 ! g #g 1 + g −g = : ! l #l l 1 − l
"ki = −
(k = g; l);
2.3. Numerical method [8; 9] The time diKerence equations for Eqs. (11) and (12) are written as kn+1 − kn @n @ n n+1 + (k ukj ) − uˆnj k = 0; Nt @xj @xj
(k = g; l);
ukin+1 − ukin @un n − uˆnj ) ki = !kn "gni + #nk "lni ; + (ukj Nt @xj
(k = g; l);
(13) (14)
where a superscript n stands for the condition at a time t = nNt, and "kin = −
1 @pn+1 1 @ n n n+1 + n ( ) − nk (ugn+1 i − ul i )
k @xi k k @xj k kij
− nk ijm (ugn j − ulnj )!mn + gi ;
(k = g; l):
Solving Eq. (14) for ukin+1 , the following is obtained. ukin+1 = − where 1 = bk
∗k n ki
Nt @pn+1 +
∗k @xi
!gn #ng +
g
l
n ki ;
= bk 'gni + ck 'lni ;
(k = g; l);
+ ck
!ln #nl +
g
l
(k = g; l);
(15)
;
(k = g; l);
808
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
ukin
n (ukj
@un uˆnj ) ki @xj
!n @ n n #n @ n n = − Nt − − nk (g gij ) + nk ( ) g g @xj l l @xj l lij +(!kn ng + #nk nl ) ijm (ugn j − ulnj )!mn − (!kn + #nk )gi ; (k = g; l); −1 1 + Nt(!gn ng + #ng nl ) −Nt(!gn ng + #ng nl ) b g cg = : bl cl Nt(!ln ng + #nl nl ) 1 − Nt(!ln ng + #nl nl ) 'kin
By taking the summation of both members of Eq. (13) for each phase and substituting Eq. (3) into the resulting equation, the conservation equation for the mass of the two-phase mixture is obtained. @ n n+1 ( u + ln uln+1 j ) = 0: @xj g gj
(16)
The #ow properties at a time t = (n + 1)Nt are calculated by solving Eqs. (13), (15) and (16) simultaneously. For the simultaneous calculations, the present method employs a fractional step method. In the rst step, the auxiliary velocity u˜ ki is estimated by neglecting the pressure gradient in Eq. (15). u˜ ki =
n ki ;
(k = g; l):
(17)
When Eq. (17) is subtracted from Eq. (15), the following is obtained: ukin+1 = u˜ ki −
1 @( ;
∗k @xi
(k = g; l);
(18)
where ( is a function satisfying pn+1 = (=Nt:
(19)
To calculate (, the following Poisson equation is derived by substituting Eq. (18) into Eq. (16). @ @xj
gn ln +
∗g
∗l
@( @ n = ( u˜ gj + ln u˜ lj ): @xj @xj g
(20)
In the second step, ukin+1 and pn+1 are calculated by substituting ( obtained from Eq. (20) into Eqs. (18) and (19), respectively. kn+1 is also calculated from Eq. (13) using ukin+1 . The abovementioned equations are solved by the Galerkin nite element method. Quadrilateral element with four nodes is employed. Fig. 2 illustrates the element. The pressure is de ned at the center of the element and assumed to be constant therein. The other variables, such as the velocities and the volumetric fractions, are de ned on the vertexes (nodes) of the element. Their values within
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
809
Fig. 2. Quadrilateral element.
the element are interpolated by linear functions. The mesh velocity is determined by the Euler’s backward scheme based on the mesh arrangements at t = nNt and (n − 1)Nt. 2.4. Boundary condition A square-section cylinder with width of B = 30 mm is forced to oscillate sinuously normal to the uniform bubbly #ow as shown in Fig. 3. The displacement of the cylinder x1 is given by the following equation: x1 = a sin(2)fc t);
(21)
where a is the oscillation amplitude and fc is the oscillation frequency. The boundary condition is listed in Table 1. The #uid velocities and the volumetric fractions are given on the inlet C0 , a traction free condition is prescribed on the outlet C1 , and a slip condition is imposed on the lateral boundary C2 . The mesh velocity uˆ is zero at these xed boundaries. On the cylinder surface, slip condition is assumed and uˆ is expressed as uˆ 1 = 2)afc cos(2)fc t);
uˆ 2 = 0:
(22)
2.5. Calculating condition Table 2 summarizes the calculating condition. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase at the inlet C0 ; g0 , is 0:05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. The dimensionless oscillation frequency of the cylinder Sc ranges from 0:05 to 0:16 at the oscillation amplitude a=B of 0:1. The #ow eld is divided into 47 × 68 meshes. The initial condition is given by the fully developed #ow for the stationary cylinder.
810
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 3. Calculating region.
Table 1 Boundary condition Inlet boundary C0 Outlet boundary C1 Lateral boundary C2 Cylinder surface C3
uki = uS ki ; ki = Ski ; uˆ i = 0 nj ( ij p − kij ) = 0; uˆ i = 0 uk1 = k2 = 0; uˆ i = 0 uk1 = uˆ 1 ; k2 = 0 on lateral surfaces uk2 = k1 = 0 on front and rear surfaces
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
811
Table 2 Calculating condition Reynolds number: Re Volumetric fraction of gas phase: g0 Oscillation amplitude: a=B Oscillation frequency: Sc
200 0.05 0.1 0:05–0:16
3. Numerical results and discussion 3.1. Unsteady change in 3uid forces When the drag and lift forces acting on the cylinder are estimated from the pressure distribution on the cylinder, the drag coe2cient CD and the lift coe2cient CL change as functions of the dimensionless time t ∗ as shown in Fig. 4. The result for the stationary cylinder is plotted in Fig. 4(a). The coe2cients CD and CL change periodically with constant periods, suggesting the periodical shed of the vortices from the cylinder.
Fig. 4. Time-variation of CD and CL .
812
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 5. Spectra of CL .
The amplitude of CL is not constant, but the change in CL is found to be superimposed on a small oscillation with a longer period. Fig. 4(b) shows the result at the oscillation frequency of Sc = 0:1. The amplitudes of CD and CL vary greatly. Large spikes are intermittently observed in the change of CD . The period of CL diKers from that of the dimensionless displacement of the cylinder x1 =B plotted by a broken line. This indicates that vortices are shed from the cylinder independent of the cylinder oscillation. The result at Sc = 0:15 is shown in Fig. 4(c). At t ∗ ¿ 14, the amplitude and period of CD and CL are almost constant, and the period of CL coincides with that of the cylinder oscillation. This suggests the appearance of the lock-in, in which the vortex shedding is synchronized with the cylinder oscillation. The amplitude of CD is much smaller than those in Figs. 4(a) and (b). The period of CD is half of that for CL , as in the case for the stationary cylinder. 3.2. Lock-in phenomena Calculating the spectra of CL by FFT, they distribute against the dimensionless frequency as shown in Fig. 5.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
813
Fig. 6. Lock-in region for Sc.
The result for the stationary cylinder is plotted in Fig. 5(a). A peak exists at a frequency of 0.129. Thus, the Strouhal number St, which is the dimensionless frequency of the vortex shedding, is found to be 0:129. Fig. 5(b) shows the result at the cylinder oscillation frequency of Sc = 0:1. Two large peaks are observed at frequencies of 0.1 and 0.135. These values correspond to the cylinder oscillation (Sc = 0:1) and the vortex shedding (St = 0:135). The St value is slightly higher than that for the stationary cylinder calculated in Fig. 5(a). The result at Sc = 0:15 is shown in Fig. 5(c). The existence of a large peak at a frequency of 0:15 indicates the occurrence of the lock-in, in which the vortex shedding frequency St coincides with the cylinder oscillation frequency Sc. To examine the lock-in region for Sc, the ratio St=Sc is plotted against Sc by the symbol o in Fig. 6. The ratio St=Sc decreases with an increment in Sc at 0:05 6 Sc 6 0:13. At Sc ¿ 0:13, it is unity, and hence the lock-in occurs. The critical value of Sc for the lock-in occurrence, (Sc)cr = 0:13, is nearly the same as the Strouhal number (St = 0:129) for the stationary cylinder. The change in St=Sc under the single-phase #ow condition (g0 = 0) is also plotted by the symbol x in Fig. 6. It is almost parallel with the result at g0 = 0:05, though the (Sc)cr value slightly lowers. When the time-averaged value of CD ; CD , and the root-mean-square of the #uctuation for CL ; CL , are calculated, they vary against Sc as presented in Fig. 7. The ratios, CD =CD0 and CL =CL0 , are plotted by the symbol o, where CD0 and CL0 are the values for the stationary cylinder. The ratio CL =CL0 is lower than unity below the lock-in state (Sc ¡ 0:13). But, it is larger than unity in the lock-in state (0:13 6 Sc ¡ 0:16), and it reaches its maximum value at Sc = 0:15. The CL value is proportional to the intensity of the vortex shed from the cylinder. The vortex, synchronized with the cylinder oscillation, is given high energy from the cylinder, resulting in an increment in the vortex intensity, as discussed in Fig. 13. Therefore, CL takes its maximum value at Sc = 0:15. The ratio CD =CD0 reaches its maximum value at Sc = 0:13. The Sc value coincides with the (Sc)cr value for the lock-in occurrence. The time-averaged distribution of the pressure on the cylinder is shown in Fig. 8. The pressure coe2cient Cp , based on the pressure p0 upstream of the cylinder, is plotted against the distance from the center of the front face of the cylinder. Though the distribution on the front (A–B) is hardly
814
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 7. Dependance of CD and CL on Sc.
Fig. 8. Pressure distribution on cylinder.
aKected by the Sc value, the cylinder oscillation at Sc = 0:13 yields the lowest pressure on the rear face (C–D). The time-averaged drag coe2cient CD depends on the pressure diKerence between the front and rear faces. Therefore, CD takes its maximum value at Sc = 0:13, as seen in Fig. 7.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
815
Fig. 9. Time-variation of volumetric velocity of liquid-phase l ul at Sc = 0:1 (asynchronous state).
3.3. Time-variation of 3ow properties The time-variations of the volumetric velocity of the liquid-phase l ul around the cylinder are presented in Figs. 9 and 10. The instantaneous distributions at four time points during one period of the cylinder oscillation are displayed. Figs. 9(a,c) and 10(a,c) are the results when the cylinder passes right and left, respectively, through the oscillation center with its maximum speed, Figs. 9(b,d) and 10(b,d) are the results when the cylinder is at rest. When the cylinder oscillation frequency is below the lock-in region, l ul distributes as shown in Fig. 9, where Sc = 0:1. Comparing Figs. 9(a) and (c), which are the distributions when the cylinder passes through the oscillation center, the #ow similarity is not observed. The vortices are shed independent of the cylinder oscillation. Fig. 10 shows the distribution at the lock-in state (Sc = 0:15). Comparing the #ow elds in Figs. 10(a) and (c), large vortices are produced just behind the cylinder, and the #ows are symmetrical with respect to the center line of the cylinder. This indicates the synchronization of the cylinder oscillation and the vortex shedding. In addition, the vortices behind the cylinder are closer to the cylinder than
816
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 10. Time-variation of volumetric velocity of liquid-phase l ul at Sc = 0:15 (synchronous state).
those in Fig. 9. This upstream migration of the vortex is also reported under single-phase #ow condition [16]. The distribution of the pressure coe2cient Cp changes with the lapse of time as shown in Figs. 11–13. Fig. 11 shows the distribution for the stationary cylinder, where the results at four time points during one period of the vortex shedding are presented. The position where Cp has the minimum value corresponds to the center of the vortex. The minimum value increases in the downstream direction, because the vortex intensity is diminished due to the #uid viscosity. The pressure is higher locally at the front of the cylinder. Figs. 12 and 13 show the distributions corresponding to the results in Figs. 9 and 10, respectively. The instantaneous distributions at four time points during one period of the cylinder oscillation are displayed. Figs. 12(a,c) and 13(a,c) are the results when the cylinder passes right and left, respectively, through the oscillation center with its maximum speed, Figs. 12(b,d) and 13(b,d) are the results when the cylinder is at rest. Fig. 12 shows the distribution when the cylinder oscillation frequency is below the lock-in region. The Sc value is 0:1. The distribution upstream of the cylinder is almost the same as that for the stationary cylinder (Fig. 11). But the Cp value heightens slightly in the wake, demonstrating the
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 11. Time-variation of pressure coe2cient Cp for stationary cylinder.
Fig. 12. Time-variation of pressure coe2cient Cp at Sc = 0:1 (asynchronous state).
817
818
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 13. Time-variation of pressure coe2cient Cp at Sc = 0:15 (synchronous state).
decrement in the vortex intensity. Therefore, the CL value at Sc = 0:1 is smaller than that for the stationary cylinder, as found in Fig. 7. The distribution at the lock-in state is shown in Fig. 13, where Sc = 0:15. When comparing with the results for the stationary cylinder (Fig. 11), the Cp value lessens in the wake and such lower Cp region spreads in the downstream direction. Because the vortex intensity is heightened by the synchronization of the vortex shedding with the cylinder oscillation. The time-variations of the volumetric fraction of the gas-phase g , corresponding to the results in Figs. 11–13, are displayed in Figs. 14 –16, respectively. Fig. 14 shows the distribution for the stationary cylinder. Since the gas-phase tends to swarm in the vortex due to the pressure gradient, g is higher in the wake [14,15]. Such region with higher g value spreads downstream of the cylinder. g is lower near the front of the cylinder because there is a positive pressure gradient toward the front. Fig. 15 shows the result at Sc = 0:1. The distribution upstream of the cylinder is almost the same as that for the stationary cylinder (Fig. 14). But, the g value in the wake is not so high. This is because the gas-phase is less able to swarm behind the cylinder owing to the gentler pressure gradient. The distribution at the lock-in state (Sc = 0:15) is shown in Fig. 16. In comparison with the result for the stationary cylinder (Fig. 14), the g value is higher in the wake and such higher region is wider downstream of the cylinder. Because the vortex intensity is higher, and accordingly the gas-phase shifts well toward such lower pressure region.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 14. Time-variation of volumetric fraction of gas-phase g for stationary cylinder.
Fig. 15. Time-variation of volumetric fraction of gas-phase g at Sc = 0:1 (asynchronous state).
819
820
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 16. Time-variation of volumetric fraction of gas-phase g at Sc = 0:15 (synchronous state).
4. Summary The bubbly #ow around a square-section cylinder, which is forced to oscillate sinuously normal to the uniform #ow, is calculated by using an incompressible two-#uid model. The width of cylinder is 30 mm. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase upstream of the cylinder is 0.05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. The numerical results are summarized as follows: (1) With an increment in the cylinder oscillation frequency Sc, the vortex shedding frequency St approaches the Sc value. When Sc is increased beyond a critical value, the lock-in occurs, in which the vortex shedding is synchronized with the cylinder oscillation. The critical value is nearly the same as the St value for the stationary cylinder. (2) The time-averaged drag coe2cient and the root-mean-square of the lift coe2cient for the cylinder take their maximum values at the lock-in state. (3) The volumetric fraction of the gas-phase g is higher in the wake of the cylinder, irrespective of the Sc value. The g value is lower below the lock-in state and higher in the lock-in state than the value for the stationary cylinder. References [1] M.J. Pettigrew, D.J. Gorman, Experimental studies on #ow-induced vibration to support steam generator design, (1977) AECL-5804.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
821
[2] F. Axisa, B. Villard, R.J. Gibert, G. Hetsroni, P. Sundheimer, Vibration of tube bundles subjected to air–water and steam–water cross-#ow, Proceedings of the ASME Winter Annual Meeting, 1984, pp. 269 –283. [3] T. Nakamura, K. Fujita, A. Tsuge, Two-phase cross-#ow-induced vibration of tube arrays (Brief review of previous studies and summary of design methods), JSME Int. J. Ser. B 36 (3) (1993) 429–438. [4] F. Hara, Air-bubble eKects on vortex-induced vibrations of a circular cylinder, Proc. ASME Symp. FIV 1 (1984) 103–113. [5] F. Hara, T. Iijima, T. Nojima, Vibration response characteristics of tandem two circular cylinders subjected to a two phase cross-#ow — eKects of pitch-to-diameter ratio, Trans. ASME, J. Pressure Vessel Tech. 114 (4) (1992) 444–452. [6] Y. Joo, V.K. Dhir, An experimental study of drag on a single tube and on a tube in an array under two-phase cross-#ow, Int. J. Multiphase Flow 20 (6) (1994) 1009–1019. [7] O.M. Gri2n, M.S. Hall, Review — vortex shedding lock-on and #ow control in bluK body wakes, Trans. ASME, J. Fluid Eng. 113 (4) (1991) 526–537. [8] T. Uchiyama, K. Minemura, Numerical simulation of bubbly #ow around an oscillating hydrofoil, Proceedings of the Asian Symposium on Multiphase Flow, Osaka, 1999, pp. 27–32. [9] T. Uchiyama, ALE nite element method for gas–liquid two-phase #ow including moving boundary based on an incompressible two-#uid model, Nucl. Eng. Des. 205 (2001) 69–82. [10] T. Belytschko, D.P. Flanagan, J.M. Kennedy, Finite element methods with user-controlled meshes for #uid-structure interaction, Comput. Methods. Appl. Mech. Eng. 33 (1982) 669–688. [11] M. Ishii, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975. [12] T. Uchiyama, Numerical analysis of air–water two-phase #ow across a staggered tube bundle using an incompressible two-#uid model, Nucl. Sci. Eng. 134 (3) (2000) 281–292. [13] T. Uchiyama, Finite element analysis of gas–liquid two-phase #ow across an in-line tube bundle, Hybrid Method Eng. 2(4) (2001), accepted and in press. [14] T. Uchiyama, Numerical simulation of gas–liquid two-phase #ow around a rectangular cylinder by the incompressible two-#uid model, Nucl. Sci. Eng. 133 (1) (1999) 92–105. [15] T. Uchiyama, Petrov–Galerkin nite element method for gas–liquid two-phase #ow based on an incompressible two-#uid model, Nucl. Eng. Des. 193 (1999) 145–157. [16] P.W. Bearman, E.D. Obasaju, An experimental study of pressure #uctuatiions on xed and oscillating square-section cylinders, J. Fluid Mech. 119 (1982) 297–321.
Finite element analysis of bubbly #ow around an oscillating square-section cylinder Tomomi Uchiyama ∗ Center for Information Media Studies, Nagoya University Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
Abstract The bubbly #ow around a square-section cylinder, which is forced to oscillate normal to a uniform #ow, is simulated by an incompressible two-#uid model. The model is solved by the Arbitrary Lagrangian–Eulerian (ALE) nite element method. The width of the cylinder is 30 mm. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase upstream of the cylinder is 0.05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. With an increment in the cylinder oscillation frequency Sc, the vortex shedding frequency St approaches the Sc value. When the Sc is increased beyond a critical value, the lock-in occurs, in which the vortex shedding is synchronized with the cylinder oscillation. The time-averaged drag coe2cient and the root-mean-square of the lift coe2cient for the cylinder take their maximum values at the lock-in state. The volumetric fraction of the gas-phase g is higher in the wake of the cylinder, irrespective of the Sc value. The g value is lower below the lock-in state and higher in the lock-in state than the value for the stationary cylinder. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Multiphase #ow; ALE method; Forced oscillation; Lock-in; Phase distribution
1. Introduction In shell-and-tube type heat exchangers and steam generators, vibrations of the tube bundle are frequently caused by a gas–liquid two-phase #ow perpendicular to the tube axis. The vibrations result in tube failure or fretting-wear, depending on the circumstances. Therefore, experimental investigations on the #uid forces acting on the tube have been thus far carried out to make clear the mechanism of the vibrations [1,2], as reviewed by Nakamura et al. [3]. With the intention of clarifying such ∗
Corresponding author. Tel.:=fax: +81-52-789-5187. E-mail address: [email protected] (T. Uchiyama).
0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 1 ) 0 0 1 0 7 - X
804
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Nomenclature a B CD CL Cp f fc fv g n p Re Sc St t t∗ u uˆ x
cylinder oscillation amplitude width of square-section cylinder drag coe2cient of cylinder lift coe2cient of cylinder pressure coe2cient =(p − p0 )=( l ul20 =2) interfacial forces cylinder oscillation frequency vortex shedding frequency acceleration of gravity unit vector normal to boundary pressure Reynolds number = l l0 ul0 B=l dimensionless oscillation frequency =fc B=l0 ul0 Strouhal number =fv B=l0 ul0 time dimensionless time =l0 ul0 t=B velocity mesh velocity orthogonal coordinates volumetric fraction viscosity density viscous stress
Subscripts o g i; j l
upstream boundary or stationary cylinder gas-phase component in direction of xi or xj liquid-phase
#ow-induced vibration in simpler systems, Hara [4] observed the #ow around a single vibrating cylinder subjected to a two-phase cross-#ow, and made clear the vibration characteristics. Hara et al. [5] also investigated the vibration phenomena of circular cylinders in tandem, subjected to a bubbly cross-#ow, to examine the eKect of the pitch-to-diameter ratio. Joo and Dhir [6] measured the drag coe2cient of a vibrating circular cylinder under a two-phase cross-#ow, and they proposed a correlation formula for the drag coe2cient. To make clear the #ow-induced vibration of a bluK body under single-phase #ow conditions, the #ow around the body and the unsteady #uid forces on the body are investigated, in which the body is forced to oscillate either normal or parallel to the uniform #ow with speci ed amplitude and frequency [7]. It is reported that the vortex shedding is synchronized with the oscillation of the body when the body oscillates normal to the #ow at a frequency close to the vortex shedding frequency.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
805
This phenomenon is called the lock-in, and the mechanism is gradually being clari ed. But, such forced-oscillation method has not been applied to study the two-phase #ow-induced vibrations, and therefore the lock-in under the two-phase #ow condition has scarcely been investigated. Since a numerical method to solve moving boundary problems related with gas–liquid two-phase #ow has not presented, the author proposed a nite element method for two-phase #ow around a moving body in prior papers [8,9]. In the method, an incompressible two-#uid model is expressed in the ALE (Arbitrary Lagrangian–Eulerian) form, and the model is solved by the Galerkin nite element method. The method has been applied to calculate the bubbly #ow around a single hydrofoil under the heaving motion [8]. The #ow around a circular cylinder, which is forced to oscillate sinuously in a quiescent two-phase mixture, has also been simulated [9]. In this study, the bubbly #ow around a square-section cylinder, which is forced to oscillate normal to a uniform #ow, is simulated by the above-mentioned ALE nite element method. The lock-in phenomena under two-phase #ow condition are discussed in relation to the #ow and the #uid forces acting on the cylinder.
2. Governing equations 2.1. Moving boundary problem Fig. 1 illustrates a moving boundary problem. A solid body is moving with velocity in a #ow eld S. It is supposed that the boundary of S consists of a moving boundary C3 , which is the body surface, and xed boundaries C0 ; C1 and C2 . The domain S transforms in accordance with the body motion. In the ALE nite element method [8–10], the mesh is rearranged at every time step to discretize ˆ The ALE method calls uˆ the mesh the domain S and accordingly it moves with a velocity u. velocity. In the case of dividing S into quadrilateral elements as shown in Fig. 1, the mesh velocity uˆ and the #uid velocity u are de ned on the nodes. 2.2. Incompressible two-3uid model It is postulated that the gas- and liquid-phases are incompressible and that the mixture #ows isothermally without phase change. A one-pressure two-#uid model [11] gives the conservation equations of the mass and momentum for each phase: @k @ (k ukj ) = 0; + @t @xj
(k = g; l);
@uki 1 @p 1 @ @uki + ukj =− + (k kij ) + fki + gi ; @t @xj
k @xi k k @xj
(1) (k = g; l);
(2)
where g + l = 1:
(3)
806
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 1. Moving boundary problem.
Subscript k denotes the gas (k = g) or liquid (k = l) phase, and subscripts k, g, l disobey the summation convention in this paper. In Eq. (2), k is the viscous stress and fk is the interfacial forces. As the forces fk , an interfacial drag force fkiD , a virtual mass force fkiV , and a lift force fkiL are taken into account [12,13]. fkiD = −k (ugi − uli ); (k = g; l); @ugi @ugi @uli @uli + ug j fkiV = −k ; − − ul j @t @xj @t @xj fkiL = − k ijm (ugj − ulj )!m ;
(k = g; l);
(4) (k = g; l);
(5) (6)
where ijm and !m in Eq. (6) are the permutation tensor and the vorticity of the liquid-phase, respectively. The coe2cients k ; k and k are expressed as 1 a l CD 1 aCD |ug − ul |; l = − |ug − ul |; (7) g = 8 g g 8 l g = ( l = g )CV ;
l = −(g =l )CV ;
(8)
g = ( l = g )CL ;
l = −(g =l )CL :
(9)
Here, a is the interfacial area concentration, and CD ; CV and CL are the drag coe2cient, the virtual mass coe2cient and the lift coe2cient, respectively. The viscous stress kij is determined by the following equation [14,15]. @ukj @uki 2 @ukm ; (k = g; l): (10) + − ij kij = k @xj @xi 3 @xm
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
807
The time diKerentials with the spatial coordinate x held constant are involved in Eqs. (1), (2) and (5). When they are rewritten by considering the mesh velocity u, ˆ the ALE description of the incompressible two-#uid model is obtained [8,9]. @k @k @ (k ukj ) − uˆ j = 0; + @t @xj @xj
(k = g; l);
@uki @uki + (ukj − uˆ j ) = ! k "g i + # k " l i ; @t @xj
(11)
(k = g; l);
(12)
where 1 @p 1 @k kij + + fkiD + fkiL + gi ;
k @xi k k @xj −1 ! g #g 1 + g −g = : ! l #l l 1 − l
"ki = −
(k = g; l);
2.3. Numerical method [8; 9] The time diKerence equations for Eqs. (11) and (12) are written as kn+1 − kn @n @ n n+1 + (k ukj ) − uˆnj k = 0; Nt @xj @xj
(k = g; l);
ukin+1 − ukin @un n − uˆnj ) ki = !kn "gni + #nk "lni ; + (ukj Nt @xj
(k = g; l);
(13) (14)
where a superscript n stands for the condition at a time t = nNt, and "kin = −
1 @pn+1 1 @ n n n+1 + n ( ) − nk (ugn+1 i − ul i )
k @xi k k @xj k kij
− nk ijm (ugn j − ulnj )!mn + gi ;
(k = g; l):
Solving Eq. (14) for ukin+1 , the following is obtained. ukin+1 = − where 1 = bk
∗k n ki
Nt @pn+1 +
∗k @xi
!gn #ng +
g
l
n ki ;
= bk 'gni + ck 'lni ;
(k = g; l);
+ ck
!ln #nl +
g
l
(k = g; l);
(15)
;
(k = g; l);
808
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
ukin
n (ukj
@un uˆnj ) ki @xj
!n @ n n #n @ n n = − Nt − − nk (g gij ) + nk ( ) g g @xj l l @xj l lij +(!kn ng + #nk nl ) ijm (ugn j − ulnj )!mn − (!kn + #nk )gi ; (k = g; l); −1 1 + Nt(!gn ng + #ng nl ) −Nt(!gn ng + #ng nl ) b g cg = : bl cl Nt(!ln ng + #nl nl ) 1 − Nt(!ln ng + #nl nl ) 'kin
By taking the summation of both members of Eq. (13) for each phase and substituting Eq. (3) into the resulting equation, the conservation equation for the mass of the two-phase mixture is obtained. @ n n+1 ( u + ln uln+1 j ) = 0: @xj g gj
(16)
The #ow properties at a time t = (n + 1)Nt are calculated by solving Eqs. (13), (15) and (16) simultaneously. For the simultaneous calculations, the present method employs a fractional step method. In the rst step, the auxiliary velocity u˜ ki is estimated by neglecting the pressure gradient in Eq. (15). u˜ ki =
n ki ;
(k = g; l):
(17)
When Eq. (17) is subtracted from Eq. (15), the following is obtained: ukin+1 = u˜ ki −
1 @( ;
∗k @xi
(k = g; l);
(18)
where ( is a function satisfying pn+1 = (=Nt:
(19)
To calculate (, the following Poisson equation is derived by substituting Eq. (18) into Eq. (16). @ @xj
gn ln +
∗g
∗l
@( @ n = ( u˜ gj + ln u˜ lj ): @xj @xj g
(20)
In the second step, ukin+1 and pn+1 are calculated by substituting ( obtained from Eq. (20) into Eqs. (18) and (19), respectively. kn+1 is also calculated from Eq. (13) using ukin+1 . The abovementioned equations are solved by the Galerkin nite element method. Quadrilateral element with four nodes is employed. Fig. 2 illustrates the element. The pressure is de ned at the center of the element and assumed to be constant therein. The other variables, such as the velocities and the volumetric fractions, are de ned on the vertexes (nodes) of the element. Their values within
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
809
Fig. 2. Quadrilateral element.
the element are interpolated by linear functions. The mesh velocity is determined by the Euler’s backward scheme based on the mesh arrangements at t = nNt and (n − 1)Nt. 2.4. Boundary condition A square-section cylinder with width of B = 30 mm is forced to oscillate sinuously normal to the uniform bubbly #ow as shown in Fig. 3. The displacement of the cylinder x1 is given by the following equation: x1 = a sin(2)fc t);
(21)
where a is the oscillation amplitude and fc is the oscillation frequency. The boundary condition is listed in Table 1. The #uid velocities and the volumetric fractions are given on the inlet C0 , a traction free condition is prescribed on the outlet C1 , and a slip condition is imposed on the lateral boundary C2 . The mesh velocity uˆ is zero at these xed boundaries. On the cylinder surface, slip condition is assumed and uˆ is expressed as uˆ 1 = 2)afc cos(2)fc t);
uˆ 2 = 0:
(22)
2.5. Calculating condition Table 2 summarizes the calculating condition. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase at the inlet C0 ; g0 , is 0:05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. The dimensionless oscillation frequency of the cylinder Sc ranges from 0:05 to 0:16 at the oscillation amplitude a=B of 0:1. The #ow eld is divided into 47 × 68 meshes. The initial condition is given by the fully developed #ow for the stationary cylinder.
810
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 3. Calculating region.
Table 1 Boundary condition Inlet boundary C0 Outlet boundary C1 Lateral boundary C2 Cylinder surface C3
uki = uS ki ; ki = Ski ; uˆ i = 0 nj ( ij p − kij ) = 0; uˆ i = 0 uk1 = k2 = 0; uˆ i = 0 uk1 = uˆ 1 ; k2 = 0 on lateral surfaces uk2 = k1 = 0 on front and rear surfaces
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
811
Table 2 Calculating condition Reynolds number: Re Volumetric fraction of gas phase: g0 Oscillation amplitude: a=B Oscillation frequency: Sc
200 0.05 0.1 0:05–0:16
3. Numerical results and discussion 3.1. Unsteady change in 3uid forces When the drag and lift forces acting on the cylinder are estimated from the pressure distribution on the cylinder, the drag coe2cient CD and the lift coe2cient CL change as functions of the dimensionless time t ∗ as shown in Fig. 4. The result for the stationary cylinder is plotted in Fig. 4(a). The coe2cients CD and CL change periodically with constant periods, suggesting the periodical shed of the vortices from the cylinder.
Fig. 4. Time-variation of CD and CL .
812
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 5. Spectra of CL .
The amplitude of CL is not constant, but the change in CL is found to be superimposed on a small oscillation with a longer period. Fig. 4(b) shows the result at the oscillation frequency of Sc = 0:1. The amplitudes of CD and CL vary greatly. Large spikes are intermittently observed in the change of CD . The period of CL diKers from that of the dimensionless displacement of the cylinder x1 =B plotted by a broken line. This indicates that vortices are shed from the cylinder independent of the cylinder oscillation. The result at Sc = 0:15 is shown in Fig. 4(c). At t ∗ ¿ 14, the amplitude and period of CD and CL are almost constant, and the period of CL coincides with that of the cylinder oscillation. This suggests the appearance of the lock-in, in which the vortex shedding is synchronized with the cylinder oscillation. The amplitude of CD is much smaller than those in Figs. 4(a) and (b). The period of CD is half of that for CL , as in the case for the stationary cylinder. 3.2. Lock-in phenomena Calculating the spectra of CL by FFT, they distribute against the dimensionless frequency as shown in Fig. 5.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
813
Fig. 6. Lock-in region for Sc.
The result for the stationary cylinder is plotted in Fig. 5(a). A peak exists at a frequency of 0.129. Thus, the Strouhal number St, which is the dimensionless frequency of the vortex shedding, is found to be 0:129. Fig. 5(b) shows the result at the cylinder oscillation frequency of Sc = 0:1. Two large peaks are observed at frequencies of 0.1 and 0.135. These values correspond to the cylinder oscillation (Sc = 0:1) and the vortex shedding (St = 0:135). The St value is slightly higher than that for the stationary cylinder calculated in Fig. 5(a). The result at Sc = 0:15 is shown in Fig. 5(c). The existence of a large peak at a frequency of 0:15 indicates the occurrence of the lock-in, in which the vortex shedding frequency St coincides with the cylinder oscillation frequency Sc. To examine the lock-in region for Sc, the ratio St=Sc is plotted against Sc by the symbol o in Fig. 6. The ratio St=Sc decreases with an increment in Sc at 0:05 6 Sc 6 0:13. At Sc ¿ 0:13, it is unity, and hence the lock-in occurs. The critical value of Sc for the lock-in occurrence, (Sc)cr = 0:13, is nearly the same as the Strouhal number (St = 0:129) for the stationary cylinder. The change in St=Sc under the single-phase #ow condition (g0 = 0) is also plotted by the symbol x in Fig. 6. It is almost parallel with the result at g0 = 0:05, though the (Sc)cr value slightly lowers. When the time-averaged value of CD ; CD , and the root-mean-square of the #uctuation for CL ; CL , are calculated, they vary against Sc as presented in Fig. 7. The ratios, CD =CD0 and CL =CL0 , are plotted by the symbol o, where CD0 and CL0 are the values for the stationary cylinder. The ratio CL =CL0 is lower than unity below the lock-in state (Sc ¡ 0:13). But, it is larger than unity in the lock-in state (0:13 6 Sc ¡ 0:16), and it reaches its maximum value at Sc = 0:15. The CL value is proportional to the intensity of the vortex shed from the cylinder. The vortex, synchronized with the cylinder oscillation, is given high energy from the cylinder, resulting in an increment in the vortex intensity, as discussed in Fig. 13. Therefore, CL takes its maximum value at Sc = 0:15. The ratio CD =CD0 reaches its maximum value at Sc = 0:13. The Sc value coincides with the (Sc)cr value for the lock-in occurrence. The time-averaged distribution of the pressure on the cylinder is shown in Fig. 8. The pressure coe2cient Cp , based on the pressure p0 upstream of the cylinder, is plotted against the distance from the center of the front face of the cylinder. Though the distribution on the front (A–B) is hardly
814
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 7. Dependance of CD and CL on Sc.
Fig. 8. Pressure distribution on cylinder.
aKected by the Sc value, the cylinder oscillation at Sc = 0:13 yields the lowest pressure on the rear face (C–D). The time-averaged drag coe2cient CD depends on the pressure diKerence between the front and rear faces. Therefore, CD takes its maximum value at Sc = 0:13, as seen in Fig. 7.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
815
Fig. 9. Time-variation of volumetric velocity of liquid-phase l ul at Sc = 0:1 (asynchronous state).
3.3. Time-variation of 3ow properties The time-variations of the volumetric velocity of the liquid-phase l ul around the cylinder are presented in Figs. 9 and 10. The instantaneous distributions at four time points during one period of the cylinder oscillation are displayed. Figs. 9(a,c) and 10(a,c) are the results when the cylinder passes right and left, respectively, through the oscillation center with its maximum speed, Figs. 9(b,d) and 10(b,d) are the results when the cylinder is at rest. When the cylinder oscillation frequency is below the lock-in region, l ul distributes as shown in Fig. 9, where Sc = 0:1. Comparing Figs. 9(a) and (c), which are the distributions when the cylinder passes through the oscillation center, the #ow similarity is not observed. The vortices are shed independent of the cylinder oscillation. Fig. 10 shows the distribution at the lock-in state (Sc = 0:15). Comparing the #ow elds in Figs. 10(a) and (c), large vortices are produced just behind the cylinder, and the #ows are symmetrical with respect to the center line of the cylinder. This indicates the synchronization of the cylinder oscillation and the vortex shedding. In addition, the vortices behind the cylinder are closer to the cylinder than
816
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 10. Time-variation of volumetric velocity of liquid-phase l ul at Sc = 0:15 (synchronous state).
those in Fig. 9. This upstream migration of the vortex is also reported under single-phase #ow condition [16]. The distribution of the pressure coe2cient Cp changes with the lapse of time as shown in Figs. 11–13. Fig. 11 shows the distribution for the stationary cylinder, where the results at four time points during one period of the vortex shedding are presented. The position where Cp has the minimum value corresponds to the center of the vortex. The minimum value increases in the downstream direction, because the vortex intensity is diminished due to the #uid viscosity. The pressure is higher locally at the front of the cylinder. Figs. 12 and 13 show the distributions corresponding to the results in Figs. 9 and 10, respectively. The instantaneous distributions at four time points during one period of the cylinder oscillation are displayed. Figs. 12(a,c) and 13(a,c) are the results when the cylinder passes right and left, respectively, through the oscillation center with its maximum speed, Figs. 12(b,d) and 13(b,d) are the results when the cylinder is at rest. Fig. 12 shows the distribution when the cylinder oscillation frequency is below the lock-in region. The Sc value is 0:1. The distribution upstream of the cylinder is almost the same as that for the stationary cylinder (Fig. 11). But the Cp value heightens slightly in the wake, demonstrating the
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 11. Time-variation of pressure coe2cient Cp for stationary cylinder.
Fig. 12. Time-variation of pressure coe2cient Cp at Sc = 0:1 (asynchronous state).
817
818
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 13. Time-variation of pressure coe2cient Cp at Sc = 0:15 (synchronous state).
decrement in the vortex intensity. Therefore, the CL value at Sc = 0:1 is smaller than that for the stationary cylinder, as found in Fig. 7. The distribution at the lock-in state is shown in Fig. 13, where Sc = 0:15. When comparing with the results for the stationary cylinder (Fig. 11), the Cp value lessens in the wake and such lower Cp region spreads in the downstream direction. Because the vortex intensity is heightened by the synchronization of the vortex shedding with the cylinder oscillation. The time-variations of the volumetric fraction of the gas-phase g , corresponding to the results in Figs. 11–13, are displayed in Figs. 14 –16, respectively. Fig. 14 shows the distribution for the stationary cylinder. Since the gas-phase tends to swarm in the vortex due to the pressure gradient, g is higher in the wake [14,15]. Such region with higher g value spreads downstream of the cylinder. g is lower near the front of the cylinder because there is a positive pressure gradient toward the front. Fig. 15 shows the result at Sc = 0:1. The distribution upstream of the cylinder is almost the same as that for the stationary cylinder (Fig. 14). But, the g value in the wake is not so high. This is because the gas-phase is less able to swarm behind the cylinder owing to the gentler pressure gradient. The distribution at the lock-in state (Sc = 0:15) is shown in Fig. 16. In comparison with the result for the stationary cylinder (Fig. 14), the g value is higher in the wake and such higher region is wider downstream of the cylinder. Because the vortex intensity is higher, and accordingly the gas-phase shifts well toward such lower pressure region.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 14. Time-variation of volumetric fraction of gas-phase g for stationary cylinder.
Fig. 15. Time-variation of volumetric fraction of gas-phase g at Sc = 0:1 (asynchronous state).
819
820
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
Fig. 16. Time-variation of volumetric fraction of gas-phase g at Sc = 0:15 (synchronous state).
4. Summary The bubbly #ow around a square-section cylinder, which is forced to oscillate sinuously normal to the uniform #ow, is calculated by using an incompressible two-#uid model. The width of cylinder is 30 mm. A bubbly #ow is assumed, and the bubble diameter is set as 1 mm. The volumetric fraction of the gas-phase upstream of the cylinder is 0.05, where the Reynolds number de ned by the volumetric velocity of the liquid-phase is 200. The numerical results are summarized as follows: (1) With an increment in the cylinder oscillation frequency Sc, the vortex shedding frequency St approaches the Sc value. When Sc is increased beyond a critical value, the lock-in occurs, in which the vortex shedding is synchronized with the cylinder oscillation. The critical value is nearly the same as the St value for the stationary cylinder. (2) The time-averaged drag coe2cient and the root-mean-square of the lift coe2cient for the cylinder take their maximum values at the lock-in state. (3) The volumetric fraction of the gas-phase g is higher in the wake of the cylinder, irrespective of the Sc value. The g value is lower below the lock-in state and higher in the lock-in state than the value for the stationary cylinder. References [1] M.J. Pettigrew, D.J. Gorman, Experimental studies on #ow-induced vibration to support steam generator design, (1977) AECL-5804.
T. Uchiyama / Finite Elements in Analysis and Design 38 (2002) 803 – 821
821
[2] F. Axisa, B. Villard, R.J. Gibert, G. Hetsroni, P. Sundheimer, Vibration of tube bundles subjected to air–water and steam–water cross-#ow, Proceedings of the ASME Winter Annual Meeting, 1984, pp. 269 –283. [3] T. Nakamura, K. Fujita, A. Tsuge, Two-phase cross-#ow-induced vibration of tube arrays (Brief review of previous studies and summary of design methods), JSME Int. J. Ser. B 36 (3) (1993) 429–438. [4] F. Hara, Air-bubble eKects on vortex-induced vibrations of a circular cylinder, Proc. ASME Symp. FIV 1 (1984) 103–113. [5] F. Hara, T. Iijima, T. Nojima, Vibration response characteristics of tandem two circular cylinders subjected to a two phase cross-#ow — eKects of pitch-to-diameter ratio, Trans. ASME, J. Pressure Vessel Tech. 114 (4) (1992) 444–452. [6] Y. Joo, V.K. Dhir, An experimental study of drag on a single tube and on a tube in an array under two-phase cross-#ow, Int. J. Multiphase Flow 20 (6) (1994) 1009–1019. [7] O.M. Gri2n, M.S. Hall, Review — vortex shedding lock-on and #ow control in bluK body wakes, Trans. ASME, J. Fluid Eng. 113 (4) (1991) 526–537. [8] T. Uchiyama, K. Minemura, Numerical simulation of bubbly #ow around an oscillating hydrofoil, Proceedings of the Asian Symposium on Multiphase Flow, Osaka, 1999, pp. 27–32. [9] T. Uchiyama, ALE nite element method for gas–liquid two-phase #ow including moving boundary based on an incompressible two-#uid model, Nucl. Eng. Des. 205 (2001) 69–82. [10] T. Belytschko, D.P. Flanagan, J.M. Kennedy, Finite element methods with user-controlled meshes for #uid-structure interaction, Comput. Methods. Appl. Mech. Eng. 33 (1982) 669–688. [11] M. Ishii, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975. [12] T. Uchiyama, Numerical analysis of air–water two-phase #ow across a staggered tube bundle using an incompressible two-#uid model, Nucl. Sci. Eng. 134 (3) (2000) 281–292. [13] T. Uchiyama, Finite element analysis of gas–liquid two-phase #ow across an in-line tube bundle, Hybrid Method Eng. 2(4) (2001), accepted and in press. [14] T. Uchiyama, Numerical simulation of gas–liquid two-phase #ow around a rectangular cylinder by the incompressible two-#uid model, Nucl. Sci. Eng. 133 (1) (1999) 92–105. [15] T. Uchiyama, Petrov–Galerkin nite element method for gas–liquid two-phase #ow based on an incompressible two-#uid model, Nucl. Eng. Des. 193 (1999) 145–157. [16] P.W. Bearman, E.D. Obasaju, An experimental study of pressure #uctuatiions on xed and oscillating square-section cylinders, J. Fluid Mech. 119 (1982) 297–321.