Ocean Engineering 28 (2001) 1545–1555
A numerical study of three-dimensional wave interaction with a square cylinder Chi-Wai Li *, Pengzhi Lin Department of Civil and Structural Engineering, The Hong Kong Polytechnic University, Kowloon, Hung Hom, Hong Kong Received 31 July 2000; accepted 8 November 2000
Abstract In this study, a three-dimensional numerical model is used to study the wave interaction with a vertical rectangular pile. The model employs the large eddy simulation (LES) method to model the effect of small-scale turbulence. The velocity and vorticity fields around the pile are presented and discussed. The drag and inertial coefficients are calculated based on the numerical computation. The calculated coefficients are found to be in a reasonable range compared with the experimental data. Additional analyses are performed to assess the relative importance of drag and initial effects, which could be quantified by the force-related Keulegan and Carpenter (KC) number: KCf=UT/(4pL). Here U is the maximum fluid particle velocity, T the wave period and L the length of structure aligned with the wave propagation direction. For small KCf, the effective drag coefficient is proportional to 1/KCf, provided the wavelength is much longer than the structural length. When wavelength is comparable to the structure dimension, the effective drag coefficient would be reduced significantly due the cancellation of forces, which has been demonstrated by numerical results. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Wave–structure interaction; Wave forces; Large eddy simulation (LES); Square cylinder; Numerical model
1. Introduction Wave–structure interaction has been the focus of studies in coastal and offshore engineering for decades. Although some of the wave–structure interaction problems * Corresponding author. Fax: +852-2334-6389. E-mail address:
[email protected] (C.-W. Li). 0029-8018/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 9 - 8 0 1 8 ( 0 1 ) 0 0 0 1 1 - 7
1546
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
can be treated by two-dimensional models (e.g., Liu et al., 1999), a majority of practical problems require the use of three-dimensional models. For example, to study the stability of piles of an offshore or coastal structure, a three-dimensional model is needed so that the complicated flow and turbulence patterns as well as the force fields induced by different wave attacks can be accurately predicted. Recently, Li and Wang (2000) proposed a quasi-three-dimensional numerical model to study the turbulent channel flow. The large eddy simulation (LES) model was used to predict the turbulence structure. In this study, the above model is extended to the fully three-dimensional by employing the s-coordinate transformation (Lin and Li, 2000). The model has the capability of solving turbulent flow around structures with arbitrary shapes. In the following paragraph, a brief description of the numerical model will be given. The model is then used to study wave interaction with a vertical rectangular pile. The flow pattern around the structure is provided and discussed. The emphasis of the study will be on the prediction of drag and inertial coefficients based on the numerical results, which has an important application to the designs of coastal and offshore structures. To facilitate the engineering design, a new method is proposed to estimate the maximum total wave force on the structure based on the limited wave parameters that are readily obtained in the field.
2. Model description The governing equations for the spatially averaged mean flow are as follows, ∂Ui ⫽0 ∂xi
(1)
∂Ui 1 ∂p¯ ∂ ∂Ui ⫹ (t ⫹Rij ) ⫹Uj ⫽⫺ ∂t ∂xj r∂xi ∂xj ij
(2)
where variables with overbars are spatially averaged quantities and Ui is the mean velocity in the ith component, p the pressure, tij the stress by molecular viscous effect, and Rij the stress due to the averaging process. The term Rij is modeled by the subgrid scale (SGS) model as follows:
冉
冊
Rij ⫽2 L2S冑2Sij Sij Sij
(3)
where Sij is the strain rate of the mean flow and LS is the characteristics length scale which equals Cs(⌬1⌬2⌬3)1/3 with Cs=0.15. The numerical model solves Eqs. (1) and (2) using the operator splitting method (Li and Yu, 1996). The s-coordinate transformation is used to map the irregular physical domain into a cube where the free surface boundary condition could be precisely applied. The details of numerical implementation are referred to in Lin and Li (2000).
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
1547
3. Numerical simulation 3.1. Problem setup In this study, a vertical pile with the dimension of 1 m×1 m is deployed in the computational domain that measures 30 m×10 m in the horizontal plane. The center of the pile is 10 m away from the left boundary and in the middle of the y-direction. The still water depth is 1 m. A linear wave train with a wave period of 4 s and wave height of 0.05 m is sent from the left boundary. A non-uniform mesh system with a total number of 130×80 is used on the horizontal plane with the finest grid ⌬x=⌬y=0.05 m arranged near the pile (see Fig. 1). In total 20 uniform grids are used in the vertical direction. A constant ⌬t=0.004 s is used that ensures numerical stability during the entire simulation. The computation is carried out up to t=40 s when 10 waves have been sent from the boundary. The radiation boundary condition is used on the downstream side and the rest of the lateral boundaries are assumed to be solid walls. 3.2. Velocity and vorticity fields on the horizontal plane In Fig. 2, the velocity field around the pile during one wave period is plotted. The velocity vectors are plotted every three computational nodes. It is observed that vortices are generated at four corners of the pile and these vortices are confined to being near the structure. This is different from the case of the vertical pile exposing to the uniform flow, where the vortices generated at the corners are carried away from the structure and form regular vortex shedding with certain frequency (Li and Lin, 2000). The major reason that the vortices are prevented from leaving the structure in this study is due to the oscillatory motion of fluid particle movement under the wave.
Fig. 1. Mesh arrangement near the pile; lines are plotted every two grid nodes for easier visibility.
1548
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
Fig. 2.
Velocity field around the pile during one wave period on the middle elevation.
In this case, the value of KC=UT/L苲0.6, suggesting that the fluid particle can never travel across the structure during one wave period. This feature is illustrated more clearly in the vorticity plot (Fig. 3), which shows that the vorticity fields are surprisingly similar during one wave period, although the velocity has changed directions. Because the velocity and vorticity fields near the structure are quite different under wave and uniform flow such as current, the resulting force field also exhibits quite different features. We shall discuss this issue in more detail in the following context. The only exception is for a very long wave (e.g., tsunami) that has a large wave period. The corresponding KC in such a case is very large and thus the force acting on the structure can be analyzed in a similar way as for the uniform flow. In this study, however, we shall focus on the cases with low KC, which are more commonly encountered in coastal and offshore regions. Another thing worth noticing is that under real wave conditions, the velocity field changes not only in time but also in space. This is obvious from Fig. 2 (t/T=7.4) where velocity may change direction in front of and behind the structure. This feature is different from most of laboratory studies employing U-tubes to simulate wave conditions, where the flow oscillates only in time rather than in space. When the value of KC becomes so small that the wavelength is comparable to the structural length, the force field induced by the real wave could be significantly different from
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
1549
Fig. 3. Vorticity field around the pile during one wave period on the middle elevation; the color from black to white represents vorticity strength from ⫺0.5 l/s to 0.5 l/s.
that in the U-tube experiments. We will discuss this effect later in more detail when the wave forces are analyzed. 3.3. Free surface profile at the central cross section of the pile When the upper boundary is set to be solid, the mean flow is basically two-dimensional. This type of problem has been investigated intensively by both numerical and experimental means (e.g., Rodi, 1997; Lyn et al., 1995). When the free surface is present, the problem will become more complicated. The pressure difference on the front and rear sides of the pile due to flow separation, which is the direct reason for net wave force on the structure, is also reflected in the abrupt change of the free surface displacement nearby. This can be confirmed in Fig. 4 that provides the free surface profile at the centerline of the basin. It is observed that during the passage of the wave crest and trough, there are time delays of wave propagation due to the presence of the pile. This causes the larger free surface difference near the pile than otherwise no structure is present. In this particular case of a small amplitude wave,
1550
Fig. 4.
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
Free surface profile at the centerline of the computational basin within one wave period.
the vertical fluid motion is relatively weak. Thus, it is an acceptable approximation to use the free surface difference to estimate the wave force on the structure. 3.4. Total wave force on the pile A more precise way of calculating wave force on the structure, however, is to integrate the pressure around the structure and take the average over the pile height to obtain the average force per unit length. If we are concerned about the maximum force, we could also take the pressure difference right on the bottom surrounding the pile to estimate the total force per unit length. It is found that both methods render almost the same numerical results. In the following study, unless otherwise mentioned, the latter method is used to calculate the net forces. It is also noted that the shear force, which is normally two orders of magnitude smaller than the pressure force, is neglected when the total force is calculated. In Fig. 5, the total forces in the x- and y-directions, Fx and Fy, are given. The force in the x-direction is orders of magnitude greater than that in the y-direction.
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
Fig. 5.
1551
Calculate total forces in the x- and y-directions.
This is another distinct feature of vertical pile subject to a uniform current. In the latter case, both forces are in the same order of magnitude and the force in the xdirection has fairly constant values over time, while the force in the y-direction oscillates strongly due the periodic vortex shedding from the structural corners (Li and Lin, 2000). However, in this study, the horizontal force is induced by the oscillatory ambient flow and thus it has the same frequency as the incident wave. In addition, since there is no vortex shedding occurring, the flow pattern is essentially symmetric about the pile in the y-direction. As a result, the force in this direction is close to zero. 3.5. Drag and inertial coefficient Because the total force is dependent upon the scale of the problem, in the engineering application, the non-dimensional drag and inertial coefficients, CD and CM, are more commonly used to compare different designs. Based on Morison et al. (1950), the total wave force acting on the marine structure could be written as the sum of drag force and inertial force, Du 1 Fx⫽ CDrAu兩u兩⫹CMrV 2 Dt
(4)
where A is the projected area of structure on the plane perpendicular to the wave propagation and V is the volume of the structure. There are two methods that could be employed to estimate the values of CD and CM. One is the so-called zero-crossing method. According to the theory, for a linear
1552
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
wave train, when u reaches its maximum value under the wave crest, Du/Dt=0. Thus, based on the calculated Fx at this moment, it is possible to estimate CD. On the other hand, when Du/Dt reaches its maximum, u=0. Thus, one can estimate the value of CM. This method is only applicable when the data are quite regular with nearly no scattering. In this study, since the data are rather smooth (see Fig. 5 for Fx), this method can be used. The values of CD and CM are found to be 2.2662 and 1.0633, respectively. They are both in the reasonable range based on the experimental data. Another method of estimating CD and CM is to use the least square method that minimizes the total errors between the theoretical results and measured or calculated data (e.g., Dean and Dalrymple, 1991). This method can handle the scattering data but the accuracy for one coefficient (e.g., CD) may degenerate when one effect (e.g., inertial force) is dominant over another effect (e.g., drag force). In this study, based on the least square method, the values of CD and CM are calculated to be 2.9269 and 1.0454, respectively. The inertial coefficient is almost the same as that estimated by the zero-crossing method. The drag coefficient, however, is about 30% larger. This is because the drag effect in this problem is much smaller than the inertial effect and thus the errors of calculating CD are relatively large. Nevertheless, when substituting two sets of coefficients into Eq. (4), we are able to recover the total force field quite accurately with both methods (see Fig. 6). From Fig. 6, it is also observed that the inertial force is dominant over the drag
Fig. 6. Comparisons of calculated Fx (solid line) and predicted Fx based on Eq. (4) (circle); the drag and inertial forces are represented by a dashed line and a dotted line, respectively.
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
1553
force in this case with the former about one order of magnitude greater than the latter. Thus, even though the values of CD are quite different from the two methods, the recovered total forces are not affected significantly. In theory, there is a 90° phase shift between the inertial and drag forces and the total force has the phase in between. In the studied problem, since the inertial force is one order of magnitude greater, the total force has almost the same phase as the inertial force. Fig. 6 also implies that the values of CD and CM cannot directly reflect the relative importance of the two effects. In order to evaluate which force is dominant over the other, the other parameters are needed to characterize the problem. 3.6. Effective drag coefficient and force related KC number Although the separate use of drag and inertial coefficients represents better physics when studying wave forces on a structure, sometimes it is difficult to implement this idea in practice. This is because the acceleration of fluid particles is a difficult quantity to measure in the field, especially for random waves in the ocean. In general, the readily accessible parameters are the maximum fluid particle velocity and the peak frequency of the wave group. Thus, the more useful formula to predict the total wave force on the structure should look like Du 1 1 Fx⫽ CDrAu兩u兩⫹CMrV ⫽ CDErAu兩u兩 2 Dt 2
(5)
where CDE is the effective drag coefficient. In order to assess the relative importance of drag and inertial forces, we make the following approximations: V苲AL and Du/Dt苲2p/Tu, where L is the length of structure aligned with the wave propagation direction. It is to be noted that when we approximate the time derivative of velocity, we based it on the linear wave theory and neglect the phase information. This is only valid when one force is dominant over the other. In this study, we have been focusing on the case where the KC number is small and the inertial force is much larger than the drag force. Grouping the drag and inertial forces, we then have:
冉
冊
4pL CDE⫽CD 1⫹ 兩u兩T
(6)
If we further define the force related KC number as: KCf⫽
冉 冊 兩u兩T 4pL
(7)
Then the effective drag coefficient can be calculated to be: 1+KCf C CDE⫽ KCf D
(8)
1554
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
In most practical cases, CD is between 1 and 3. Substituting the value of CD into Eq. (8), we are able estimate the effective drag coefficient with the given KCf number that could be easily derived from the basic wave parameters and structural dimension. Further substituting CDE into Eq. (5), we are able to estimate the maximum wave force acting on the structure that could be used as the reference value for engineering designs. In order to verify the above proposal, we ran a few more numerical experiments with different wave periods and different wave amplitudes such that the KCf will change within a wide range. The effective drag coefficients for different cases are obtained by using Eq. (5) based on the calculated force and maximum particle velocity (under the wave crest). The values of CDE are plotted in Fig. 7, together with the curve based on Eq. (8) using CD=1.1. It is found that for the cases of relatively long waves (kL⬍1, where k is the wave number), the calculated points are in very good agreement with the proposed curve. However, as the wavelength becomes shorter (kL⬎1), the calculated CDE deviates from the curve dramatically. This is because as the wavelength reduces, the structure may be exposed to the quickly varying flow field and certain cancellation of positive and negative forces becomes possible. In particular, when the wavelength is equivalent to the structural length, almost complete force cancellation could result, which causes a very small value of CDE. Apparently, when such a short wave is studied, the additional parameter kL should also be considered in Eq. (8). This problem will be investigated in a future study. 4. Conclusion Wave interaction with a vertical rectangular pile is investigated numerically in this study. It is found that the velocity and vorticity fields near the structure exhibit
Fig. 7. Calculated (symbols) and predicted (line) effective drag coefficient CDE as a function of force related KCf; circle: kL=0.52; cross: KL=6.29, 2.82, 1.88, 1.20, and 0.52 from left to right.
C.-W. Li, P. Lin / Ocean Engineering 28 (2001) 1545–1555
1555
quite different features from the case of a pile exposed to a uniform current. For small KC numbers, the inertial force normally dominates over the drag force by orders of magnitude. To facilitate the engineering design, an effective drag coefficient CDE is proposed that is the function of the force related KC number. By using CDE, it is possible to estimate the maximum wave force, given the readily accessible wave parameters. It is found that as long as the wavelength is comparably larger than the structural dimension, the proposed way of calculating CDE is fairly accurate. Future work will include the effect of wavelength–structural dimension (characterized by kL) into the formulation.
Acknowledgements The research work herein was supported by a research grant from the Hong Kong Polytechnic University (G-YW28).
References Dean, R.G., Dalrymple, R.A., 1991. Water Wave Mechanics for Engineers and Scientists. In: Advanced Series on Ocean Engineering, vol. 2. World Scientific Publishing Co. Ltd. Li, C.W., Lin, P., 2000. Simulation of three-dimensional flow and turbulence structures around a rectangular pier in channel flow. 8th Int. Symp. Stochastic Hydraul., Peking, PR China pp. 173–184. Li, C.W., Wang, J.H., 2000. Large eddy simulation of free surface shallow water flow. Int. J. Num. Met. Fluids 34 (1), 31–46. Li, C.W., Yu, T.S., 1996. Numerical investigation of turbulent shallow recirculating flows by a quasithree-dimensional k⫺e model. Int. J. Num. Met. Fluids 23, 485–501. Lin, P., Li, C.W., 2000. A s-coordinate three-dimensional numerical model for surface wave propagation. Int. J. Num. Met. Fluids (submitted for publication). Liu, P.L.-F., Lin, P., Chang, K.-A., Sakakiyama, T., 1999. Numerical modeling of wave interaction with porous structures. J. Waterway, Port, Coastal and Ocean Engng 125 (6), 322–330. Lyn, D.A., Einav, S., Rodi, W., Park, J.-H., 1995. A laser-Doppler velocimetry study of ensemble-averaged characteristics of turbulent near wake of a square cylinder. J. Fluid Mech. 304, 285–319. Morison, J.R., O’Brien, M.P., Johnson, J.W., Schaaf, S.A., 1950. The force exerted by surface waves on piles. Petrol. Trans. AIME 189, 149–154. Rodi, W., 1997. Comparison of LES and RANS calculations of the flow around bluff bodies. J. Wind Engng. Indust. Aerodyn. 69-71, 55–75.