- Email: [email protected]
Numerical calculation of hydrodynamic characteristics of tidal currents for submarine excavation engineering in coastal area
Accepted Manuscript Numerical calculation of hydrodynamic characteristics of tidal current for submarine excavation engineering in coastal area Jian-hua Li, Liang-sheng Zhu, Shan-ju Zhang PII:
S1674-2370(16)30017-5
DOI:
10.1016/j.wse.2016.06.005
Reference:
WSE 60
To appear in:
Water Science and Engineering
Received Date: 30 December 2015 Accepted Date: 1 April 2016
Please cite this article as: Li, J.-h., Zhu, L.-s., Zhang, S.-j., Numerical calculation of hydrodynamic characteristics of tidal current for submarine excavation engineering in coastal area, Water Science and Engineering (2016), doi: 10.1016/j.wse.2016.06.005. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Numerical calculation of hydrodynamic characteristics of tidal current for submarine excavation engineering in coastal area Jian-hua LI, Liang-sheng ZHU*, Shan-ju ZHANG College of Civil and Transportation Engineering, South China University of Technology, Guangzhou 510640, China Received 30 December 2015; accepted 1 April 2016
M AN U
SC
RI PT
Abstract: In coastal areas with complicated flow movement, deposition and scour readily occur in submarine excavation projects. In this study, a small-scale model, with a high resolution in the vertical direction, was used to simulate the tidal current around a submarine excavation project. The finite volume method was used to solve Navier-Stokes equations and the Reynolds stress transport equation, and the entire process of the tidal current was simulated, with unstructured meshes generated in the irregular shape area, and structured meshes generated in other water areas. The meshes near the bottom and free surface were densified with a minimum layer thickness of 0.05 m. The volume of fluid method was used to track the free surface, the volume fraction of cells on the upstream boundary was obtained from the volume fraction of adjacent cells, and that on the downstream boundary was determined by the water level process. The numerical results agree with the observed data, and some conclusions can be drawn: after the foundation trench excavation, the flow velocity decreases quite a bit through the foundation trench, with reverse flow occurring on the lee slope in the foundation trench; the swirling flow impedes inflow, leading to the occurrence of dammed water above the foundation trench; the turbulent motion is stronger during ebbing than in other tidal stages, the range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the maximum value of Reynolds shear stress occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle. The numerical calculation method shows a strong performance in simulation of the hydrodynamic characteristics of tidal currents in the foundation trench, providing a basis for submarine engineering construction in coastal areas. Keywords: Small-scale model; Tidal current; Hydrodynamic characteristic; Coastal area; Submarine excavation engineering; Reynolds stress model
1. Introduction
AC C
EP
TE D
Submarine excavation projects should be constructed precisely when they involve complicated construction conditions, such as channel excavation engineering, submarine pipeline engineering, and submarine immersed tube tunnel engineering. Back silting and slope damage usually occur in excavation projects, owing to deposition and scour. This problem has crucial implications for engineering construction. For example, the 15th immersed tube of the Hong Kong-Zhuhai-Macau Bridge could not be placed in the foundation trench due to serious sedimentation (Li, 2015). When immersed tubes were sunk and installed, the lateral force, drifting force, and mooring force acting on the tubes changed significantly because of the hydrodynamic change (Lu et al., 2014). Therefore, investigation of the hydrodynamic characteristics for submarine excavation projects in coastal areas is highly significant. Tidal action is one of the most important impact factors in a coastal area. The horizontal scale of a tidal current is much greater than its vertical scale. Based on the static pressure hypothesis, three-dimensional shallow water equations have been derived using the water level variable instead of the pressure variable (Liu and He, 2000). According to previous studies, three-dimensional large-scale hydrodynamic models, such as the Princeton ocean model (POM) (Ohashi and Sheng, 2013; Gao et al., 2013), the environmental fluid dynamics code model (EFDCM) (Zhou et al., 2014; Ren et al., 2015), the estuarine, coastal, and ocean model (ECOM) (Ningsih and Azhar, 2013), the coupled hydrodynamical ecological model for regional and shelf seas (COHERENS) (Shi et al., 2010; Tuomi et al., 2012), and the finite-volume coastal ocean model (FVCOM) (Chen et al., 2007), have been commonly used to study the spatial and temporal variability of tidal current movement. Their commonly used turbulence submodels are the Smagorinsky, the Mellor and Yamada, and the k -ε models (Abbott, 1997). The sigma coordinate mode along the vertical direction has also been used to fit both the free surface and bed topography (Berntsen and Furnes, 2005). ————————————— This work was supported by the National Natural Science Foundation of China (Grant No. 41406005). *Corresponding author E-mail address: [email protected] (Liang-sheng ZHU) 1
2.1. Governing equations
TE D
2. Mathematic model descriptions
M AN U
SC
RI PT
The large-scale hydrodynamic models mentioned MANUSCRIPT above have been limited to two-dimensional or weakly ACCEPTED three-dimensional flows over a large calculation range, meaning that they are only meshed at a low resolution in the vertical direction, considering the effect of computational efficiency. If the terrain changes drastically, the sigma coordinate mode cannot achieve a high resolution in both shallow and deep waters (Gräwe et al., 2015; Berntsen et al., 2015). It is also difficult to densify the local water area in the vertical direction (Tao and Zhang, 2007). Therefore, a small-scale model should be used to study the tidal current in a study region around submarine excavation projects. Small-scale models have been used to solve the Navier-Stokes equations and analyze flow structures in complicated calculation domains (Tang et al., 2015). They have the advantages of smaller grids, more layers, and higher resolutions in the vertical direction (Tang el al., 2014). Using the small-scale model in the simulation, the characteristics of unidirectional flow around the submarine structures were focused, and the currents at the upstream boundary were considered the inflow to the calculation domains (Wu and Wang, 2009; Lubin et al., 2010; Nielsen et al., 2013). However, the results could not reflect the entire process of reciprocating motion due to tidal currents, even though the hydrodynamic characteristics during ebbing and flooding were also studied. Lefebvre et al. (2014) investigated the dynamic behavior of the bedform-influenced flow field and roughness lengths over a large bedform during a tidal cycle by establishment of a small-scale model with more than 50 layers in the vertical direction. However, their study essentially simulated unidirectional flow in ebbing and flooding, respectively, and the hydrodynamic characteristics for transition were not studied. In this study, we aimed at obtaining the hydrodynamic characteristics of tidal currents for submarine excavation projects. The finite volume method was used to solve the Navier-Stokes equations and the Reynolds stress transport equation. With the foundation trench excavation in the Hong Kong-Zhuhai-Macau Bridge Project in the Pearl River Estuary as an example, the processes of tidal movement after the foundation trench excavation were simulated, and the numerical results were compared with the observed data for the purposes of analyzing the reasonability of the numerical calculation method used in this study.
The basic governing equations for mass and momentum conservation (Pope, 2000; Zhao et al., 2013) are as follows: ∂ρ ∂ (1) + ( ρ ui ) = 0 ∂t ∂xi ∂ui − ρ ui′u ′j + Si (2) µ ∂x j where ui is the velocity component in the x i direction, and i = 1, 2, 3 ; p is the pressure; − ρ ui′u ′j is the Reynolds stress; t is time; ρ is the density of fluid; µ is the dynamic viscosity coefficient; and S i is the momentum source term. Considering anisotropics of turbulent viscosity, the Reynolds stress model (Versteeg and Malalasekera, 2007) was used in this study. The Reynolds stress transport equation and turbulent dissipation rate transport equation (Pope, 2000; Duan et al., 2015) are shown below: ∂ ∂ ρ u i′u ′j + ρu r u i′u ′j = DT ,ij + D L,ij + Pij + G ij + φ ij + ε ij + Fij + S RS (3) ∂t ∂x r
(
AC C
EP
∂ ∂ ∂p ∂ ( ρ ui ) + ( ρ ui u j ) = − + ∂t ∂x j ∂xi ∂x j
)
(
∂ ∂ ∂ ( ρε ) + ( ρε ui ) = ∂t ∂xi ∂xi
)
µt µ + σε
∂ε 1 ε ε2 ρ + C P + C G − C + Sε ( ) ε1 ii ε 3 ii ε2 k k ∂xi 2
(4)
where DT,ij is the turbulent diffusion term; DL,ij is the molecular diffusion term; Pij is the stress production term; Gij is the buoyancy production term; φij is the pressure strain term; ε ij is the dissipation term; Fij is the production term according to the system rotation; k is the turbulent kinetic energy; ε is the turbulent dissipation rate; µ t is the turbulent viscosity; S RS and S ε are the source terms; σ ε is the Prandtl number 2
for ε with a value of 1.3; and C ε 1ACCEPTED , C ε 2 , and C ε 3MANUSCRIPT are constants with the values of 1.44, 1.92, and 0, respectively. DL,ij , Pij , and Fij only contain two-order correlation terms, and DT,ij . Gij , φij , and ε ij contain uncertain correlation terms, which can be calculated as follows: ∂ µ t ∂ui′u ′j DT,ij = (5) ∂xr σ k ∂xr Gij = 0 (6) 32 ε 2 1 ε 3 3 Ck − φij = −C1 ρ ui′u ′j − kδ ij − C2 Pij − Prr δ ij + C1′ρ ur′um′ nr nmδ ij − ui′ur′ n j nr − u ′j ur′ ni nr k 3 3 k 2 2 εd
M AN U
SC
RI PT
(7) Ck 3 2 1 3 1 3 1 C2′ C2 Prm − Prr δ rm nr nmδ ij − C2 Pir − Prr δ ir n j nr − C2 Pjr − Prr δ jr ni nr 3 2 3 2 3 εd 2 ε ij = ρεδ ij (8) 3 where σ k is the Prandtl number for k with a value of 1.0, δ ij is the Kronecker function, n r is the unit normal vector, and d is the distance to the wall. The values of C1 , C2 , C1′ , and C2′ are 1.8, 0.6, 0.5, and 0.3, respectively; and C = Cµ3 4 κ , where κ is the von Kármán constant with a value of 0.4187, and C µ = 0.09 . The air-water interface was considered in this study, and the volume of fluid (VOF) method was used to track the free surface in the model (Rui et al., 2015). In a cell, Fw and Fa are the volume fractions of the water phase and air phase, respectively, and 0 ≤ Fw (or Fa ) ≤ 1 . They are in the relation as follows: Fw + Fa = 1 (9)
2.2. Calculation of hydrodynamic characteristics
(
)
ε
TE D
Drastic turbulence takes place in the vicinity of submarine projects. The hydrodynamic characteristics can be analyzed by calculating the characteristic parameters (Launder and Spalding, 1974), including the turbulent kinetic energy k , the turbulent viscosity µ t , the Reynolds shear stresses τ xy , τ xz , and τ yz , and the wall shear stress τ w , as follows: 1 k = u1′u1′ + u 2′ u 2′ + u 3′ u 3′ (10) 2 k2 µ t = ρ Cµ (11)
(
12 ρκ uWA Cµ1 4 kWA
12 ln E ρ Cµ1 4 k WA d WA µ
)
(12)
(13)
AC C
τw =
EP
τ xy = τ yx = − ρ u1′u 2′ τ xz = τ zx = − ρ u1′u 3′ τ yz = τ zy = − ρ u 2′ u 3′
where uWA , k WA , and d WA are the mean velocity, the turbulent kinetic energy, and the distance to the wall of the wall adjacent cell, and E is an empirical constant with a value of 9.793.
2.3. Numerical calculation scheme 2.3.1. Numerical solution method
The finite volume method was used to solve Reynolds-averaged Navier-Stokes equations, combined with the Reynolds stress model (Pope, 2000). Considering the characteristics of tidal currents and computational efficiency, the pressure-implicit with splitting of operators (PISO) algorithm was used to calculate the pressure and velocity (Jang et al., 1986). The PISO algorithm is based on the higher degree of the approximate relation between the corrections for pressure and velocity. It can maintain a stable calculation with a larger time step for both momentum and pressure (Issa, 1986). The discretization schemes for the diffusion term and the convection term were carried out using a central difference scheme and a second-order wind scheme, respectively. The discretization schemes for the turbulent 3
kinetic energy, turbulent dissipation rate, and Reynolds stress were implemented using a first-order upwind ACCEPTED MANUSCRIPT scheme. The numerical oscillation was limited because the upwind scheme could maintain the stability of numerical calculation. The stability condition (Courant et al., 1967) in the calculation process is as follows:
CFL =
∑λ
i max
i
2V
Ai
≤1
(14)
where CFL is the Courant number, V is the cell volume, Ai is the area of cell face i, and λi max is the maximum of local eigenvalues.
RI PT
2.3.2. Mesh generation
AC C
EP
TE D
M AN U
SC
With the foundation trench excavation in the Hong Kong-Zhuhai-Macau Bridge Project in the Pearl River Estuary as an example, hybrid meshes were generated in the calculation domain. The mesh in the vertical direction was densified at a high resolution in a local region. Unstructured meshes were generated to fit the irregular shape of the structural sections and drastic terrain. Structured meshes were generated in other water areas and at the wall boundary. The meshes near the bottom and free surface were densified. Fig. 1 and Fig. 2(a) present the calculation domains and meshes before and after the foundation trench excavation. The wedge mesh was used for the area of the foundation trench, and the hexahedral mesh was used for the other parts. Along the vertical direction, the meshes near the bottom and free surface were densified with a minimum layer thickness of 0.05 m. According to the observed hydrologic data, there is a mainly reciprocating current in the area between the outlets of the Lingding Estuary and Lantau Island. Therefore, the x direction of the calculation domain was set along the north-to-south direction. The width of the foundation trench with its baseline perpendicular to the flow direction was 50 m. A composite slope was set on both sides with lower segment and upper segment gradients of 1:3 and 1:7, respectively (Cao and Xiao, 2012). Its bottom elevation was −45 m, and the bottom width along the flow direction was 42 m. The size and meshes of the foundation trench after excavation can be seen in Fig. 2(b).
Fig. 1. Schematic diagram of calculation domain and mesh generation before foundation trench excavation
4
RI PT
ACCEPTED MANUSCRIPT
Fig. 2. Schematic diagrams of calculation domain and mesh generation after foundation trench excavation
2.3.3. Boundary conditions
M AN U
SC
The varying processes of velocity and water level were set as the upstream and downstream boundary conditions, respectively. The x direction of the domain was set to be along the dominant flow direction, so the lateral boundaries along the x direction were set as symmetry boundaries. The other boundaries were set as wall boundaries. The time series of the velocity and water level data were obtained from the simulated results of the tidal currents in the Pearl River Estuary using a large-scale model. On the water level boundary, the pressure is in accordance with the assumption of static pressure distribution. It can be calculated with Ps = Pa + ρ g (η − z ) (15) where Ps is the absolute pressure on the water level boundary, Pa is the standard atmospheric pressure, g is the gravitational acceleration, η is the position of the free surface, and z is the vertical coordinate value. The volume fractions of cells on the upstream boundary were obtained from adjacent cells at the previous time step. For the water phase, the volume fraction for surface cell i on the upstream boundary is as follows: FwU ( i, t ′ ) = Fwat ′−1 (16)
EP
F ( i, t ′ ) = 0 wD FwD (i, t ′) = 1 F (i, t ′) = 1 − ( z0 − hwl ) wD 2 ∆z
TE D
where t ′ is the time step, and Fwat ′−1 is the volume fraction of the water phase of the adjacent cell at the previous time step. FwU (i,0) can be determined by the initial flow field. The volume fractions of cells on the downstream boundary were obtained based on the water level process. For the water phase, if ( x0 , y0 , z0 ) is the centroid of the free surface on the downstream boundary, the volume fraction for surface cell i on the downstream boundary is determined as follows:
zi > z0
(17)
AC C
zi < z0
zi = z0
where z i is the vertical coordinate of the centroid of cell i , hwl is the water level at the current time step, and ∆z is the grid spacing in the z direction.
3. Application of numerical calculation method 3.1. Validation In this study, the precision of the numerical calculation method was validated with the observed data at an observation station in the Lingdingyang Estuary under natural conditions. The observed data were obtained by one-day observation (from 12:00 on August 6 to 12:00 on August 7, 2008) using an acoustic wave and current meter. The corresponding monitoring point was set at the center of the calculation domain, located at (505, −22.5).
5
ACCEPTED MANUSCRIPT
3.1.1. Validation of water level
M AN U
SC
RI PT
The water level process calculated at the monitoring point matches the observed values, as shown in Fig. 3(a). The relative error is less than 0.4%. For the numerical results, the highest water level is 11.90 m and the lowest is 10.20 m, which are the same as the observed values. During flooding and ebbing, the calculated values are slightly higher than the observed values, indicating that the water level process obtained from the numerical simulation is reasonable.
Fig. 3. Comparisons of calculated and observed values of water level and depth-averaged velocity at monitoring point
3.1.2. Validation of depth-averaged velocity
EP
TE D
The time series of the calculated results of the depth-averaged velocity at the monitoring point match the observed values, as shown in Fig. 3(b). The calculated and observed results have the same period, and their maximum and minimum values appear at the same time. The average values of the calculated and observed results of the depth-averaged velocity during flooding (from 20:00 on August 6 to 01:00 on August 7) are 0.35 m/s and 0.39 m/s, respectively, with a relative error of 10.26%. The average values of calculated and observed results of the depth-average velocity during ebbing (from 01:00 to 07:00 on August 7) are 0.44 m/s and 0.51 m/s, respectively, with a relative error of 13.73%. At maximum flooding (at 00:00 on August 7), the calculated and observed values are 0.51 m/s and 0.55 m/s, respectively, with a relative error of 7.27%. At maximum ebbing (at 06:00 on August 7), the calculated and observed values are 0.63 m/s and 0.71 m/s, respectively, with a
AC C
relative error of 11.27%. Fig. 3(b) shows that the calculated results are slightly lower than the observed values. The error is caused by the discharge boundary data obtained from the calculated results of the Pearl River Estuary with the large-scale tidal current model. 3.1.3. Validation of vertical distribution of flow velocity Fig. 4 shows the vertical distributions of flow velocity (u) at the monitoring point at different times during flooding and ebbing. Figs. 4(a) through (d) show the velocity distributions during flooding, and Figs. 4(e) through (h) show the velocity distributions during ebbing. In the simulation, the changes in temperature and salinity were not considered. The distributions follow the logarithmic law (Ni et al., 2009), and the calculated results show errors as compared with the observed values, which are caused by the stratification and mixing phenomena. At 22:00 and 23:00 on August 6 and 02:00 and 03:00 on August 7, the distributions match the observed data, and the stratification is weak. However, when the stratification strengthens, the deviations between the calculated results and observed data grow. Moreover, the deviations during ebbing are relatively small.
6
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 4. Calculated results and observed values of vertical distributions of flow velocity at monitoring point
3.2. Engineering application
3.2.1. Hydrodynamic characteristics under natural conditions
EP
TE D
During a tidal cycle, the duration of ebbing is longer than that of flooding, and the water level rises slowly during flooding and drops rapidly during ebbing. Fig. 3 shows that the time of the maximum value of water level lags behind that of the maximum value of flow velocity during a tidal cycle. During flooding, the calculated maximum value of flow velocity is 0.51 m/s, occurring before the high tide level. In the initial stage of flooding, the flow velocity continues to decrease and reaches 0 after 21:00 on August 6. This is due to land runoff and the inertial force of ebbing. Then, it reverses and begins to increase. During ebbing, the calculated
AC C
maximum value of flow velocity is 0.63 m/s, occurring before the low tide level. In the initial stage of ebbing, the flow velocity continues to decrease and reaches 0 after 01:00 on August 7. Then, it reverses and begins to increase. The transition of flow velocity of tidal currents in the lower layer precedes that in the upper layer (Fig. 5), due to the inertia of water particle movement (Xie et al., 2011). In general, the flow velocity in the upper layer is larger than that in the lower layer. It is easier for water particles to overcome the inertia in the lower layer than in the upper layer. Therefore, it is easier for flow to reverse its direction in the lower layer.
7
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 5. Vertical distributions of flow velocity at monitoring point during transition
3.2.2. Hydrodynamic characteristics after foundation trench excavation
AC C
EP
TE D
After the foundation trench excavation, the area of flow section increases, and the flow slows down through the foundation trench. Fig. 6 presents the velocity distributions in the section of y = −22.5 m. The vertical distributions of flow velocity show wedge-shaped distributions above the lee slope. The velocities at the monitoring point decrease by 60% and 52% at maximum flooding and maximum ebbing after excavation, respectively, in contrast to those before excavation.
Fig. 6. Velocity distributions in section of y = −22.5 m at maximum flooding and maximum ebbing
Fig. 7 presents the flow fields in the section of y = −22.5 m in the foundation trench, with Figs. 7(a) through (f) showing the flow fields during flooding, and Figs. 7(g) through (l) showing the flow fields during ebbing. It indicates that at 22:00 on August 6 and 02:00 on August 7, the flow slows down without swirling flow
8
TE D
M AN U
SC
RI PT
occurring during the transition. Then,ACCEPTED swirling flow occurs due to the increasing velocity after the transition. MANUSCRIPT Furthermore, reverse flow mainly occurs on the lee slope. When the reverse flow moves upward along the lee slope, dammed water occurs above the foundation trench. The maximum change of water level is 0.036 m according to the water level difference before and after excavation. It can be concluded from Fig. 8 that the excavation has little influence on the water level. The results of flow fields match the results of Cao and Xiao (2012), which were obtained based on a two-dimensional horizontal numerical model and field tests.
AC C
EP
Fig. 7. Flow fields in section of y = −22.5 m in foundation trench
Fig. 8. Variations in water level before and after excavation
Intense turbulence takes place in the foundation trench due to the occurrence of swirling flow. According to Eqs. (10) through (13), the hydrodynamic characteristics were obtained. The maximum turbulent viscosities at maximum flooding and maximum ebbing are 100 kg/(m·s) after foundation trench excavation, while the corresponding values under natural conditions are 13.6 kg/(m·s) and 55.8 kg/(m·s), respectively, indicating that the excavation of the foundation trench intensifies the turbulence in this area. The comparison of turbulent characteristics in different tidal stages shows that the turbulent flow is stronger during ebbing than in other tidal stages. As shown in Fig. 9, the range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the strong turbulent motion mainly takes place in the water above the stoss slope. However, the turbulent viscosity decreases during the transition, as shown in Figs. 9(c) and (d). 9
RI PT
ACCEPTED MANUSCRIPT
Fig. 9. Distributions of turbulent viscosity µ t at maximum flooding, maximum ebbing, and trannsition in section of y = −22.5 m
EP
TE D
M AN U
SC
Fig. 10 shows that the values of Reynolds shear stress τ xz at maximum flooding and maximum ebbing are larger than those during the transition, and the maximum value of τ xz occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle, demonstrating that there exists intensive exchange of water mass in different layers, and the shear velocity shows a large degree of variation due to the occurrence of eddies.
Fig. 10. Distributions of Reynolds shear stress τ xz at maximum flooding, maximum ebbing, and transition in section of y = −22.5 m
AC C
Fig. 11 presents the wall shear stress distribution along the flow direction at maximum flooding, maximum ebbing, and the transition. The values of wall shear stress at maximum ebbing are largest due to the larger values of flow velocity. During the transition, the values of wall shear stress are small. Furthermore, the wall shear stress is higher on the north side, implying that the frictional resistance is larger on the north side, and the sediment incipient motion and transport occur more readily on the south side. The results match the phenomenon of the thickness of deposition on the north side being greater than that on the south side. Also, the hydrodynamic characteristics obtained in this study agree with a previous study by Xin et al. (2012), in which the hydrodynamic characteristics were obtained based on trial tunnel experiments over a long-period observation.
10
RI PT
ACCEPTED MANUSCRIPT
Fig. 11. Wall shear stress distribution along flow direction at maximum flooding, maximum ebbing, and transition
4. Conclusions
References
TE D
M AN U
SC
Based on the Navier-Stokes equations and the Reynolds stress model, a three-dimensional small-scale numerical model of tidal currents for submarine excavation projects was established. The results match the observed data and show a strong performance of the numerical model in simulation of the hydrodynamic characteristics of tidal currents. The following conclusions are drawn: (1) After the foundation trench excavation, the flow velocity decreases quite a bit through the foundation trench, and reverse flow occurs on the lee slope in the foundation trench. With the action of swirling flow, the reverse flow moves upward along the lee slope, and the inflow is impeded, leading to the occurrence of dammed water above the foundation trench. However, the excavation has little influence on the water level. (2) The turbulent motion is stronger during ebbing than in other tidal stages. The range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the maximum value of Reynolds shear stress occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle. The hydrodynamic characteristics obtained in this study demonstrate the phenomenon of the thickness of deposition on the north side being greater than that on the south side. The flow field characteristics and hydrodynamic characteristics of tidal currents around the submarine excavation project are reasonably described by the numerical calculation method.
AC C
EP
Abbott, M.B., 1997. Range of tidal flow modeling. Journal of Hydraulic Engineering 123(4), 257-277. http://dx.doi.org/10.1061/(ASCE)0733-9429(1997)123:4(257). Berntsen, J., Furnes, G., 2005. Internal pressure errors in sigma-coordinate ocean models: Sensitivity of the growth of the flow to the time stepping method and possible non-hydrostatic effects. Continental Shelf Research 25(7-8), 829-848. http://dx.doi.org/ 10.1016/j.csr.2004.09.025. Berntsen, J., Thiem, Ø., Avlesen, H., 2015. Internal pressure gradient errors in σ -coordinate ocean models in high resolution fjord studies. Ocean Modelling 92, 42-55. http://dx.doi.org/10.1016/j.ocemod.2015.05.009. Cao, H.J., Xiao, L.B., 2012. Back silting of trench excavation for the immersed tunnel in Hong Kong-Zhuhai-Macao Bridge construction. Port and Waterway Engineering (1), 12-17 (in Chinese). Chen, C.S., Huang, H.S., Beardsley, R.C., Liu, H.D., Xu, Q.C., Cowles, G., 2007. A finite volume numerical approach for coastal ocean circulation studies: Comparisons with finite difference models. Journal of Geophysical Research 112(C03018), 1-34. http://dx.doi.org/10.1029/2006JC003485. Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11(2), 215-234. http://dx.doi.org/10.1147/rd.112.0215. Duan, D.H., Zhang, J.B., Gao, Z.X., 2015. Investigation on spanwise wall oscillation in turbulent channel flow based on Reynolds stress model (RSM). Procedia Engineering 99, 665-670. http://dx.doi.org/10.1016/j.proeng.2014.12.587 Gao, J.S., Xue, H.J., Chai, F., and Shi, M.C. 2013. Modeling the circulation in the gulf of Tonkin, South China Sea. Ocean Dynamics 63(8), 979-993. http://dx.doi.org/10.1007/s10236-013-0636-y. Gräwe, U., Holtermann, P., Klingbeil, K., Burchard, H., 2015. Advantages of vertically adaptive coordinates in numerical models of stratified shelf seas. Ocean Modelling 92, 56-68. http://dx.doi.org/10.1016/j.ocemod.2015.05.008. Issa, R.I., 1986. Solution of implicitly discretized fluid flow equations by operator splitting. Journal of Computational Physics 62, 11
AC C
EP
TE D
M AN U
SC
RI PT
40-65. http://dx.doi.org/10.1016/0021-9991(86)90099-9. ACCEPTED MANUSCRIPT Jang, D.S., Jetli, R., Acharya, S. 1986. Comparison of the PISO, SIMPLER, and SIMPLEC algorithms for the treatment of the pressure-velocity coupling in steady flow problems. Numerical Heat Transfer 10(3), 209-228. http://dx.doi.org/10.1080/10407798608552506. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3(2), 269-289. http://dx.doi.org/10.1016/0045-7825(74)90029-2. Lefebvre, A., Paarlberg, A.J., Ernstsen, V.B., Winter, C., 2014. Flow separation and roughness lengths over large bedforms in a tidal environment: A numerical investigation. Continental Shelf Research 91, 57-69. http://dx.doi.org/10.1016/j.csr.2014.09.001. Li, Z.L., 2015. Three times installation of 15th immersed tube in Hong Kong-Zhuhai-Macau submarine tunnel engineering in Lingdingyang Bay. Construction and Architecture (19), 13-15 (in Chinese). Liu, H., He, Y.S., 2000. Recent development in 3-D hydrodynamic mathematical model studies for estuaries. The Ocean Engineering 18(2), 87-93 (in Chinese). Lu, W.Q., Wu, W.G., Su, L.W., Chen, K.Q., Peng, S., 2014. Hydrodynamic experimental study on immersed tube tunnel element of Hong Kong-Zhuhai-Macao Bridge. China Civil Engineering Journal 47(3), 138-144 (in Chinese). Lubin, P., Glockner, S., Chanson, H., 2010. Numerical simulation of a weak breaking tidal bore. Mechanics Research Communications 37(1), 119-121. http://dx.doi.org/10.1016/j.mechrescom.2009.09.008. Ni, Z.H., Song, Z.Y., Wu, L.C., 2009. Study on the double-logarithmic profile of tidal flow velocity in the near-bed layers. Acta Oceanologica Sinica 28(6), 84-92. Nielsen, A.W., Liu, X.F., Sumer, B.M., Fredsøe, J., 2013. Flow and bed shear stresses in scour protections around a pile in a current. Coastal Engineering 72, 20-38. http://dx.doi.org/10.1016/j.coastaleng.2012.09.001. Ningsih, N.S., Azhar, M.A., 2013. Modelling of hydrodynamic circulation in Benoa Bay, Bali. Journal of Marine Science and Technology 18(2), 203-212. http://dx.doi.org/10.1007/s00773-012-0195-9. Ohashi, K., Sheng, J.Y., 2013. Influence of St. Lawrence River discharge on the circulation and hydrography in Canadian Atlantic waters. Continental Shelf Research 58, 32-49. http://dx.doi.org/10.1016/j.csr.2013.03.005. Pope, S.B., 2000. Turbulent Flow. Cambridge University Press, Cambridge, pp. 83-95. Ren, L., Nash, S., Hartnett, M., 2015. Observation and modeling of tide- and wind-induced surface currents in Galway Bay. Water Science and Engineering 8(4), 345-352. http://dx.doi.org/10.1016/j.wse.2015.12.001. Rui, B.F., Falcão, D.S., Oliveira, V.B., Pinto, A.M.F.R., 2015. Numerical simulations of two-phase flow in proton exchange membrane fuel cells using the volume of fluid method: A review. Journal of Power Sources 277, 329-342. http://dx.doi.org/10.1016/j.jpowsour.2014.11.124. Shi, J.Z., Li, C., Dou, X.P., 2010. Three-dimensional modeling of tidal circulation within the north and south passages of the partially-mixed Changjiang River Estuary, China. Journal of Hydrodynamics, Ser. B 22(5), 656-661. http://dx.doi.org/10.1016/S1001-6058(10)60010-8. Tang, G.Q., Chen, C.Q., Zhao, M., Lu, L., 2015. Numerical simulation of flow past twin near-wall circular cylinders in tandem arrangement at low Reynolds number. Water Science and Engineering 8(4), 315-325. http://dx.doi.org/10.1016/j.wse.2015.06.002. Tang, H.S., Qu, K., Wu, X.G., 2014. An overset grid method for integration of fully 3D fluid dynamics and geophysics fluid dynamics models to simulate multiphysics coastal ocean flows. Journal of Computational Physics 273, 548-571. http://dx.doi.org/10.1016/j.jcp.2014.05.010. Tao, J.F., Zhang, C.K., 2007. Review on several vertical coordinate modes of 3-D flow numerical model in estuarine and coastal waters. The Ocean Engineering 25(1), 133-142 (in Chinese). Tuomi, L., Kai, M., Lehmann, A., 2012. The performance of the parameterisations of vertical turbulence in the 3D modelling of hydrodynamics in the Baltic Sea. Continental Shelf Research 50-51, 64-79. http://dx.doi.org/10.1016/j.csr.2012.08.007. Versteeg, H.K., Malalasekera, W., 2007. An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Second ed. Pearson Education Limited, Harlow. Wu, Z.C., Wang, D.Z., 2009. Numerical solution for tidal flow in opening channel using combined MacCormack-finite analysis method. Journal of Hydrodynamics, Ser. B 21(4), 505-511. http://dx.doi.org/10.1016/S1001-6058(08)60177-8. Xie, D.F., Pan, C.H., Wu, X.G., 2011. Three-dimensional numerical modeling of tidal bore in Qiantang based on FVCOM. The Ocean Engineering 29(1), 47-52 (in Chinese). Xin, W.J., Jia, Y.S., He, J., 2012. Back silting analysis of the trial dredged-trough for immersed tube tunnel of Hong Kong-Zhuhai-Macao Bridge. Hydro-Science and Engineering (2), 71-78 (in Chinese). Zhao, Y.P., Bi, C.W., Dong, G.H., Gui, F.K., Cui, Y., Guan, C.T., Xu, T.J., 2013. Numerical simulation of the flow around fishing plane nets using the porous media model. Ocean Engineering 62, 25-37. http://dx.doi.org/10.1016/j.oceaneng.2013.01.009. Zhou, J.T., Falconer, R.A., Lin, B.L., 2014. Refinements to the EFDC model for predicting the hydro-environmental impacts of a barrage across the Severn Estuary. Renewable Energy 62, 490-505. http://dx.doi.org/10.1016/j.renene.2013.08.012.
12
S1674-2370(16)30017-5
DOI:
10.1016/j.wse.2016.06.005
Reference:
WSE 60
To appear in:
Water Science and Engineering
Received Date: 30 December 2015 Accepted Date: 1 April 2016
Please cite this article as: Li, J.-h., Zhu, L.-s., Zhang, S.-j., Numerical calculation of hydrodynamic characteristics of tidal current for submarine excavation engineering in coastal area, Water Science and Engineering (2016), doi: 10.1016/j.wse.2016.06.005. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Numerical calculation of hydrodynamic characteristics of tidal current for submarine excavation engineering in coastal area Jian-hua LI, Liang-sheng ZHU*, Shan-ju ZHANG College of Civil and Transportation Engineering, South China University of Technology, Guangzhou 510640, China Received 30 December 2015; accepted 1 April 2016
M AN U
SC
RI PT
Abstract: In coastal areas with complicated flow movement, deposition and scour readily occur in submarine excavation projects. In this study, a small-scale model, with a high resolution in the vertical direction, was used to simulate the tidal current around a submarine excavation project. The finite volume method was used to solve Navier-Stokes equations and the Reynolds stress transport equation, and the entire process of the tidal current was simulated, with unstructured meshes generated in the irregular shape area, and structured meshes generated in other water areas. The meshes near the bottom and free surface were densified with a minimum layer thickness of 0.05 m. The volume of fluid method was used to track the free surface, the volume fraction of cells on the upstream boundary was obtained from the volume fraction of adjacent cells, and that on the downstream boundary was determined by the water level process. The numerical results agree with the observed data, and some conclusions can be drawn: after the foundation trench excavation, the flow velocity decreases quite a bit through the foundation trench, with reverse flow occurring on the lee slope in the foundation trench; the swirling flow impedes inflow, leading to the occurrence of dammed water above the foundation trench; the turbulent motion is stronger during ebbing than in other tidal stages, the range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the maximum value of Reynolds shear stress occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle. The numerical calculation method shows a strong performance in simulation of the hydrodynamic characteristics of tidal currents in the foundation trench, providing a basis for submarine engineering construction in coastal areas. Keywords: Small-scale model; Tidal current; Hydrodynamic characteristic; Coastal area; Submarine excavation engineering; Reynolds stress model
1. Introduction
AC C
EP
TE D
Submarine excavation projects should be constructed precisely when they involve complicated construction conditions, such as channel excavation engineering, submarine pipeline engineering, and submarine immersed tube tunnel engineering. Back silting and slope damage usually occur in excavation projects, owing to deposition and scour. This problem has crucial implications for engineering construction. For example, the 15th immersed tube of the Hong Kong-Zhuhai-Macau Bridge could not be placed in the foundation trench due to serious sedimentation (Li, 2015). When immersed tubes were sunk and installed, the lateral force, drifting force, and mooring force acting on the tubes changed significantly because of the hydrodynamic change (Lu et al., 2014). Therefore, investigation of the hydrodynamic characteristics for submarine excavation projects in coastal areas is highly significant. Tidal action is one of the most important impact factors in a coastal area. The horizontal scale of a tidal current is much greater than its vertical scale. Based on the static pressure hypothesis, three-dimensional shallow water equations have been derived using the water level variable instead of the pressure variable (Liu and He, 2000). According to previous studies, three-dimensional large-scale hydrodynamic models, such as the Princeton ocean model (POM) (Ohashi and Sheng, 2013; Gao et al., 2013), the environmental fluid dynamics code model (EFDCM) (Zhou et al., 2014; Ren et al., 2015), the estuarine, coastal, and ocean model (ECOM) (Ningsih and Azhar, 2013), the coupled hydrodynamical ecological model for regional and shelf seas (COHERENS) (Shi et al., 2010; Tuomi et al., 2012), and the finite-volume coastal ocean model (FVCOM) (Chen et al., 2007), have been commonly used to study the spatial and temporal variability of tidal current movement. Their commonly used turbulence submodels are the Smagorinsky, the Mellor and Yamada, and the k -ε models (Abbott, 1997). The sigma coordinate mode along the vertical direction has also been used to fit both the free surface and bed topography (Berntsen and Furnes, 2005). ————————————— This work was supported by the National Natural Science Foundation of China (Grant No. 41406005). *Corresponding author E-mail address: [email protected] (Liang-sheng ZHU) 1
2.1. Governing equations
TE D
2. Mathematic model descriptions
M AN U
SC
RI PT
The large-scale hydrodynamic models mentioned MANUSCRIPT above have been limited to two-dimensional or weakly ACCEPTED three-dimensional flows over a large calculation range, meaning that they are only meshed at a low resolution in the vertical direction, considering the effect of computational efficiency. If the terrain changes drastically, the sigma coordinate mode cannot achieve a high resolution in both shallow and deep waters (Gräwe et al., 2015; Berntsen et al., 2015). It is also difficult to densify the local water area in the vertical direction (Tao and Zhang, 2007). Therefore, a small-scale model should be used to study the tidal current in a study region around submarine excavation projects. Small-scale models have been used to solve the Navier-Stokes equations and analyze flow structures in complicated calculation domains (Tang et al., 2015). They have the advantages of smaller grids, more layers, and higher resolutions in the vertical direction (Tang el al., 2014). Using the small-scale model in the simulation, the characteristics of unidirectional flow around the submarine structures were focused, and the currents at the upstream boundary were considered the inflow to the calculation domains (Wu and Wang, 2009; Lubin et al., 2010; Nielsen et al., 2013). However, the results could not reflect the entire process of reciprocating motion due to tidal currents, even though the hydrodynamic characteristics during ebbing and flooding were also studied. Lefebvre et al. (2014) investigated the dynamic behavior of the bedform-influenced flow field and roughness lengths over a large bedform during a tidal cycle by establishment of a small-scale model with more than 50 layers in the vertical direction. However, their study essentially simulated unidirectional flow in ebbing and flooding, respectively, and the hydrodynamic characteristics for transition were not studied. In this study, we aimed at obtaining the hydrodynamic characteristics of tidal currents for submarine excavation projects. The finite volume method was used to solve the Navier-Stokes equations and the Reynolds stress transport equation. With the foundation trench excavation in the Hong Kong-Zhuhai-Macau Bridge Project in the Pearl River Estuary as an example, the processes of tidal movement after the foundation trench excavation were simulated, and the numerical results were compared with the observed data for the purposes of analyzing the reasonability of the numerical calculation method used in this study.
The basic governing equations for mass and momentum conservation (Pope, 2000; Zhao et al., 2013) are as follows: ∂ρ ∂ (1) + ( ρ ui ) = 0 ∂t ∂xi ∂ui − ρ ui′u ′j + Si (2) µ ∂x j where ui is the velocity component in the x i direction, and i = 1, 2, 3 ; p is the pressure; − ρ ui′u ′j is the Reynolds stress; t is time; ρ is the density of fluid; µ is the dynamic viscosity coefficient; and S i is the momentum source term. Considering anisotropics of turbulent viscosity, the Reynolds stress model (Versteeg and Malalasekera, 2007) was used in this study. The Reynolds stress transport equation and turbulent dissipation rate transport equation (Pope, 2000; Duan et al., 2015) are shown below: ∂ ∂ ρ u i′u ′j + ρu r u i′u ′j = DT ,ij + D L,ij + Pij + G ij + φ ij + ε ij + Fij + S RS (3) ∂t ∂x r
(
AC C
EP
∂ ∂ ∂p ∂ ( ρ ui ) + ( ρ ui u j ) = − + ∂t ∂x j ∂xi ∂x j
)
(
∂ ∂ ∂ ( ρε ) + ( ρε ui ) = ∂t ∂xi ∂xi
)
µt µ + σε
∂ε 1 ε ε2 ρ + C P + C G − C + Sε ( ) ε1 ii ε 3 ii ε2 k k ∂xi 2
(4)
where DT,ij is the turbulent diffusion term; DL,ij is the molecular diffusion term; Pij is the stress production term; Gij is the buoyancy production term; φij is the pressure strain term; ε ij is the dissipation term; Fij is the production term according to the system rotation; k is the turbulent kinetic energy; ε is the turbulent dissipation rate; µ t is the turbulent viscosity; S RS and S ε are the source terms; σ ε is the Prandtl number 2
for ε with a value of 1.3; and C ε 1ACCEPTED , C ε 2 , and C ε 3MANUSCRIPT are constants with the values of 1.44, 1.92, and 0, respectively. DL,ij , Pij , and Fij only contain two-order correlation terms, and DT,ij . Gij , φij , and ε ij contain uncertain correlation terms, which can be calculated as follows: ∂ µ t ∂ui′u ′j DT,ij = (5) ∂xr σ k ∂xr Gij = 0 (6) 32 ε 2 1 ε 3 3 Ck − φij = −C1 ρ ui′u ′j − kδ ij − C2 Pij − Prr δ ij + C1′ρ ur′um′ nr nmδ ij − ui′ur′ n j nr − u ′j ur′ ni nr k 3 3 k 2 2 εd
M AN U
SC
RI PT
(7) Ck 3 2 1 3 1 3 1 C2′ C2 Prm − Prr δ rm nr nmδ ij − C2 Pir − Prr δ ir n j nr − C2 Pjr − Prr δ jr ni nr 3 2 3 2 3 εd 2 ε ij = ρεδ ij (8) 3 where σ k is the Prandtl number for k with a value of 1.0, δ ij is the Kronecker function, n r is the unit normal vector, and d is the distance to the wall. The values of C1 , C2 , C1′ , and C2′ are 1.8, 0.6, 0.5, and 0.3, respectively; and C = Cµ3 4 κ , where κ is the von Kármán constant with a value of 0.4187, and C µ = 0.09 . The air-water interface was considered in this study, and the volume of fluid (VOF) method was used to track the free surface in the model (Rui et al., 2015). In a cell, Fw and Fa are the volume fractions of the water phase and air phase, respectively, and 0 ≤ Fw (or Fa ) ≤ 1 . They are in the relation as follows: Fw + Fa = 1 (9)
2.2. Calculation of hydrodynamic characteristics
(
)
ε
TE D
Drastic turbulence takes place in the vicinity of submarine projects. The hydrodynamic characteristics can be analyzed by calculating the characteristic parameters (Launder and Spalding, 1974), including the turbulent kinetic energy k , the turbulent viscosity µ t , the Reynolds shear stresses τ xy , τ xz , and τ yz , and the wall shear stress τ w , as follows: 1 k = u1′u1′ + u 2′ u 2′ + u 3′ u 3′ (10) 2 k2 µ t = ρ Cµ (11)
(
12 ρκ uWA Cµ1 4 kWA
12 ln E ρ Cµ1 4 k WA d WA µ
)
(12)
(13)
AC C
τw =
EP
τ xy = τ yx = − ρ u1′u 2′ τ xz = τ zx = − ρ u1′u 3′ τ yz = τ zy = − ρ u 2′ u 3′
where uWA , k WA , and d WA are the mean velocity, the turbulent kinetic energy, and the distance to the wall of the wall adjacent cell, and E is an empirical constant with a value of 9.793.
2.3. Numerical calculation scheme 2.3.1. Numerical solution method
The finite volume method was used to solve Reynolds-averaged Navier-Stokes equations, combined with the Reynolds stress model (Pope, 2000). Considering the characteristics of tidal currents and computational efficiency, the pressure-implicit with splitting of operators (PISO) algorithm was used to calculate the pressure and velocity (Jang et al., 1986). The PISO algorithm is based on the higher degree of the approximate relation between the corrections for pressure and velocity. It can maintain a stable calculation with a larger time step for both momentum and pressure (Issa, 1986). The discretization schemes for the diffusion term and the convection term were carried out using a central difference scheme and a second-order wind scheme, respectively. The discretization schemes for the turbulent 3
kinetic energy, turbulent dissipation rate, and Reynolds stress were implemented using a first-order upwind ACCEPTED MANUSCRIPT scheme. The numerical oscillation was limited because the upwind scheme could maintain the stability of numerical calculation. The stability condition (Courant et al., 1967) in the calculation process is as follows:
CFL =
∑λ
i max
i
2V
Ai
≤1
(14)
where CFL is the Courant number, V is the cell volume, Ai is the area of cell face i, and λi max is the maximum of local eigenvalues.
RI PT
2.3.2. Mesh generation
AC C
EP
TE D
M AN U
SC
With the foundation trench excavation in the Hong Kong-Zhuhai-Macau Bridge Project in the Pearl River Estuary as an example, hybrid meshes were generated in the calculation domain. The mesh in the vertical direction was densified at a high resolution in a local region. Unstructured meshes were generated to fit the irregular shape of the structural sections and drastic terrain. Structured meshes were generated in other water areas and at the wall boundary. The meshes near the bottom and free surface were densified. Fig. 1 and Fig. 2(a) present the calculation domains and meshes before and after the foundation trench excavation. The wedge mesh was used for the area of the foundation trench, and the hexahedral mesh was used for the other parts. Along the vertical direction, the meshes near the bottom and free surface were densified with a minimum layer thickness of 0.05 m. According to the observed hydrologic data, there is a mainly reciprocating current in the area between the outlets of the Lingding Estuary and Lantau Island. Therefore, the x direction of the calculation domain was set along the north-to-south direction. The width of the foundation trench with its baseline perpendicular to the flow direction was 50 m. A composite slope was set on both sides with lower segment and upper segment gradients of 1:3 and 1:7, respectively (Cao and Xiao, 2012). Its bottom elevation was −45 m, and the bottom width along the flow direction was 42 m. The size and meshes of the foundation trench after excavation can be seen in Fig. 2(b).
Fig. 1. Schematic diagram of calculation domain and mesh generation before foundation trench excavation
4
RI PT
ACCEPTED MANUSCRIPT
Fig. 2. Schematic diagrams of calculation domain and mesh generation after foundation trench excavation
2.3.3. Boundary conditions
M AN U
SC
The varying processes of velocity and water level were set as the upstream and downstream boundary conditions, respectively. The x direction of the domain was set to be along the dominant flow direction, so the lateral boundaries along the x direction were set as symmetry boundaries. The other boundaries were set as wall boundaries. The time series of the velocity and water level data were obtained from the simulated results of the tidal currents in the Pearl River Estuary using a large-scale model. On the water level boundary, the pressure is in accordance with the assumption of static pressure distribution. It can be calculated with Ps = Pa + ρ g (η − z ) (15) where Ps is the absolute pressure on the water level boundary, Pa is the standard atmospheric pressure, g is the gravitational acceleration, η is the position of the free surface, and z is the vertical coordinate value. The volume fractions of cells on the upstream boundary were obtained from adjacent cells at the previous time step. For the water phase, the volume fraction for surface cell i on the upstream boundary is as follows: FwU ( i, t ′ ) = Fwat ′−1 (16)
EP
F ( i, t ′ ) = 0 wD FwD (i, t ′) = 1 F (i, t ′) = 1 − ( z0 − hwl ) wD 2 ∆z
TE D
where t ′ is the time step, and Fwat ′−1 is the volume fraction of the water phase of the adjacent cell at the previous time step. FwU (i,0) can be determined by the initial flow field. The volume fractions of cells on the downstream boundary were obtained based on the water level process. For the water phase, if ( x0 , y0 , z0 ) is the centroid of the free surface on the downstream boundary, the volume fraction for surface cell i on the downstream boundary is determined as follows:
zi > z0
(17)
AC C
zi < z0
zi = z0
where z i is the vertical coordinate of the centroid of cell i , hwl is the water level at the current time step, and ∆z is the grid spacing in the z direction.
3. Application of numerical calculation method 3.1. Validation In this study, the precision of the numerical calculation method was validated with the observed data at an observation station in the Lingdingyang Estuary under natural conditions. The observed data were obtained by one-day observation (from 12:00 on August 6 to 12:00 on August 7, 2008) using an acoustic wave and current meter. The corresponding monitoring point was set at the center of the calculation domain, located at (505, −22.5).
5
ACCEPTED MANUSCRIPT
3.1.1. Validation of water level
M AN U
SC
RI PT
The water level process calculated at the monitoring point matches the observed values, as shown in Fig. 3(a). The relative error is less than 0.4%. For the numerical results, the highest water level is 11.90 m and the lowest is 10.20 m, which are the same as the observed values. During flooding and ebbing, the calculated values are slightly higher than the observed values, indicating that the water level process obtained from the numerical simulation is reasonable.
Fig. 3. Comparisons of calculated and observed values of water level and depth-averaged velocity at monitoring point
3.1.2. Validation of depth-averaged velocity
EP
TE D
The time series of the calculated results of the depth-averaged velocity at the monitoring point match the observed values, as shown in Fig. 3(b). The calculated and observed results have the same period, and their maximum and minimum values appear at the same time. The average values of the calculated and observed results of the depth-averaged velocity during flooding (from 20:00 on August 6 to 01:00 on August 7) are 0.35 m/s and 0.39 m/s, respectively, with a relative error of 10.26%. The average values of calculated and observed results of the depth-average velocity during ebbing (from 01:00 to 07:00 on August 7) are 0.44 m/s and 0.51 m/s, respectively, with a relative error of 13.73%. At maximum flooding (at 00:00 on August 7), the calculated and observed values are 0.51 m/s and 0.55 m/s, respectively, with a relative error of 7.27%. At maximum ebbing (at 06:00 on August 7), the calculated and observed values are 0.63 m/s and 0.71 m/s, respectively, with a
AC C
relative error of 11.27%. Fig. 3(b) shows that the calculated results are slightly lower than the observed values. The error is caused by the discharge boundary data obtained from the calculated results of the Pearl River Estuary with the large-scale tidal current model. 3.1.3. Validation of vertical distribution of flow velocity Fig. 4 shows the vertical distributions of flow velocity (u) at the monitoring point at different times during flooding and ebbing. Figs. 4(a) through (d) show the velocity distributions during flooding, and Figs. 4(e) through (h) show the velocity distributions during ebbing. In the simulation, the changes in temperature and salinity were not considered. The distributions follow the logarithmic law (Ni et al., 2009), and the calculated results show errors as compared with the observed values, which are caused by the stratification and mixing phenomena. At 22:00 and 23:00 on August 6 and 02:00 and 03:00 on August 7, the distributions match the observed data, and the stratification is weak. However, when the stratification strengthens, the deviations between the calculated results and observed data grow. Moreover, the deviations during ebbing are relatively small.
6
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 4. Calculated results and observed values of vertical distributions of flow velocity at monitoring point
3.2. Engineering application
3.2.1. Hydrodynamic characteristics under natural conditions
EP
TE D
During a tidal cycle, the duration of ebbing is longer than that of flooding, and the water level rises slowly during flooding and drops rapidly during ebbing. Fig. 3 shows that the time of the maximum value of water level lags behind that of the maximum value of flow velocity during a tidal cycle. During flooding, the calculated maximum value of flow velocity is 0.51 m/s, occurring before the high tide level. In the initial stage of flooding, the flow velocity continues to decrease and reaches 0 after 21:00 on August 6. This is due to land runoff and the inertial force of ebbing. Then, it reverses and begins to increase. During ebbing, the calculated
AC C
maximum value of flow velocity is 0.63 m/s, occurring before the low tide level. In the initial stage of ebbing, the flow velocity continues to decrease and reaches 0 after 01:00 on August 7. Then, it reverses and begins to increase. The transition of flow velocity of tidal currents in the lower layer precedes that in the upper layer (Fig. 5), due to the inertia of water particle movement (Xie et al., 2011). In general, the flow velocity in the upper layer is larger than that in the lower layer. It is easier for water particles to overcome the inertia in the lower layer than in the upper layer. Therefore, it is easier for flow to reverse its direction in the lower layer.
7
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Fig. 5. Vertical distributions of flow velocity at monitoring point during transition
3.2.2. Hydrodynamic characteristics after foundation trench excavation
AC C
EP
TE D
After the foundation trench excavation, the area of flow section increases, and the flow slows down through the foundation trench. Fig. 6 presents the velocity distributions in the section of y = −22.5 m. The vertical distributions of flow velocity show wedge-shaped distributions above the lee slope. The velocities at the monitoring point decrease by 60% and 52% at maximum flooding and maximum ebbing after excavation, respectively, in contrast to those before excavation.
Fig. 6. Velocity distributions in section of y = −22.5 m at maximum flooding and maximum ebbing
Fig. 7 presents the flow fields in the section of y = −22.5 m in the foundation trench, with Figs. 7(a) through (f) showing the flow fields during flooding, and Figs. 7(g) through (l) showing the flow fields during ebbing. It indicates that at 22:00 on August 6 and 02:00 on August 7, the flow slows down without swirling flow
8
TE D
M AN U
SC
RI PT
occurring during the transition. Then,ACCEPTED swirling flow occurs due to the increasing velocity after the transition. MANUSCRIPT Furthermore, reverse flow mainly occurs on the lee slope. When the reverse flow moves upward along the lee slope, dammed water occurs above the foundation trench. The maximum change of water level is 0.036 m according to the water level difference before and after excavation. It can be concluded from Fig. 8 that the excavation has little influence on the water level. The results of flow fields match the results of Cao and Xiao (2012), which were obtained based on a two-dimensional horizontal numerical model and field tests.
AC C
EP
Fig. 7. Flow fields in section of y = −22.5 m in foundation trench
Fig. 8. Variations in water level before and after excavation
Intense turbulence takes place in the foundation trench due to the occurrence of swirling flow. According to Eqs. (10) through (13), the hydrodynamic characteristics were obtained. The maximum turbulent viscosities at maximum flooding and maximum ebbing are 100 kg/(m·s) after foundation trench excavation, while the corresponding values under natural conditions are 13.6 kg/(m·s) and 55.8 kg/(m·s), respectively, indicating that the excavation of the foundation trench intensifies the turbulence in this area. The comparison of turbulent characteristics in different tidal stages shows that the turbulent flow is stronger during ebbing than in other tidal stages. As shown in Fig. 9, the range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the strong turbulent motion mainly takes place in the water above the stoss slope. However, the turbulent viscosity decreases during the transition, as shown in Figs. 9(c) and (d). 9
RI PT
ACCEPTED MANUSCRIPT
Fig. 9. Distributions of turbulent viscosity µ t at maximum flooding, maximum ebbing, and trannsition in section of y = −22.5 m
EP
TE D
M AN U
SC
Fig. 10 shows that the values of Reynolds shear stress τ xz at maximum flooding and maximum ebbing are larger than those during the transition, and the maximum value of τ xz occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle, demonstrating that there exists intensive exchange of water mass in different layers, and the shear velocity shows a large degree of variation due to the occurrence of eddies.
Fig. 10. Distributions of Reynolds shear stress τ xz at maximum flooding, maximum ebbing, and transition in section of y = −22.5 m
AC C
Fig. 11 presents the wall shear stress distribution along the flow direction at maximum flooding, maximum ebbing, and the transition. The values of wall shear stress at maximum ebbing are largest due to the larger values of flow velocity. During the transition, the values of wall shear stress are small. Furthermore, the wall shear stress is higher on the north side, implying that the frictional resistance is larger on the north side, and the sediment incipient motion and transport occur more readily on the south side. The results match the phenomenon of the thickness of deposition on the north side being greater than that on the south side. Also, the hydrodynamic characteristics obtained in this study agree with a previous study by Xin et al. (2012), in which the hydrodynamic characteristics were obtained based on trial tunnel experiments over a long-period observation.
10
RI PT
ACCEPTED MANUSCRIPT
Fig. 11. Wall shear stress distribution along flow direction at maximum flooding, maximum ebbing, and transition
4. Conclusions
References
TE D
M AN U
SC
Based on the Navier-Stokes equations and the Reynolds stress model, a three-dimensional small-scale numerical model of tidal currents for submarine excavation projects was established. The results match the observed data and show a strong performance of the numerical model in simulation of the hydrodynamic characteristics of tidal currents. The following conclusions are drawn: (1) After the foundation trench excavation, the flow velocity decreases quite a bit through the foundation trench, and reverse flow occurs on the lee slope in the foundation trench. With the action of swirling flow, the reverse flow moves upward along the lee slope, and the inflow is impeded, leading to the occurrence of dammed water above the foundation trench. However, the excavation has little influence on the water level. (2) The turbulent motion is stronger during ebbing than in other tidal stages. The range with the maximum value of turbulent viscosity, occurring on the south side of the foundation trench at maximum ebbing, is greater than those in other tidal stages in a tidal cycle, and the maximum value of Reynolds shear stress occurs on the south side of the foundation trench at maximum ebbing in a tidal cycle. The hydrodynamic characteristics obtained in this study demonstrate the phenomenon of the thickness of deposition on the north side being greater than that on the south side. The flow field characteristics and hydrodynamic characteristics of tidal currents around the submarine excavation project are reasonably described by the numerical calculation method.
AC C
EP
Abbott, M.B., 1997. Range of tidal flow modeling. Journal of Hydraulic Engineering 123(4), 257-277. http://dx.doi.org/10.1061/(ASCE)0733-9429(1997)123:4(257). Berntsen, J., Furnes, G., 2005. Internal pressure errors in sigma-coordinate ocean models: Sensitivity of the growth of the flow to the time stepping method and possible non-hydrostatic effects. Continental Shelf Research 25(7-8), 829-848. http://dx.doi.org/ 10.1016/j.csr.2004.09.025. Berntsen, J., Thiem, Ø., Avlesen, H., 2015. Internal pressure gradient errors in σ -coordinate ocean models in high resolution fjord studies. Ocean Modelling 92, 42-55. http://dx.doi.org/10.1016/j.ocemod.2015.05.009. Cao, H.J., Xiao, L.B., 2012. Back silting of trench excavation for the immersed tunnel in Hong Kong-Zhuhai-Macao Bridge construction. Port and Waterway Engineering (1), 12-17 (in Chinese). Chen, C.S., Huang, H.S., Beardsley, R.C., Liu, H.D., Xu, Q.C., Cowles, G., 2007. A finite volume numerical approach for coastal ocean circulation studies: Comparisons with finite difference models. Journal of Geophysical Research 112(C03018), 1-34. http://dx.doi.org/10.1029/2006JC003485. Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11(2), 215-234. http://dx.doi.org/10.1147/rd.112.0215. Duan, D.H., Zhang, J.B., Gao, Z.X., 2015. Investigation on spanwise wall oscillation in turbulent channel flow based on Reynolds stress model (RSM). Procedia Engineering 99, 665-670. http://dx.doi.org/10.1016/j.proeng.2014.12.587 Gao, J.S., Xue, H.J., Chai, F., and Shi, M.C. 2013. Modeling the circulation in the gulf of Tonkin, South China Sea. Ocean Dynamics 63(8), 979-993. http://dx.doi.org/10.1007/s10236-013-0636-y. Gräwe, U., Holtermann, P., Klingbeil, K., Burchard, H., 2015. Advantages of vertically adaptive coordinates in numerical models of stratified shelf seas. Ocean Modelling 92, 56-68. http://dx.doi.org/10.1016/j.ocemod.2015.05.008. Issa, R.I., 1986. Solution of implicitly discretized fluid flow equations by operator splitting. Journal of Computational Physics 62, 11
AC C
EP
TE D
M AN U
SC
RI PT
40-65. http://dx.doi.org/10.1016/0021-9991(86)90099-9. ACCEPTED MANUSCRIPT Jang, D.S., Jetli, R., Acharya, S. 1986. Comparison of the PISO, SIMPLER, and SIMPLEC algorithms for the treatment of the pressure-velocity coupling in steady flow problems. Numerical Heat Transfer 10(3), 209-228. http://dx.doi.org/10.1080/10407798608552506. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3(2), 269-289. http://dx.doi.org/10.1016/0045-7825(74)90029-2. Lefebvre, A., Paarlberg, A.J., Ernstsen, V.B., Winter, C., 2014. Flow separation and roughness lengths over large bedforms in a tidal environment: A numerical investigation. Continental Shelf Research 91, 57-69. http://dx.doi.org/10.1016/j.csr.2014.09.001. Li, Z.L., 2015. Three times installation of 15th immersed tube in Hong Kong-Zhuhai-Macau submarine tunnel engineering in Lingdingyang Bay. Construction and Architecture (19), 13-15 (in Chinese). Liu, H., He, Y.S., 2000. Recent development in 3-D hydrodynamic mathematical model studies for estuaries. The Ocean Engineering 18(2), 87-93 (in Chinese). Lu, W.Q., Wu, W.G., Su, L.W., Chen, K.Q., Peng, S., 2014. Hydrodynamic experimental study on immersed tube tunnel element of Hong Kong-Zhuhai-Macao Bridge. China Civil Engineering Journal 47(3), 138-144 (in Chinese). Lubin, P., Glockner, S., Chanson, H., 2010. Numerical simulation of a weak breaking tidal bore. Mechanics Research Communications 37(1), 119-121. http://dx.doi.org/10.1016/j.mechrescom.2009.09.008. Ni, Z.H., Song, Z.Y., Wu, L.C., 2009. Study on the double-logarithmic profile of tidal flow velocity in the near-bed layers. Acta Oceanologica Sinica 28(6), 84-92. Nielsen, A.W., Liu, X.F., Sumer, B.M., Fredsøe, J., 2013. Flow and bed shear stresses in scour protections around a pile in a current. Coastal Engineering 72, 20-38. http://dx.doi.org/10.1016/j.coastaleng.2012.09.001. Ningsih, N.S., Azhar, M.A., 2013. Modelling of hydrodynamic circulation in Benoa Bay, Bali. Journal of Marine Science and Technology 18(2), 203-212. http://dx.doi.org/10.1007/s00773-012-0195-9. Ohashi, K., Sheng, J.Y., 2013. Influence of St. Lawrence River discharge on the circulation and hydrography in Canadian Atlantic waters. Continental Shelf Research 58, 32-49. http://dx.doi.org/10.1016/j.csr.2013.03.005. Pope, S.B., 2000. Turbulent Flow. Cambridge University Press, Cambridge, pp. 83-95. Ren, L., Nash, S., Hartnett, M., 2015. Observation and modeling of tide- and wind-induced surface currents in Galway Bay. Water Science and Engineering 8(4), 345-352. http://dx.doi.org/10.1016/j.wse.2015.12.001. Rui, B.F., Falcão, D.S., Oliveira, V.B., Pinto, A.M.F.R., 2015. Numerical simulations of two-phase flow in proton exchange membrane fuel cells using the volume of fluid method: A review. Journal of Power Sources 277, 329-342. http://dx.doi.org/10.1016/j.jpowsour.2014.11.124. Shi, J.Z., Li, C., Dou, X.P., 2010. Three-dimensional modeling of tidal circulation within the north and south passages of the partially-mixed Changjiang River Estuary, China. Journal of Hydrodynamics, Ser. B 22(5), 656-661. http://dx.doi.org/10.1016/S1001-6058(10)60010-8. Tang, G.Q., Chen, C.Q., Zhao, M., Lu, L., 2015. Numerical simulation of flow past twin near-wall circular cylinders in tandem arrangement at low Reynolds number. Water Science and Engineering 8(4), 315-325. http://dx.doi.org/10.1016/j.wse.2015.06.002. Tang, H.S., Qu, K., Wu, X.G., 2014. An overset grid method for integration of fully 3D fluid dynamics and geophysics fluid dynamics models to simulate multiphysics coastal ocean flows. Journal of Computational Physics 273, 548-571. http://dx.doi.org/10.1016/j.jcp.2014.05.010. Tao, J.F., Zhang, C.K., 2007. Review on several vertical coordinate modes of 3-D flow numerical model in estuarine and coastal waters. The Ocean Engineering 25(1), 133-142 (in Chinese). Tuomi, L., Kai, M., Lehmann, A., 2012. The performance of the parameterisations of vertical turbulence in the 3D modelling of hydrodynamics in the Baltic Sea. Continental Shelf Research 50-51, 64-79. http://dx.doi.org/10.1016/j.csr.2012.08.007. Versteeg, H.K., Malalasekera, W., 2007. An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Second ed. Pearson Education Limited, Harlow. Wu, Z.C., Wang, D.Z., 2009. Numerical solution for tidal flow in opening channel using combined MacCormack-finite analysis method. Journal of Hydrodynamics, Ser. B 21(4), 505-511. http://dx.doi.org/10.1016/S1001-6058(08)60177-8. Xie, D.F., Pan, C.H., Wu, X.G., 2011. Three-dimensional numerical modeling of tidal bore in Qiantang based on FVCOM. The Ocean Engineering 29(1), 47-52 (in Chinese). Xin, W.J., Jia, Y.S., He, J., 2012. Back silting analysis of the trial dredged-trough for immersed tube tunnel of Hong Kong-Zhuhai-Macao Bridge. Hydro-Science and Engineering (2), 71-78 (in Chinese). Zhao, Y.P., Bi, C.W., Dong, G.H., Gui, F.K., Cui, Y., Guan, C.T., Xu, T.J., 2013. Numerical simulation of the flow around fishing plane nets using the porous media model. Ocean Engineering 62, 25-37. http://dx.doi.org/10.1016/j.oceaneng.2013.01.009. Zhou, J.T., Falconer, R.A., Lin, B.L., 2014. Refinements to the EFDC model for predicting the hydro-environmental impacts of a barrage across the Severn Estuary. Renewable Energy 62, 490-505. http://dx.doi.org/10.1016/j.renene.2013.08.012.
12