# Study on Impact of Turbine Location on Hydrodynamics in Tidal Farm

Authors: Hongqiang Zhang; Daming Li; Yanqing Li; Ting Yang; Shan Luo; Shunfa Tian; Shilong Bu
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

This work aims at investigating the impact of different tidal turbine locations on hydrodynamics in near-field and far-field flow; thus, three- and two-dimensional (3- & 2-D) models were exploited in combination and applied in a case study of Putuo-Hulu Islands’ tidal farm. We present a method for the simplification of tidal turbine which, based on the energy equation, determines turbine’s equivalent roughness by calculating resistance loss in flow passage. A 3-D turbine model of near-field flow in the frame of Computational Fluid Dynamics was constructed, and the flow velocity distributions in 7 combinations of “impeller rotating speed - inflow velocity” were simulated. Also a 2-D tidal model of far-field flow was established, and Finite Element Method was adopted to solve the 2-D shallow water circulation equations; thus, the impact of tidal turbine could be simulated by utilizing different location’s compositive roughness. The results show that the impact of turbine location on hydrodynamics is depending on the open degree of sea area, channel trend, and bathymetric and geographic features and that the farther distance from the turbine, the less impact on the flow field. Overall, the impact of turbine location on far-field flow is not significant, and the flow velocity varies below 8% relative to the velocity prior to turbine installation.


---

## Body

## 1. Introduction

Tidal current power is a green and renewable energy, and the prospect of its development is considerable. Many countries with abundant marine resources are promoting the research on tidal current power generation technology and building up test field of tidal current generation.In near-field hydrodynamic study of tidal turbine, Betz’ Law and Blade Element Momentum (BEM) theory were adopted to generalize the turbine. Bai et al. [1] calculated the eddy effect of the tidal turbine by Computational Fluid Dynamics (CFD) method based on Betz’ Law. O’Doherty et al. [2] simulated the effect of tidal turbine on hydrodynamic characteristics of flow field under different arrangement modes.In far-field hydrodynamic study of tidal turbine, momentum method was used to generalize tidal turbine, that is, to add resistance loss item of turbine to the dynamic control equation. Couch and Bryden [3] added the bottom friction coefficient to a one-dimensional (1-D) model to simulate the effects of tidal turbine. Bluden and Bahaj [4, 5] added the bottom friction coefficient in the model and utilized it to simulate the array of tidal turbine; the results showed that the effect of array on initial flow rate is about 5% to 10%. Bryden and Couch [6] established a 1-D model to study the potential impact of tidal current development on the local tidal current field under steady conditions. Karsten et al. [7] constructed a two-dimensional (2-D) model to study the tidal current field of Fendi Bay by increasing the bottom friction. Myers and Bahaj [8] replaced turbines with perforated discs, and studies showed that the test method is feasible. Geórgenes et al. [9] investigated the relationship between tidal distortion and salinity by utilizing classical methods of comparison of three cross-channel circulation characteristics. Gaurav and Tate [10] applied the adaptive hydraulics numerical code to study tidal propagation in Lower Columbia River estuary. Jeyar et al. investigated the hydrodynamic impacts of whole lagoon due to tidal waves using a 2-D, depth-averaged hydrodynamic model based on the shallow water equations. Gillibrand et al. [11] adopted a three-dimensional (3-D) hydrodynamic model to consider the potential effects of energy extraction by an array of tidal turbines on ambient near-bed velocity field and local bed shear stress. Yamamoto et al. [12] exploited ANSYS to select the best shape of oscillator, in order to efficiently harness the tidal flow energy in a cost-efficient manner.Although the study of tidal current power started late in China, the scientists also made some new breakthroughs in the research and development of energy-harvesting device, and the study of flow field also has a lot of achievements. For example, Xin et al. [13] predicted the hydrodynamic characteristics of double rotor turbines arranged in series. An [14] made the 3-D steady-state numerical simulation of the horizontal axis turbine by CFD method, and mutual effect rule between devices under different arrangement conditions was obtained. Yuan [15] utilized a perforated disc to simulate the tidal turbine by increasing drag coefficient to replace device group. Hou [16] exploited Finite Volume Community Ocean Model to simulate the effect of tidal turbine on 2-D hydrodynamic environment. Li [17] used CFD method to study reasonable arrangement distance of the array of horizontal axis tidal turbine. Ma [18] established blade model according to BEM drawing on the experience of wind energy generation and applied CFD method in studying the hydrodynamic characteristics of the impellers under different working conditions. Ding et al. [19] found local surge current velocity in each depth, with a magnitude of the same order as astronomic tidal currents.At home and abroad, studies on the effect of tidal turbine on hydrodynamic characteristics of large-scale far-fields which are combined with the integrated simulation of near-field hydrodynamic characteristics of turbine are few, especially in the aspects of turbine simplification and bottom friction determination. But in the previous studies, generalized bottom frictions were fixed values and did not change with the variation of flow velocity. Due to the large subjectivity of the existing research methods and the poor adaptability of the simplification method under different working conditions, the research methods could be developed in the direction of “near-field & far-field” integrated simulation. Therefore, in order to clarify the effect on flow field because of turbine installation, and then rationally arrange the berthages of tidal turbine in test field, “Putuo-Hulu Islands Tidal Power Demonstration Project” in Zhoushan City, China, was regarded as the study object. Firstly, we constructed a 3-D turbine model of near-field flow in the frame of CFD and simulated the flow velocity distributions in 7 combinations of “impeller rotating speed - inflow velocity”. Next, we presented a simplification method of tidal turbine, based on the energy equation, to determine turbine’s equivalent roughness by calculating resistance loss in flow passage. Then, we established a 2-D tidal model of far-field flow and applied Finite Element Method (FEM) in solving the 2-D shallow water circulation equations. Finally, the impact of tidal location on hydrodynamics could be simulated by utilizing different location’s compositive roughness. We combined 3-D turbine model and 2-D tidal model and simulated and analyzed the impact of turbine location on near-field and far-field hydrodynamics in tidal farm.


## 2. Study Area

### 2.1. Tidal Farm Introduction

The site of project is located in the channel between Putuo Island and Hulu Island in Zhoushan City, China, and the geographical location of the tidal farm is shown in Figure1. The channel is about 6.6 km long and 1.6 km wide, and the water depth is about 20 m to 80 m (1985 National Elevation Benchmark). The offshore topography on both sides of channel is steeper, and the central site is relatively flat.Figure 1
Geographical location of tidal farm.Putuo Island and Hulu Island lie in a sea area where islands are numerous, and water depth changes complexly. In the sea area, the tidal current is characterized by the irregular semidiurnal tide, the fluctuation of tidal current is obviously affected by topography, and shallow sea effect is obvious. Velocity of tidal current between islands is relatively large, the mainly flow pattern is reciprocating flow, and the direction is roughly parallel to the channel. Meanwhile, the direction of tidal current is basically consistent with the tidal wave propagation [20]. In the channel, the tidal current is reciprocating flow along NW-SE direction, the vertical average velocity is about 0.67 m/s to 0.81 m/s, and the maximum velocity of tidal current on surface is 2.29 m/s. The channel is a sub-rich area of tidal current resources, and under the cover of Hulu Island, the extreme sea conditions of open sea have little effect on channel.


### 2.2. Tidal Turbine Location

Based on the compositive analysis of the bathymetric data, engineering geology, and tidal current resource conditions in this sea area, combining with the existing measured data, it is proposed to set the measuring points B1 to B5 as the berthages of the tidal turbine. The tidal turbine location in tidal farm is shown in Figure2. The white parts in Figure 2 represent the islands or land (i.e., correspond to Putuo Island and Hulu Island in Figure 1), and the rest represents the sea area. Figure 2 is an enlarged view of the red frame of Figure 1.Figure 2
Tidal turbine location.


## 2.1. Tidal Farm Introduction

The site of project is located in the channel between Putuo Island and Hulu Island in Zhoushan City, China, and the geographical location of the tidal farm is shown in Figure1. The channel is about 6.6 km long and 1.6 km wide, and the water depth is about 20 m to 80 m (1985 National Elevation Benchmark). The offshore topography on both sides of channel is steeper, and the central site is relatively flat.Figure 1
Geographical location of tidal farm.Putuo Island and Hulu Island lie in a sea area where islands are numerous, and water depth changes complexly. In the sea area, the tidal current is characterized by the irregular semidiurnal tide, the fluctuation of tidal current is obviously affected by topography, and shallow sea effect is obvious. Velocity of tidal current between islands is relatively large, the mainly flow pattern is reciprocating flow, and the direction is roughly parallel to the channel. Meanwhile, the direction of tidal current is basically consistent with the tidal wave propagation [20]. In the channel, the tidal current is reciprocating flow along NW-SE direction, the vertical average velocity is about 0.67 m/s to 0.81 m/s, and the maximum velocity of tidal current on surface is 2.29 m/s. The channel is a sub-rich area of tidal current resources, and under the cover of Hulu Island, the extreme sea conditions of open sea have little effect on channel.


## 2.2. Tidal Turbine Location

Based on the compositive analysis of the bathymetric data, engineering geology, and tidal current resource conditions in this sea area, combining with the existing measured data, it is proposed to set the measuring points B1 to B5 as the berthages of the tidal turbine. The tidal turbine location in tidal farm is shown in Figure2. The white parts in Figure 2 represent the islands or land (i.e., correspond to Putuo Island and Hulu Island in Figure 1), and the rest represents the sea area. Figure 2 is an enlarged view of the red frame of Figure 1.Figure 2
Tidal turbine location.


## 3. Methodology

### 3.1. Establishment of 3-D Near-Field Hydrodynamic Model

#### 3.1.1. Turbine Model and Grids Generation

There are two different carrier forms of tidal turbine with horizontal axis, variable pitch, and double impellers. One is for column type with single pile; that is, the structure is supported by a single pile and the impeller can be up and down according to the work needs. The other is for floating; the unit is supported by the floating platform, and the entire device is fixed by four anchor chains. Both the two carrier forms are double impellers with a rotary diameter (D0) of 16.3 m, and the equivalent diameter (DE=2D0) [21] of an impeller is about 23 m.The role of the carrier is ignored in the 3- and 2-D numerical simulation, and only the effect of the double impellers on the flow field is considered. The rotation direction of the impellers is opposite, the spacing between impellers is 3 m, and the hub diameter is 0.1D0. Then, the 3-D modeling software is used to simplify and establish the model of double impellers. The turbine model is shown in Figure 3.Figure 3
3-D turbine model.According to the calculation needs, a cylinder flow field is designed, as shown in Figure4(a). The length of the flow field is 300 m (X direction), the diameter is 100 m (Y and Z direction), the longitudinal axis is parallel to the inflow direction, and the center of the double impellers is located in the longitudinal axis which is away from the upstream inlet 60 m and 240 m from the downstream outlet. This 3-D region is divided by unstructured tetrahedral mesh, in which the green parts represent the rotating domains of impellers and other parts are the nonrotating domains, as shown in Figure 4(b). The grids of rotating domains where the impellers are located are refined, and the grids far away from the impellers are sparse. The maximum grid size on the surface of the impellers (i.e., nearby the walls) is 0.02 m, the maximum grid size of the rotating domains is 0.1 m, and the maximum grid size of the nonrotating domains is 0.6 m. And then, a total of 3,489,648 grids are generated.Figure 4
Computational domain and grids: (a) schematic of computational domain and (b) 3-D grids.
(a) (b)The density of grids in the range of 240 m from rotating domain to nonrotating domain in the cylinder flow field decreases gradually, and the maximum size of the grids increases from 0.1 m to 0.6 m. Moreover, there is no too large or too small size mutation. All these conform to the requirements of renormalization group (RNG) k-ε turbulence model for the calculating grids. In the simulation of Fluent, the calculation equations of RNG k-ε turbulence model are as follows:(1)∂∂tρk+∂∂xiρkui=∂∂xjαkμeff∂k∂xj+Gk+Gb-ρε-YM+Sk(2)∂∂tρε+∂∂xiρεui=∂∂xjαεμeff∂ε∂xj+C1εεkGk+C3εGb-C2ερε2k-Rε+Sεwhere k is turbulence kinetic energy, ε is energy dissipating rate, xi is spatial scale, ui is flow velocity, αk and αε are the Prandtl numbers’ reciprocals of effective turbulent of k and ε, respectively, μeff is viscosity of effective vortex, Gk is turbulence kinetic energy generated by average velocity gradient, Gb is turbulence kinetic energy generated by flotage, YM is influence of turbulent pulsation expansion on overall dissipating rate, and Sk and Sε are user-defined source items.This paper selects 3 points which are away from the impeller surface 50 m, 100 m, and 150 m on the axis of impeller in Condition 1 (inflow velocity is 1.843 m/s), and the three points are recorded as Point 1, 2, and 3, respectively. The changes of velocity of the 3 points under different numbers of grids are shown in Figure5.Figure 5
Variations in velocity with grid numbers.From Figure5, it is known that when the grid number is greater than 3 million, the variations of velocity of the 3 points are very small with the increase of grid number, thus indicating that the grid number has little effect on the simulation results. Therefore, when the grid number is higher than 3 million, the simulation results almost do not change with the variations of the grid number, and results are undoubtedly convergent. So the grid density has little effect on the simulation results; that is, the grid number is almost independent of the simulation results.


#### 3.1.2. Relationship between Impeller’s Rotating Speed and Inflow Velocity

Since the turbine model in this paper is constructed by imitating the shape of SeaGen impeller of tidal turbine, its radiusr=8 m, rated flow velocity vr=2.25 m/s, and rated revolving speed of turbine nr=14.3 rad/min. From those, we could obtain that the rated angular speed ωr=2πnr=1.5 rad/s and the design tip speed ratio λ=rωr/vr=5.3. Therefore, the design tip speed ratio of impeller λ is taken as 5.3 to establish the corresponding relationship between the impeller’s rotating speed and the inflow velocity, as shown in Table 1.Table 1
Relationship between impeller’s rotating speed and inflow velocity.
Inflow Velocity (m/s) Parameters Impeller Rotating Position per Second Relationship Schematic between Impeller’s Rotating  Speed/Angular Speed and Inflow Velocity Angular Speed (rad/s) Rotating Speed (r/min) 0.614 0.4 3.82     0.921 0.6 5.73     1.229 0.8 7.64     1.536 1.0 9.55     1.843 1.2 11.46     2.150 1.4 13.37     2.457 1.6 15.28


#### 3.1.3. Solver Setting

The CFD software, Fluent, is applied in numerical simulation, and the specific setting of the CFD solver is as follows:( 1 ) Select the 3-D single-precision solver.( 2 ) Select the pressure-based correction method, and select the steady state for the time type.( 3 ) Select RNG k-ε turbulence model and Multiple Reference Frame (MRF) model.( 4 ) The inflow velocity of the velocity-inlet is vi, because the outlet velocity or pressure cannot be known before the problem is solved. Free outflow is used at the outlet boundary, the side wall is a stationary wall, and the impellers are moving walls.( 5 ) 1000 iterations each time step.


### 3.2. Turbine Simplification

For the flow passage where the turbine is located, the flux of any flow section of this flow passage is equal when it is assumed that there is no flow crossing. According to the continuity equation, the acreage of flow sections of the downstream flow passage will be greater than the upstream because the flow velocity of the downstream of the turbine will decrease, as shown in Figure6.Figure 6
Flow passage parameters.The turbine causes partial energy loss of the water flow; it can be expressed by the frictional resistance losshf between flow sections along upstream to downstream of the flow passage, and the corresponding equivalent roughness nE is obtained, then the roughness n can be obtained by combining the original natural roughness n0 of the grid, so as to realize the simplification of the effect of the turbine on the hydrodynamic force, that is, the local roughness correction method.The cross section of the impeller is taken as the reference plane, two flow sections are taken in the upstream L1 and downstream L2, respectively, and the water body surrounded by the two flow sections and the side wall of the flow tube is studied. It is assumed that the average flow velocity of the flow section 1-1 at a certain time isv1, the pressure is p1, and the average flow velocity of the flow section 2-2 is v2, the pressure is p2. The energy equation of the constant pipe flow of incompressible fluid is established, as shown in (3), and only the resistance loss along the path is considered. In (3), α1=α2=1 and z1=z2=0, because the central axis is regarded as the datum plane.(3)z1+p1ρg+α1v122g=z2+p2ρg+α2v222g+hf1-2And then, (3) can be written as (4).(4)hf1-2=p1-p2ρg+v12-v222gAccording to Darcy Formula,hf1-2 can be expressed as (5)hf1-2=λL4Rv122gThe approximate relationship between the resistance coefficient (λ) and the Chézy Coefficient (C) is shown as (6)λ=8gC2The relationship betweenC and nE is obtained by the Manning Formula, as shown in the following.(7)C=1nER1/6The equivalent roughness (nE) of tidal current power generation equipment is obtained by synthesizing (3) to (7), as shown in (8)nE=hf1-2R4/3Lv12=p1-p2ρg+v12-v222gR4/3Lv12The wall shear stress formula is shown as(9)τ=ρgRJ=ρgv2C2=ρgn2v2R1/6The shear stressτ of the seabed is considered to be composed of two parts. One is the equivalent shear stress Δτ produced by turbine’s layout, the corresponding equivalent roughness is nE, and the flow velocity is v2. The other one is the shear stress of the sea itself τ0 when the flow velocity is v2, and the corresponding roughness of seabed is n0. According to the superposition principle, the shear stress of the seabed can be expressed as (10)τ=τ0+Δτ=ρgRJ=ρgn02v22R1/6+ρgnE2v22R1/6=ρgn2v22R1/6The compositive roughness of the turbine after it is placed is as shown in(11)n=n02+nE2


### 3.3. Establishment of 2-D Far-Field Hydrodynamic Model

#### 3.3.1. Finite Element Discretization of Governing Equation

The 3-D N-S equation is integrated along the depth direction (Z direction) to obtain the 2-D shallow water circulation equation, which is used as the governing equation. Meanwhile, the role of wind is ignored, as shown in the following:(12)∂z∂t+∂uh∂x+∂vh∂y=0(13)∂u∂t+u∂u∂x+v∂u∂y+g∂z∂x=-gu2+v2C2hu+fv+Mx(14)∂v∂t+u∂v∂x+v∂v∂y+g∂z∂y=-gu2+v2C2hv-fu+Mywhere Mx=AH(∂2u/∂x2+∂2u/∂y2), My=AH(∂2v/∂x2+∂2v/∂y2) are vortex viscosity terms, AH is vortex viscosity coefficient and its value ranges from 0.4 to 0.8, z is tide level, h is water depth, z=h+z0, z0 is bottom elevation, t is time, u and v are components of velocity in x and y direction, respectively, C is Chézy coefficient, f is Coriolis force coefficient, and g is gravitational acceleration.In this paper, FEM is used to discretize the 2-D shallow water circulation equation. The computational area is divided by the uniform triangle grids, the number of unstructured grids isM, and the number of nodes is N. In addition, φn (n=1,2,…,N) is set as the basis function of the solution function of the finite element equation, and φpe is set as the shape function of local units. Then the solution of the governing equation can be expressed as (15)ut,x,y=∑n=1Nuntφnx,y=∑e=1M∑r=13φpreurevt,x,y=∑n=1Nvntφnx,y=∑e=1M∑r=13φprevrezt,x,y=∑n=1Nzntφnx,y=∑e=1M∑r=13φprezrewhere the part after the first equal sign is the solution function under the condition that the grid nodes are arranged by ordinal number in the whole coordinate, the part after the second equal sign is the solution function under the condition that the grids are arranged by ordinal number in the local coordinate. The solution of arbitrary units in (15) is given by (16)uet,x,y=∑r=13φprex,yuretvet,x,y=∑r=13φprex,yvretzet,x,y=∑r=13φprex,yzretEquation (16) is substituted into (12) to (14), he(x,y,t)=ze(x,y,t)-z0(x,y), then continuity Equation (17) and motion Equation (18) of triangle grids are obtained by Galerkin method.(17)∫Ω∑r=13φpre∂zre∂t+∑j=13∂φpje∂xuje∑r=13φprehre+∑j=13∂φpje∂xhje∑r=13φpreure+∑j=13∂φpje∂yvje∑r=13φprehre+∑j=13∂φpje∂yhje∑r=13φprevreφqkedΩ=0k=1,2,3(18)∫Ω∑r=13φpre∂ure∂t+∑j=13∂φpje∂xuje∑r=13φpreure+∑j=13∂φpje∂yuje∑r=13φprevre+g∑r=13∂φpre∂xzre+g∑j=13φpjeuje2+∑j=13φpjevje2∑j=13cj2φpjehje∑r=13φpreure-∑r=13frφprevreφqkedΩΩ=0k=1,2,3A r k = ∫ Ω φ p ( r ) e φ q ( k ) e d Ω, Brk=∑j=13(uje∂φpje/∂x+vje∂φpje/∂y)Ark, Crk=∑j=13(hje∂φpje/∂x)Ark, Drk=∑j=13(hje∂φpje/∂y)Ark are substituted into Equation (17), and Erk=∑j=13(uje∂φpje/∂x)Ark, Erk′=∑j=13(vje∂φpje/∂x)Ark, Frk=∑j=13(uje∂φpje/∂y)Ark, Frk′=∑j=13(vje∂φpje/∂y)Ark, Hrk=g∫Ω(∂φpre/∂x)φq(k)edΩ, Hrk′=g∫Ω∂φpre/∂yφq(k)edΩ, Jrk=∫Ωfrφp(r)eφq(k)edΩ, Irk = g∫Ωφpreφqke((∑j=13φpjeuje)2+(∑j=13φpjevje)2/ ∑j=13cj2φpjehje)dΩ are substituted into Equation (18). Then, finite element equation of element matrix Equation (19) is obtained by changing the order of integration summation.(19)Az˙+Bh+Cu+Dv=0Au˙+Eu+Fv+Hz+Iu-Jv=0Av˙+E′u+F′v+H′z+Iv+Ju=0where z˙=∂z1/∂t∂z2/∂t∂z3/∂te, u˙=∂u1/∂t∂u2/∂t∂u3/∂te, v˙=∂v1/∂t∂v2/∂t∂v3/∂te, z=z1z2z3e, h=h1h2h3e, u=u1u2u3e, v=v1v2v3e, A=A11A12A13A21A22A23A31A32A33, B, C, D, E, F, H, I, J, E′, F′, H′ are similar to A, and they are 3×3 matrixes which correspond to the elements corresponding to the unit.After the unit partition is determined, the number of the local unit nodes in whole coordinate is determined to form the whole matrix (increase the subscriptT to represent). Moreover, the time differential is discretized by the forward difference, and the moment to be sought is expressed by the superscript t+Δt, and the remainder is t. Then (20) is obtained.(20)ATzt+Δt=ATzt-ΔtBTht+CTut+DTvtATut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtATvt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtEquation (20) can be rewritten as (21) by using mass concentration method, diagonalizing 3×3 matrix starting from unit matrix, replacing matrix AT at time t+Δt by diagonal matrix A0, and then solving ut+Δt and vt+Δt of the equation and finally ut+Δt and vt+Δt are substituted into the equation zt+Δt.(21)A0ut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtA0vt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtA0zt+Δt=ATzt-ΔtBTht+CTut+Δt+DTvt+Δt


#### 3.3.2. Grids Generation and Submarine Topography

The study area mainly includes the northeastern coast of Putuo Island and Hulu Island, the quadrilateral open sea facing the East China Sea. The study area is about 10 km from west to east, about 8 km from south to north, covering about 74 km2. The unstructured grids are used to divide the computational area in the model, and the density of the grids increases gradually from the open sea area to the island coastline. Finally, the number of unstructured grids’ nodes is 7,970 and the number of unstructured grids is 15,264. The grids generation of tidal model is shown in Figure 7. According to the data, the contour map of submarine elevation is obtained through the transformation of coordinate and datum plane, the interpolation of grid nodes, and a series of processing, as shown in Figure 8.Figure 7
Grids generation of tidal model.Figure 8
Submarine elevation contour.Comparing Figures7 and 8, it could be found that the topography near the open boundary is relatively flat and the simulation of tidal wave is more stable. Therefore, increasing the grid size and reducing the resolution is reasonable but also can improve the model calculation efficiency. Increasing the density of the grid at the waterway and the island coast can improve the accuracy of the operation and fully reflect the submarine topography.The type of the open boundary of the model is the tide boundary; that is, the water level process of each grid node of the determined boundary is given along the open boundary. The water level process is obtained by harmonic analysis method and it is predicted by the global ocean tide model [22], and eight tidal constituents M2, S2, N2, K2, K1, O1, P1, Q1 are taken into account.


#### 3.3.3. Validation of 2-D Model

The marine hydrological data is obtained by the Second Institute of Oceanography of the State Oceanic Administration in August 2013, such as tide and tidal wave. The observation period was 10:00 on August 16th, 2013 to 17:00 on August 17th, 2013 (small tide), 14:00 on August 19th to 15:00 on August 19th (middle tide), and 10:00 on August 23th to 11:00 on August 24th (spring tide). In addition, the temporary tide level station of Hulu Island also provided the hourly tide level data from August 5th, 2013 to September 4th, 2013. The measuring points’ location is shown in Figure2, in which 0 is Hulu Island tide test point, B1 to B5 are the test points for the flow velocity and flow direction of tidal farm.Tidal model is validated by a long series pattern, and the model is calculated from 3:00 on August 16th, 2013 to 14:00 on August 24th, which lasted 204 hours. The measured data of the tidal level at the temporary test points on the western shore of Hulu Island are compared with the calculated results, as shown in Figure9.Figure 9
Validation of tidal level.By contrast, the calculated results are in good agreement with the measured results, the phase of the calculated results are consistent with the measured results, and the relative error of the tidal height is less than 10%, which indicates that the model is satisfactory for the verification of tide level and can simulate the change of tide level in the sea area more accurately. Taking the middle tide process for example, the validation results of flow velocity and flow direction of each measuring point are shown in Figure10.Figure 10
Validation of flow velocity and direction: (a), (c), (e), (g) and (i) are validation results of flow velocity; (b), (d), (f), (h) and (j) are validation results of flow direction.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)On the whole, the variation trend of flow velocity and flow direction of measuring points is basically the same as measured value. From the point of time corresponding to flow peak and the change of flow direction, the two are also relatively close. Relative error between the simulated value and the calculated value of flow velocity is less than 10%, and the error can be controlled within 18% in addition to some special moments.From the foregoing, the velocity distribution of the tidal farm in the computational area is reasonable, the calculated results of flow velocity and flow direction are consistent with the measured data, and the variation trend is also similar, which shows that the mathematical model of the tidal current is basically reliable and can be used in the hydrodynamic study of the tidal farm.


## 3.1. Establishment of 3-D Near-Field Hydrodynamic Model

### 3.1.1. Turbine Model and Grids Generation

There are two different carrier forms of tidal turbine with horizontal axis, variable pitch, and double impellers. One is for column type with single pile; that is, the structure is supported by a single pile and the impeller can be up and down according to the work needs. The other is for floating; the unit is supported by the floating platform, and the entire device is fixed by four anchor chains. Both the two carrier forms are double impellers with a rotary diameter (D0) of 16.3 m, and the equivalent diameter (DE=2D0) [21] of an impeller is about 23 m.The role of the carrier is ignored in the 3- and 2-D numerical simulation, and only the effect of the double impellers on the flow field is considered. The rotation direction of the impellers is opposite, the spacing between impellers is 3 m, and the hub diameter is 0.1D0. Then, the 3-D modeling software is used to simplify and establish the model of double impellers. The turbine model is shown in Figure 3.Figure 3
3-D turbine model.According to the calculation needs, a cylinder flow field is designed, as shown in Figure4(a). The length of the flow field is 300 m (X direction), the diameter is 100 m (Y and Z direction), the longitudinal axis is parallel to the inflow direction, and the center of the double impellers is located in the longitudinal axis which is away from the upstream inlet 60 m and 240 m from the downstream outlet. This 3-D region is divided by unstructured tetrahedral mesh, in which the green parts represent the rotating domains of impellers and other parts are the nonrotating domains, as shown in Figure 4(b). The grids of rotating domains where the impellers are located are refined, and the grids far away from the impellers are sparse. The maximum grid size on the surface of the impellers (i.e., nearby the walls) is 0.02 m, the maximum grid size of the rotating domains is 0.1 m, and the maximum grid size of the nonrotating domains is 0.6 m. And then, a total of 3,489,648 grids are generated.Figure 4
Computational domain and grids: (a) schematic of computational domain and (b) 3-D grids.
(a) (b)The density of grids in the range of 240 m from rotating domain to nonrotating domain in the cylinder flow field decreases gradually, and the maximum size of the grids increases from 0.1 m to 0.6 m. Moreover, there is no too large or too small size mutation. All these conform to the requirements of renormalization group (RNG) k-ε turbulence model for the calculating grids. In the simulation of Fluent, the calculation equations of RNG k-ε turbulence model are as follows:(1)∂∂tρk+∂∂xiρkui=∂∂xjαkμeff∂k∂xj+Gk+Gb-ρε-YM+Sk(2)∂∂tρε+∂∂xiρεui=∂∂xjαεμeff∂ε∂xj+C1εεkGk+C3εGb-C2ερε2k-Rε+Sεwhere k is turbulence kinetic energy, ε is energy dissipating rate, xi is spatial scale, ui is flow velocity, αk and αε are the Prandtl numbers’ reciprocals of effective turbulent of k and ε, respectively, μeff is viscosity of effective vortex, Gk is turbulence kinetic energy generated by average velocity gradient, Gb is turbulence kinetic energy generated by flotage, YM is influence of turbulent pulsation expansion on overall dissipating rate, and Sk and Sε are user-defined source items.This paper selects 3 points which are away from the impeller surface 50 m, 100 m, and 150 m on the axis of impeller in Condition 1 (inflow velocity is 1.843 m/s), and the three points are recorded as Point 1, 2, and 3, respectively. The changes of velocity of the 3 points under different numbers of grids are shown in Figure5.Figure 5
Variations in velocity with grid numbers.From Figure5, it is known that when the grid number is greater than 3 million, the variations of velocity of the 3 points are very small with the increase of grid number, thus indicating that the grid number has little effect on the simulation results. Therefore, when the grid number is higher than 3 million, the simulation results almost do not change with the variations of the grid number, and results are undoubtedly convergent. So the grid density has little effect on the simulation results; that is, the grid number is almost independent of the simulation results.


### 3.1.2. Relationship between Impeller’s Rotating Speed and Inflow Velocity

Since the turbine model in this paper is constructed by imitating the shape of SeaGen impeller of tidal turbine, its radiusr=8 m, rated flow velocity vr=2.25 m/s, and rated revolving speed of turbine nr=14.3 rad/min. From those, we could obtain that the rated angular speed ωr=2πnr=1.5 rad/s and the design tip speed ratio λ=rωr/vr=5.3. Therefore, the design tip speed ratio of impeller λ is taken as 5.3 to establish the corresponding relationship between the impeller’s rotating speed and the inflow velocity, as shown in Table 1.Table 1
Relationship between impeller’s rotating speed and inflow velocity.
Inflow Velocity (m/s) Parameters Impeller Rotating Position per Second Relationship Schematic between Impeller’s Rotating  Speed/Angular Speed and Inflow Velocity Angular Speed (rad/s) Rotating Speed (r/min) 0.614 0.4 3.82     0.921 0.6 5.73     1.229 0.8 7.64     1.536 1.0 9.55     1.843 1.2 11.46     2.150 1.4 13.37     2.457 1.6 15.28


### 3.1.3. Solver Setting

The CFD software, Fluent, is applied in numerical simulation, and the specific setting of the CFD solver is as follows:( 1 ) Select the 3-D single-precision solver.( 2 ) Select the pressure-based correction method, and select the steady state for the time type.( 3 ) Select RNG k-ε turbulence model and Multiple Reference Frame (MRF) model.( 4 ) The inflow velocity of the velocity-inlet is vi, because the outlet velocity or pressure cannot be known before the problem is solved. Free outflow is used at the outlet boundary, the side wall is a stationary wall, and the impellers are moving walls.( 5 ) 1000 iterations each time step.


## 3.1.1. Turbine Model and Grids Generation

There are two different carrier forms of tidal turbine with horizontal axis, variable pitch, and double impellers. One is for column type with single pile; that is, the structure is supported by a single pile and the impeller can be up and down according to the work needs. The other is for floating; the unit is supported by the floating platform, and the entire device is fixed by four anchor chains. Both the two carrier forms are double impellers with a rotary diameter (D0) of 16.3 m, and the equivalent diameter (DE=2D0) [21] of an impeller is about 23 m.The role of the carrier is ignored in the 3- and 2-D numerical simulation, and only the effect of the double impellers on the flow field is considered. The rotation direction of the impellers is opposite, the spacing between impellers is 3 m, and the hub diameter is 0.1D0. Then, the 3-D modeling software is used to simplify and establish the model of double impellers. The turbine model is shown in Figure 3.Figure 3
3-D turbine model.According to the calculation needs, a cylinder flow field is designed, as shown in Figure4(a). The length of the flow field is 300 m (X direction), the diameter is 100 m (Y and Z direction), the longitudinal axis is parallel to the inflow direction, and the center of the double impellers is located in the longitudinal axis which is away from the upstream inlet 60 m and 240 m from the downstream outlet. This 3-D region is divided by unstructured tetrahedral mesh, in which the green parts represent the rotating domains of impellers and other parts are the nonrotating domains, as shown in Figure 4(b). The grids of rotating domains where the impellers are located are refined, and the grids far away from the impellers are sparse. The maximum grid size on the surface of the impellers (i.e., nearby the walls) is 0.02 m, the maximum grid size of the rotating domains is 0.1 m, and the maximum grid size of the nonrotating domains is 0.6 m. And then, a total of 3,489,648 grids are generated.Figure 4
Computational domain and grids: (a) schematic of computational domain and (b) 3-D grids.
(a) (b)The density of grids in the range of 240 m from rotating domain to nonrotating domain in the cylinder flow field decreases gradually, and the maximum size of the grids increases from 0.1 m to 0.6 m. Moreover, there is no too large or too small size mutation. All these conform to the requirements of renormalization group (RNG) k-ε turbulence model for the calculating grids. In the simulation of Fluent, the calculation equations of RNG k-ε turbulence model are as follows:(1)∂∂tρk+∂∂xiρkui=∂∂xjαkμeff∂k∂xj+Gk+Gb-ρε-YM+Sk(2)∂∂tρε+∂∂xiρεui=∂∂xjαεμeff∂ε∂xj+C1εεkGk+C3εGb-C2ερε2k-Rε+Sεwhere k is turbulence kinetic energy, ε is energy dissipating rate, xi is spatial scale, ui is flow velocity, αk and αε are the Prandtl numbers’ reciprocals of effective turbulent of k and ε, respectively, μeff is viscosity of effective vortex, Gk is turbulence kinetic energy generated by average velocity gradient, Gb is turbulence kinetic energy generated by flotage, YM is influence of turbulent pulsation expansion on overall dissipating rate, and Sk and Sε are user-defined source items.This paper selects 3 points which are away from the impeller surface 50 m, 100 m, and 150 m on the axis of impeller in Condition 1 (inflow velocity is 1.843 m/s), and the three points are recorded as Point 1, 2, and 3, respectively. The changes of velocity of the 3 points under different numbers of grids are shown in Figure5.Figure 5
Variations in velocity with grid numbers.From Figure5, it is known that when the grid number is greater than 3 million, the variations of velocity of the 3 points are very small with the increase of grid number, thus indicating that the grid number has little effect on the simulation results. Therefore, when the grid number is higher than 3 million, the simulation results almost do not change with the variations of the grid number, and results are undoubtedly convergent. So the grid density has little effect on the simulation results; that is, the grid number is almost independent of the simulation results.


## 3.1.2. Relationship between Impeller’s Rotating Speed and Inflow Velocity

Since the turbine model in this paper is constructed by imitating the shape of SeaGen impeller of tidal turbine, its radiusr=8 m, rated flow velocity vr=2.25 m/s, and rated revolving speed of turbine nr=14.3 rad/min. From those, we could obtain that the rated angular speed ωr=2πnr=1.5 rad/s and the design tip speed ratio λ=rωr/vr=5.3. Therefore, the design tip speed ratio of impeller λ is taken as 5.3 to establish the corresponding relationship between the impeller’s rotating speed and the inflow velocity, as shown in Table 1.Table 1
Relationship between impeller’s rotating speed and inflow velocity.
Inflow Velocity (m/s) Parameters Impeller Rotating Position per Second Relationship Schematic between Impeller’s Rotating  Speed/Angular Speed and Inflow Velocity Angular Speed (rad/s) Rotating Speed (r/min) 0.614 0.4 3.82     0.921 0.6 5.73     1.229 0.8 7.64     1.536 1.0 9.55     1.843 1.2 11.46     2.150 1.4 13.37     2.457 1.6 15.28


## 3.1.3. Solver Setting

The CFD software, Fluent, is applied in numerical simulation, and the specific setting of the CFD solver is as follows:( 1 ) Select the 3-D single-precision solver.( 2 ) Select the pressure-based correction method, and select the steady state for the time type.( 3 ) Select RNG k-ε turbulence model and Multiple Reference Frame (MRF) model.( 4 ) The inflow velocity of the velocity-inlet is vi, because the outlet velocity or pressure cannot be known before the problem is solved. Free outflow is used at the outlet boundary, the side wall is a stationary wall, and the impellers are moving walls.( 5 ) 1000 iterations each time step.


## 3.2. Turbine Simplification

For the flow passage where the turbine is located, the flux of any flow section of this flow passage is equal when it is assumed that there is no flow crossing. According to the continuity equation, the acreage of flow sections of the downstream flow passage will be greater than the upstream because the flow velocity of the downstream of the turbine will decrease, as shown in Figure6.Figure 6
Flow passage parameters.The turbine causes partial energy loss of the water flow; it can be expressed by the frictional resistance losshf between flow sections along upstream to downstream of the flow passage, and the corresponding equivalent roughness nE is obtained, then the roughness n can be obtained by combining the original natural roughness n0 of the grid, so as to realize the simplification of the effect of the turbine on the hydrodynamic force, that is, the local roughness correction method.The cross section of the impeller is taken as the reference plane, two flow sections are taken in the upstream L1 and downstream L2, respectively, and the water body surrounded by the two flow sections and the side wall of the flow tube is studied. It is assumed that the average flow velocity of the flow section 1-1 at a certain time isv1, the pressure is p1, and the average flow velocity of the flow section 2-2 is v2, the pressure is p2. The energy equation of the constant pipe flow of incompressible fluid is established, as shown in (3), and only the resistance loss along the path is considered. In (3), α1=α2=1 and z1=z2=0, because the central axis is regarded as the datum plane.(3)z1+p1ρg+α1v122g=z2+p2ρg+α2v222g+hf1-2And then, (3) can be written as (4).(4)hf1-2=p1-p2ρg+v12-v222gAccording to Darcy Formula,hf1-2 can be expressed as (5)hf1-2=λL4Rv122gThe approximate relationship between the resistance coefficient (λ) and the Chézy Coefficient (C) is shown as (6)λ=8gC2The relationship betweenC and nE is obtained by the Manning Formula, as shown in the following.(7)C=1nER1/6The equivalent roughness (nE) of tidal current power generation equipment is obtained by synthesizing (3) to (7), as shown in (8)nE=hf1-2R4/3Lv12=p1-p2ρg+v12-v222gR4/3Lv12The wall shear stress formula is shown as(9)τ=ρgRJ=ρgv2C2=ρgn2v2R1/6The shear stressτ of the seabed is considered to be composed of two parts. One is the equivalent shear stress Δτ produced by turbine’s layout, the corresponding equivalent roughness is nE, and the flow velocity is v2. The other one is the shear stress of the sea itself τ0 when the flow velocity is v2, and the corresponding roughness of seabed is n0. According to the superposition principle, the shear stress of the seabed can be expressed as (10)τ=τ0+Δτ=ρgRJ=ρgn02v22R1/6+ρgnE2v22R1/6=ρgn2v22R1/6The compositive roughness of the turbine after it is placed is as shown in(11)n=n02+nE2


## 3.3. Establishment of 2-D Far-Field Hydrodynamic Model

### 3.3.1. Finite Element Discretization of Governing Equation

The 3-D N-S equation is integrated along the depth direction (Z direction) to obtain the 2-D shallow water circulation equation, which is used as the governing equation. Meanwhile, the role of wind is ignored, as shown in the following:(12)∂z∂t+∂uh∂x+∂vh∂y=0(13)∂u∂t+u∂u∂x+v∂u∂y+g∂z∂x=-gu2+v2C2hu+fv+Mx(14)∂v∂t+u∂v∂x+v∂v∂y+g∂z∂y=-gu2+v2C2hv-fu+Mywhere Mx=AH(∂2u/∂x2+∂2u/∂y2), My=AH(∂2v/∂x2+∂2v/∂y2) are vortex viscosity terms, AH is vortex viscosity coefficient and its value ranges from 0.4 to 0.8, z is tide level, h is water depth, z=h+z0, z0 is bottom elevation, t is time, u and v are components of velocity in x and y direction, respectively, C is Chézy coefficient, f is Coriolis force coefficient, and g is gravitational acceleration.In this paper, FEM is used to discretize the 2-D shallow water circulation equation. The computational area is divided by the uniform triangle grids, the number of unstructured grids isM, and the number of nodes is N. In addition, φn (n=1,2,…,N) is set as the basis function of the solution function of the finite element equation, and φpe is set as the shape function of local units. Then the solution of the governing equation can be expressed as (15)ut,x,y=∑n=1Nuntφnx,y=∑e=1M∑r=13φpreurevt,x,y=∑n=1Nvntφnx,y=∑e=1M∑r=13φprevrezt,x,y=∑n=1Nzntφnx,y=∑e=1M∑r=13φprezrewhere the part after the first equal sign is the solution function under the condition that the grid nodes are arranged by ordinal number in the whole coordinate, the part after the second equal sign is the solution function under the condition that the grids are arranged by ordinal number in the local coordinate. The solution of arbitrary units in (15) is given by (16)uet,x,y=∑r=13φprex,yuretvet,x,y=∑r=13φprex,yvretzet,x,y=∑r=13φprex,yzretEquation (16) is substituted into (12) to (14), he(x,y,t)=ze(x,y,t)-z0(x,y), then continuity Equation (17) and motion Equation (18) of triangle grids are obtained by Galerkin method.(17)∫Ω∑r=13φpre∂zre∂t+∑j=13∂φpje∂xuje∑r=13φprehre+∑j=13∂φpje∂xhje∑r=13φpreure+∑j=13∂φpje∂yvje∑r=13φprehre+∑j=13∂φpje∂yhje∑r=13φprevreφqkedΩ=0k=1,2,3(18)∫Ω∑r=13φpre∂ure∂t+∑j=13∂φpje∂xuje∑r=13φpreure+∑j=13∂φpje∂yuje∑r=13φprevre+g∑r=13∂φpre∂xzre+g∑j=13φpjeuje2+∑j=13φpjevje2∑j=13cj2φpjehje∑r=13φpreure-∑r=13frφprevreφqkedΩΩ=0k=1,2,3A r k = ∫ Ω φ p ( r ) e φ q ( k ) e d Ω, Brk=∑j=13(uje∂φpje/∂x+vje∂φpje/∂y)Ark, Crk=∑j=13(hje∂φpje/∂x)Ark, Drk=∑j=13(hje∂φpje/∂y)Ark are substituted into Equation (17), and Erk=∑j=13(uje∂φpje/∂x)Ark, Erk′=∑j=13(vje∂φpje/∂x)Ark, Frk=∑j=13(uje∂φpje/∂y)Ark, Frk′=∑j=13(vje∂φpje/∂y)Ark, Hrk=g∫Ω(∂φpre/∂x)φq(k)edΩ, Hrk′=g∫Ω∂φpre/∂yφq(k)edΩ, Jrk=∫Ωfrφp(r)eφq(k)edΩ, Irk = g∫Ωφpreφqke((∑j=13φpjeuje)2+(∑j=13φpjevje)2/ ∑j=13cj2φpjehje)dΩ are substituted into Equation (18). Then, finite element equation of element matrix Equation (19) is obtained by changing the order of integration summation.(19)Az˙+Bh+Cu+Dv=0Au˙+Eu+Fv+Hz+Iu-Jv=0Av˙+E′u+F′v+H′z+Iv+Ju=0where z˙=∂z1/∂t∂z2/∂t∂z3/∂te, u˙=∂u1/∂t∂u2/∂t∂u3/∂te, v˙=∂v1/∂t∂v2/∂t∂v3/∂te, z=z1z2z3e, h=h1h2h3e, u=u1u2u3e, v=v1v2v3e, A=A11A12A13A21A22A23A31A32A33, B, C, D, E, F, H, I, J, E′, F′, H′ are similar to A, and they are 3×3 matrixes which correspond to the elements corresponding to the unit.After the unit partition is determined, the number of the local unit nodes in whole coordinate is determined to form the whole matrix (increase the subscriptT to represent). Moreover, the time differential is discretized by the forward difference, and the moment to be sought is expressed by the superscript t+Δt, and the remainder is t. Then (20) is obtained.(20)ATzt+Δt=ATzt-ΔtBTht+CTut+DTvtATut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtATvt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtEquation (20) can be rewritten as (21) by using mass concentration method, diagonalizing 3×3 matrix starting from unit matrix, replacing matrix AT at time t+Δt by diagonal matrix A0, and then solving ut+Δt and vt+Δt of the equation and finally ut+Δt and vt+Δt are substituted into the equation zt+Δt.(21)A0ut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtA0vt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtA0zt+Δt=ATzt-ΔtBTht+CTut+Δt+DTvt+Δt


### 3.3.2. Grids Generation and Submarine Topography

The study area mainly includes the northeastern coast of Putuo Island and Hulu Island, the quadrilateral open sea facing the East China Sea. The study area is about 10 km from west to east, about 8 km from south to north, covering about 74 km2. The unstructured grids are used to divide the computational area in the model, and the density of the grids increases gradually from the open sea area to the island coastline. Finally, the number of unstructured grids’ nodes is 7,970 and the number of unstructured grids is 15,264. The grids generation of tidal model is shown in Figure 7. According to the data, the contour map of submarine elevation is obtained through the transformation of coordinate and datum plane, the interpolation of grid nodes, and a series of processing, as shown in Figure 8.Figure 7
Grids generation of tidal model.Figure 8
Submarine elevation contour.Comparing Figures7 and 8, it could be found that the topography near the open boundary is relatively flat and the simulation of tidal wave is more stable. Therefore, increasing the grid size and reducing the resolution is reasonable but also can improve the model calculation efficiency. Increasing the density of the grid at the waterway and the island coast can improve the accuracy of the operation and fully reflect the submarine topography.The type of the open boundary of the model is the tide boundary; that is, the water level process of each grid node of the determined boundary is given along the open boundary. The water level process is obtained by harmonic analysis method and it is predicted by the global ocean tide model [22], and eight tidal constituents M2, S2, N2, K2, K1, O1, P1, Q1 are taken into account.


### 3.3.3. Validation of 2-D Model

The marine hydrological data is obtained by the Second Institute of Oceanography of the State Oceanic Administration in August 2013, such as tide and tidal wave. The observation period was 10:00 on August 16th, 2013 to 17:00 on August 17th, 2013 (small tide), 14:00 on August 19th to 15:00 on August 19th (middle tide), and 10:00 on August 23th to 11:00 on August 24th (spring tide). In addition, the temporary tide level station of Hulu Island also provided the hourly tide level data from August 5th, 2013 to September 4th, 2013. The measuring points’ location is shown in Figure2, in which 0 is Hulu Island tide test point, B1 to B5 are the test points for the flow velocity and flow direction of tidal farm.Tidal model is validated by a long series pattern, and the model is calculated from 3:00 on August 16th, 2013 to 14:00 on August 24th, which lasted 204 hours. The measured data of the tidal level at the temporary test points on the western shore of Hulu Island are compared with the calculated results, as shown in Figure9.Figure 9
Validation of tidal level.By contrast, the calculated results are in good agreement with the measured results, the phase of the calculated results are consistent with the measured results, and the relative error of the tidal height is less than 10%, which indicates that the model is satisfactory for the verification of tide level and can simulate the change of tide level in the sea area more accurately. Taking the middle tide process for example, the validation results of flow velocity and flow direction of each measuring point are shown in Figure10.Figure 10
Validation of flow velocity and direction: (a), (c), (e), (g) and (i) are validation results of flow velocity; (b), (d), (f), (h) and (j) are validation results of flow direction.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)On the whole, the variation trend of flow velocity and flow direction of measuring points is basically the same as measured value. From the point of time corresponding to flow peak and the change of flow direction, the two are also relatively close. Relative error between the simulated value and the calculated value of flow velocity is less than 10%, and the error can be controlled within 18% in addition to some special moments.From the foregoing, the velocity distribution of the tidal farm in the computational area is reasonable, the calculated results of flow velocity and flow direction are consistent with the measured data, and the variation trend is also similar, which shows that the mathematical model of the tidal current is basically reliable and can be used in the hydrodynamic study of the tidal farm.


## 3.3.1. Finite Element Discretization of Governing Equation

The 3-D N-S equation is integrated along the depth direction (Z direction) to obtain the 2-D shallow water circulation equation, which is used as the governing equation. Meanwhile, the role of wind is ignored, as shown in the following:(12)∂z∂t+∂uh∂x+∂vh∂y=0(13)∂u∂t+u∂u∂x+v∂u∂y+g∂z∂x=-gu2+v2C2hu+fv+Mx(14)∂v∂t+u∂v∂x+v∂v∂y+g∂z∂y=-gu2+v2C2hv-fu+Mywhere Mx=AH(∂2u/∂x2+∂2u/∂y2), My=AH(∂2v/∂x2+∂2v/∂y2) are vortex viscosity terms, AH is vortex viscosity coefficient and its value ranges from 0.4 to 0.8, z is tide level, h is water depth, z=h+z0, z0 is bottom elevation, t is time, u and v are components of velocity in x and y direction, respectively, C is Chézy coefficient, f is Coriolis force coefficient, and g is gravitational acceleration.In this paper, FEM is used to discretize the 2-D shallow water circulation equation. The computational area is divided by the uniform triangle grids, the number of unstructured grids isM, and the number of nodes is N. In addition, φn (n=1,2,…,N) is set as the basis function of the solution function of the finite element equation, and φpe is set as the shape function of local units. Then the solution of the governing equation can be expressed as (15)ut,x,y=∑n=1Nuntφnx,y=∑e=1M∑r=13φpreurevt,x,y=∑n=1Nvntφnx,y=∑e=1M∑r=13φprevrezt,x,y=∑n=1Nzntφnx,y=∑e=1M∑r=13φprezrewhere the part after the first equal sign is the solution function under the condition that the grid nodes are arranged by ordinal number in the whole coordinate, the part after the second equal sign is the solution function under the condition that the grids are arranged by ordinal number in the local coordinate. The solution of arbitrary units in (15) is given by (16)uet,x,y=∑r=13φprex,yuretvet,x,y=∑r=13φprex,yvretzet,x,y=∑r=13φprex,yzretEquation (16) is substituted into (12) to (14), he(x,y,t)=ze(x,y,t)-z0(x,y), then continuity Equation (17) and motion Equation (18) of triangle grids are obtained by Galerkin method.(17)∫Ω∑r=13φpre∂zre∂t+∑j=13∂φpje∂xuje∑r=13φprehre+∑j=13∂φpje∂xhje∑r=13φpreure+∑j=13∂φpje∂yvje∑r=13φprehre+∑j=13∂φpje∂yhje∑r=13φprevreφqkedΩ=0k=1,2,3(18)∫Ω∑r=13φpre∂ure∂t+∑j=13∂φpje∂xuje∑r=13φpreure+∑j=13∂φpje∂yuje∑r=13φprevre+g∑r=13∂φpre∂xzre+g∑j=13φpjeuje2+∑j=13φpjevje2∑j=13cj2φpjehje∑r=13φpreure-∑r=13frφprevreφqkedΩΩ=0k=1,2,3A r k = ∫ Ω φ p ( r ) e φ q ( k ) e d Ω, Brk=∑j=13(uje∂φpje/∂x+vje∂φpje/∂y)Ark, Crk=∑j=13(hje∂φpje/∂x)Ark, Drk=∑j=13(hje∂φpje/∂y)Ark are substituted into Equation (17), and Erk=∑j=13(uje∂φpje/∂x)Ark, Erk′=∑j=13(vje∂φpje/∂x)Ark, Frk=∑j=13(uje∂φpje/∂y)Ark, Frk′=∑j=13(vje∂φpje/∂y)Ark, Hrk=g∫Ω(∂φpre/∂x)φq(k)edΩ, Hrk′=g∫Ω∂φpre/∂yφq(k)edΩ, Jrk=∫Ωfrφp(r)eφq(k)edΩ, Irk = g∫Ωφpreφqke((∑j=13φpjeuje)2+(∑j=13φpjevje)2/ ∑j=13cj2φpjehje)dΩ are substituted into Equation (18). Then, finite element equation of element matrix Equation (19) is obtained by changing the order of integration summation.(19)Az˙+Bh+Cu+Dv=0Au˙+Eu+Fv+Hz+Iu-Jv=0Av˙+E′u+F′v+H′z+Iv+Ju=0where z˙=∂z1/∂t∂z2/∂t∂z3/∂te, u˙=∂u1/∂t∂u2/∂t∂u3/∂te, v˙=∂v1/∂t∂v2/∂t∂v3/∂te, z=z1z2z3e, h=h1h2h3e, u=u1u2u3e, v=v1v2v3e, A=A11A12A13A21A22A23A31A32A33, B, C, D, E, F, H, I, J, E′, F′, H′ are similar to A, and they are 3×3 matrixes which correspond to the elements corresponding to the unit.After the unit partition is determined, the number of the local unit nodes in whole coordinate is determined to form the whole matrix (increase the subscriptT to represent). Moreover, the time differential is discretized by the forward difference, and the moment to be sought is expressed by the superscript t+Δt, and the remainder is t. Then (20) is obtained.(20)ATzt+Δt=ATzt-ΔtBTht+CTut+DTvtATut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtATvt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtEquation (20) can be rewritten as (21) by using mass concentration method, diagonalizing 3×3 matrix starting from unit matrix, replacing matrix AT at time t+Δt by diagonal matrix A0, and then solving ut+Δt and vt+Δt of the equation and finally ut+Δt and vt+Δt are substituted into the equation zt+Δt.(21)A0ut+Δt=ATut-ΔtHTzt+ET+ITut+FT-JTvtA0vt+Δt=ATvt-ΔtH′Tzt+E′T+JTut+F′T+ITvtA0zt+Δt=ATzt-ΔtBTht+CTut+Δt+DTvt+Δt


## 3.3.2. Grids Generation and Submarine Topography

The study area mainly includes the northeastern coast of Putuo Island and Hulu Island, the quadrilateral open sea facing the East China Sea. The study area is about 10 km from west to east, about 8 km from south to north, covering about 74 km2. The unstructured grids are used to divide the computational area in the model, and the density of the grids increases gradually from the open sea area to the island coastline. Finally, the number of unstructured grids’ nodes is 7,970 and the number of unstructured grids is 15,264. The grids generation of tidal model is shown in Figure 7. According to the data, the contour map of submarine elevation is obtained through the transformation of coordinate and datum plane, the interpolation of grid nodes, and a series of processing, as shown in Figure 8.Figure 7
Grids generation of tidal model.Figure 8
Submarine elevation contour.Comparing Figures7 and 8, it could be found that the topography near the open boundary is relatively flat and the simulation of tidal wave is more stable. Therefore, increasing the grid size and reducing the resolution is reasonable but also can improve the model calculation efficiency. Increasing the density of the grid at the waterway and the island coast can improve the accuracy of the operation and fully reflect the submarine topography.The type of the open boundary of the model is the tide boundary; that is, the water level process of each grid node of the determined boundary is given along the open boundary. The water level process is obtained by harmonic analysis method and it is predicted by the global ocean tide model [22], and eight tidal constituents M2, S2, N2, K2, K1, O1, P1, Q1 are taken into account.


## 3.3.3. Validation of 2-D Model

The marine hydrological data is obtained by the Second Institute of Oceanography of the State Oceanic Administration in August 2013, such as tide and tidal wave. The observation period was 10:00 on August 16th, 2013 to 17:00 on August 17th, 2013 (small tide), 14:00 on August 19th to 15:00 on August 19th (middle tide), and 10:00 on August 23th to 11:00 on August 24th (spring tide). In addition, the temporary tide level station of Hulu Island also provided the hourly tide level data from August 5th, 2013 to September 4th, 2013. The measuring points’ location is shown in Figure2, in which 0 is Hulu Island tide test point, B1 to B5 are the test points for the flow velocity and flow direction of tidal farm.Tidal model is validated by a long series pattern, and the model is calculated from 3:00 on August 16th, 2013 to 14:00 on August 24th, which lasted 204 hours. The measured data of the tidal level at the temporary test points on the western shore of Hulu Island are compared with the calculated results, as shown in Figure9.Figure 9
Validation of tidal level.By contrast, the calculated results are in good agreement with the measured results, the phase of the calculated results are consistent with the measured results, and the relative error of the tidal height is less than 10%, which indicates that the model is satisfactory for the verification of tide level and can simulate the change of tide level in the sea area more accurately. Taking the middle tide process for example, the validation results of flow velocity and flow direction of each measuring point are shown in Figure10.Figure 10
Validation of flow velocity and direction: (a), (c), (e), (g) and (i) are validation results of flow velocity; (b), (d), (f), (h) and (j) are validation results of flow direction.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)On the whole, the variation trend of flow velocity and flow direction of measuring points is basically the same as measured value. From the point of time corresponding to flow peak and the change of flow direction, the two are also relatively close. Relative error between the simulated value and the calculated value of flow velocity is less than 10%, and the error can be controlled within 18% in addition to some special moments.From the foregoing, the velocity distribution of the tidal farm in the computational area is reasonable, the calculated results of flow velocity and flow direction are consistent with the measured data, and the variation trend is also similar, which shows that the mathematical model of the tidal current is basically reliable and can be used in the hydrodynamic study of the tidal farm.


## 4. Results and Discussions

### 4.1. 3-D Near-Field Hydrodynamic Characteristics

Based on CFD results, velocity distributions in center plane of double impellers in each combination are shown in Table2. Because blade’s diameter is larger, linear velocity near blade is much greater than velocity of inflow and other flow field. In order to reflect effect of impellers on near-field flow and compare similarities and differences between flow fields better, area with large flow velocity near blade will be ignored, and area with velocity distributions from 0 to 1.14vi is emphasized when drawing contours.Table 2
Velocity distributions on impellers’ center planes in different combinations.
v i - ω i Combinations Flow Velocity Contours on Impellers’ Center Planes ① Inflow Velocityv1 0.614 (m/s)  Angular Speedω1 0.4 (rad/s)  ② Inflow Velocityv2 0.921 (m/s)  Angular Speedω2 0.6 (rad/s)  ③ Inflow Velocityv3 1.229 (m/s)  Angular Speedω3 0.8 (rad/s)  ④ Inflow Velocityv4 1.536 (m/s)  Angular Speedω4 1.0 (rad/s)  ⑤ Inflow Velocityv5 1.843 (m/s)  Angular Speedω5 1.2 (rad/s)  ⑥ Inflow Velocityv6 2.150 (m/s)  Angular Speedω6 1.4 (rad/s)  ⑦ Inflow Velocityv7 2.457 (m/s)  Angular Speedω7 1.6 (rad/s)Flow velocity profiles of eachvi-ωi combination have some common features, which are manifested in the upstream, downstream, and inner and outer sides of impellers, respectively. As for the upstream of impellers, the upstream inflow is blocked because of the installation and operation of turbine, which is shown as a flow velocity deceleration zone which is convex to inflow in the upstream of impeller in flow velocity profiles. After the flow running through the impellers, the flow velocity decreased significantly because of energy loss, and a strip of wake flow with the length of about 240 m (about 10DE) is formed downstream of each impeller. The wake of impellers gradually closes together in the downstream and merges into one. The width of the wake narrows along the forward direction of X-axis and disappears gradually in the distance. The gradient distribution of the change of flow velocity in the strip of wake flow is not well-distributed, and the closer the velocity of the impeller is, the more severe the loss is. On the outside at downstream of the impellers, due to turbine’s layout and the widening of the strip of wake flow, the water area on both sides decreases and the flow velocity increases slightly (relative to the inflow velocity). Similarly, there will be a “separation zone” on the inside at downstream of the impellers, where the flow velocity is larger than the wake flow.With differentvi-ωi combinations, the flow velocity distributions are different. The greater the inflow velocity is, the faster the impellers rotate, and the greater the hindrance function of the impellers on the upstream is. As for the length of the strip of wake flow, the green contours are taken as an example, the larger the impellers’ angular speed, the longer the length of the green strip of wake flow, and then the position where the strip of wake flow narrows obviously moves downstream, which results in the length of the area where the flow velocity increased (dark red) in the outside of the double impellers is longer in the X direction, and the space on both sides of the wake is gradually filled in the Y direction.The flow velocity of the wake at the impeller axis in these 7 combinations has been restored to 96.9%, 94.1%, 94.0%, 91.1%, 91.2%, 88.3%, and 88.3% of their inflow velocity, respectively. In this paper, we take the investigation of Yang as a criterion to determine whether the computational range is normal enough; that is, when the flow velocity is restored to 80% ~ 90% of the inlet velocity, it is determined that it meets the energy generation efficiency and the work needs [23]. The velocity of the wake can be restored to more than 80% of the inflow velocity, so the computational range of the 3-D model is sufficient.


### 4.2. Turbine Simplification Results

Equivalent roughness is calculated according to (8), then compositive roughness is obtained from (11) by combining with natural roughness of tidal farm which is 0.025, the results shown in Table 3.Table 3
Simulated results in different combinations.
Combinations Pressurep1 (Pa) Pressurep2 (Pa) Energy Losshf Resistance Coefficientλ Chézy Coefficient C Equivalent RoughnessnE Compositive Roughnessn No. v ω (m/s) (rad/s) 1 0.614 0.4 101097 101079 0.015 1.598 7.005 0.274 0.275 2 0.921 0.6 100777 100737 0.036 1.687 6.817 0.282 0.283 3 1.229 0.8 100301 100228 0.068 1.760 6.674 0.288 0.289 4 1.536 1.0 99664 99551 0.109 1.809 6.583 0.292 0.293 5 1.843 1.2 98854 98693 0.160 1.848 6.514 0.295 0.296 6 2.15 1.4 97836 97616 0.222 1.887 6.445 0.298 0.299 7 2.457 1.6 96593 96304 0.296 1.923 6.385 0.301 0.302Figure11 is drawn according to the correspondence of each combination, in which the blue points represent the calculated results of compositive roughness, and the black line represents the trend line which is obtained by regression of scatter points. The function is expressed as (22)n=0.01899·ln⁡v+0.28441Figure 11
Relationship between inflow velocity and compositive roughness.Compositive roughness of turbine increases logarithmically with the increasing of inflow velocity.Equation (22) is applied to that the impeller’s rotating speed is changed according to the design tip speed ratio λr= 5.3 with the flow velocity changing, and the equivalent roughness is calculated by the method mentioned above. Equation (22) can be used for the calculation of 2-D tidal current field but should be used according to the actual situation. Based on the “roughness increment - compositive roughness” relationship expressed by (11) and the “compositive roughness - inflow velocity” relationship expressed by (22), the “bridge” which connects 3-D model and 2-D model is constructed. Thus the local roughness correction, that is, the “roughness increment - compositive roughness” relationship is obtained through the 3-D model to simulate the resistance loss (i.e., CFD simulation results of hf) from upstream to downstream, and the resistance loss is converted into equivalent roughness of the turbine by (8). Through regression of the simulated results, the “compositive roughness - inflow velocity” relationship is obtained and substituted into the grids where turbines are located for calculation. The roughness is changed with the velocity in the 2-D flow field, which reflects the effect of the turbine on the flow field; thus, the connection of 3-D model and 2-D model can be achieved.


### 4.3. 2-D Far-Field Hydrodynamic Characteristics

In this paper, 5 different installation locations of the tidal turbine are considered, and the comparison of the flow field before and after the installation is simulated in the process of spring tide under the 5 conditions. Taking the center of the double impeller as the center of the circle, the radius 100 m, 200 m and 300 m is taken as the exploratory circumference to count and compare the changes of flow field. Take the B4 berthage as an example, and the exploratory circumference as shown in Figure12.Figure 12
Exploratory circumference schematic.Difference of flow velocity is the velocity after installation minus the velocity before installation at the same time, and maximum difference of flow velocity represents maximum value of deceleration in flow velocity at a certain time. Distribution of maximum difference of flow velocity near turbine was shown in Figure13. To reflect deceleration effect of device better, Figure 13 only makes the contour lines of maximum difference of flow velocity in range from -0.015 m/s to -0.05 m/s. Since maximum difference of flow velocity does not appear to be the same at all times, area enclosed by contour lines is an extreme case, which means that maximal affecting range of the turbine on flow field at any time will not exceed area shown in Figure 13.Figure 13
Distribution of maximum difference of flow velocity in (a) B1 location; (b) B2 location; (c) B3 location; (d) B4 location; and (e) B5 location.
(a) (b) (c) (d) (e)It can be seen from Figure13 that turbine location at both end and middle of channel has a certain effect on surrounding flow field, and a certain range of velocity deceleration zones which are “8-shaped” areas appears in facing flow side and reverse flow side of turbine, and flow velocity on two sides of turbine increases slightly. In the process of spring tide, difference of flow velocity in near flow field has largest change, and affecting area is biggest when B4 berthage is installed with turbine, B3 berthage is next to B4 berthage, B5 berthage is followed by B3 berthage, B1 berthage is second to B5 berthage, and B2 berthage is smallest. B4 and B3 berthage are in narrow section of channel, and B2 berthage is located in open section of channel. This indicates that the effect of installation of turbine on flow field in narrow sea area is relatively larger, and the effect of installation of turbine on flow field in open sea area is relatively smaller. In addition, water depth of B2 berthage is largest in the 5 berthages, but the effect on the flow field is the smallest. Affecting range of B4 and B3 berthage in channel on flow field is substantially consistent with waterway direction, indicating that the effect of installation of turbine on flow field is also related to water depth, topography, and channel direction.Hourly changing process of extreme value of deceleration and relative difference of flow velocity in different exploratory circumference with radius of 100 m, 200 m, and 300 m are respectively shown in Figures14, 15, and 16. Wherein maximum of deceleration represents maximum value of decrease of flow velocity at every point on circumference at a certain time, definition of minimum of deceleration is just the opposite. The relative difference of flow velocity represents the ratio of the maximum of deceleration to the maximum of flow velocity before the installation of the turbine.Figure 14
Variations in flow velocity of circumference with radius of 100 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)Figure 15
Variations in flow velocity of circumference with radius of 200 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)Figure 16
Variations in flow velocity of circumference with radius of 300 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)As can be seen from these figures, the extremal curves of deceleration at each berthage show a certain periodicity, which is the result of the periodicity of the tidal current. The curves of B2 berthage are special, and the curve patterns of other berthages are similar. Most of the deceleration minimum curves of all circumferences are located above the horizontal axis and most of the deceleration maximum curves are below the horizontal axis. This shows that there is a certain degree of deceleration and growth of flow velocity appearing near the turbine, which is also consistent with the results of the 3-D near-field simulation. The peak value of the deceleration maximum curves of each turbine generally decreases with the increase of the circumferential radius, which means that the farther away from the turbine, the smaller the effect of the turbine on the flow field. From the curves of relative difference of flow velocity, the effect of the turbine on the far-field of flow field is smaller as a whole, and the maximum change in flow velocity does not exceed 8% of the flow velocity before the installation of the turbine. The degree of deceleration and growth of all berthages were compared as shown in Table4.Table 4
Comparison of deceleration and growth degree.
Circumferential Radius Comparison of Deceleration Degree Comparison of Growth Degree 100 m B4> B3> B1> B5> B2 B3> B2> B1> B5> B4 200 m B4> B3> B1> B5> B2 B3> B2> B1> B4> B5 300 m B4> B3> B1> B5> B2 B4> B2> B3> B1> B5As can be seen from Table4, the results of comparison of the deceleration degree (B4> B3> B1> B5> B2) are approximately the same as the comparison of the maximum value distribution of the flow velocity difference (B4> B3> B5> B1> B2). For the growth degree, B4 berthage changes greatly. There is no growth in its 100 m range, but the degree of growth near 300 m is the largest. In order to better reflect the specific changes in the affected area, the change of the extremum of flow velocity difference is formed by different location under the spring tide, as shown in Table 5.Table 5
Extremum of flow velocity variation.
Location Radius 100 m Radius 200 m Radius 300 m Deceleration Growth Relative Difference Deceleration Growth Relative Difference Deceleration Growth Relative Difference (m/s) (m/s) (%) (m/s) (m/s) (%) (m/s) (m/s) (%) B1 -0.023 0.006 1.758 -0.022 0.007 1.754 -0.016 0.006 1.391 B2 -0.001 0.017 0.419 -0.001 0.011 0.226 -0.003 0.009 0.281 B3 -0.064 0.034 4.754 -0.043 0.024 3.354 -0.034 0.007 3.147 B4 -0.091 0.002 7.461 -0.088 0.01 6.877 -0.053 0.024 3.947 B5 -0.019 0.002 1.503 -0.012 0.003 0.953 -0.009 0.003 0.714It can be seen from Table5 that B4 location has the largest effect on the flow field and the affected area is also the most extensive, and B2 location has the least effect on the flow field. The deceleration maximum of B1, B3, B4, and B5 decreases with the increase of radius, the change of B2 location is smaller, and there is fluctuation around the radius of 300 m. The growth maximum of B2 and B3 decreases with the increase of radius, and the growth changes of B1, B4, and B5 are smaller and the fluctuation also appears. This shows that when the degree of deceleration and growth of flow velocity is greater, the value of deceleration and growth will decrease with the increase of radius. When the deceleration and growth of velocity are less than 0.01 m/s, the value of deceleration and growth will only fluctuate and there is no obvious law.To sum up, for the installation of tidal turbine in tidal farm, the following conclusions can be obtained by combining the consideration of the bathymetric data and the tidal resource in the engineering sea area:( 1 ) Compared to other berthages, the effect of test point B4 on the flow field is the largest after turbine’s layout; this point can be considered to avoid in the layout of berthages.( 2 ) Tidal current power resources of B1, B3, and B5 test points are more abundant; the degree of effect after turbine’s layout is also within the acceptable range.( 3 ) Combining the offshore distance of each point, test points B1 and B3 are the ideal layout points for the five berthages through synthetical consideration.


## 4.1. 3-D Near-Field Hydrodynamic Characteristics

Based on CFD results, velocity distributions in center plane of double impellers in each combination are shown in Table2. Because blade’s diameter is larger, linear velocity near blade is much greater than velocity of inflow and other flow field. In order to reflect effect of impellers on near-field flow and compare similarities and differences between flow fields better, area with large flow velocity near blade will be ignored, and area with velocity distributions from 0 to 1.14vi is emphasized when drawing contours.Table 2
Velocity distributions on impellers’ center planes in different combinations.
v i - ω i Combinations Flow Velocity Contours on Impellers’ Center Planes ① Inflow Velocityv1 0.614 (m/s)  Angular Speedω1 0.4 (rad/s)  ② Inflow Velocityv2 0.921 (m/s)  Angular Speedω2 0.6 (rad/s)  ③ Inflow Velocityv3 1.229 (m/s)  Angular Speedω3 0.8 (rad/s)  ④ Inflow Velocityv4 1.536 (m/s)  Angular Speedω4 1.0 (rad/s)  ⑤ Inflow Velocityv5 1.843 (m/s)  Angular Speedω5 1.2 (rad/s)  ⑥ Inflow Velocityv6 2.150 (m/s)  Angular Speedω6 1.4 (rad/s)  ⑦ Inflow Velocityv7 2.457 (m/s)  Angular Speedω7 1.6 (rad/s)Flow velocity profiles of eachvi-ωi combination have some common features, which are manifested in the upstream, downstream, and inner and outer sides of impellers, respectively. As for the upstream of impellers, the upstream inflow is blocked because of the installation and operation of turbine, which is shown as a flow velocity deceleration zone which is convex to inflow in the upstream of impeller in flow velocity profiles. After the flow running through the impellers, the flow velocity decreased significantly because of energy loss, and a strip of wake flow with the length of about 240 m (about 10DE) is formed downstream of each impeller. The wake of impellers gradually closes together in the downstream and merges into one. The width of the wake narrows along the forward direction of X-axis and disappears gradually in the distance. The gradient distribution of the change of flow velocity in the strip of wake flow is not well-distributed, and the closer the velocity of the impeller is, the more severe the loss is. On the outside at downstream of the impellers, due to turbine’s layout and the widening of the strip of wake flow, the water area on both sides decreases and the flow velocity increases slightly (relative to the inflow velocity). Similarly, there will be a “separation zone” on the inside at downstream of the impellers, where the flow velocity is larger than the wake flow.With differentvi-ωi combinations, the flow velocity distributions are different. The greater the inflow velocity is, the faster the impellers rotate, and the greater the hindrance function of the impellers on the upstream is. As for the length of the strip of wake flow, the green contours are taken as an example, the larger the impellers’ angular speed, the longer the length of the green strip of wake flow, and then the position where the strip of wake flow narrows obviously moves downstream, which results in the length of the area where the flow velocity increased (dark red) in the outside of the double impellers is longer in the X direction, and the space on both sides of the wake is gradually filled in the Y direction.The flow velocity of the wake at the impeller axis in these 7 combinations has been restored to 96.9%, 94.1%, 94.0%, 91.1%, 91.2%, 88.3%, and 88.3% of their inflow velocity, respectively. In this paper, we take the investigation of Yang as a criterion to determine whether the computational range is normal enough; that is, when the flow velocity is restored to 80% ~ 90% of the inlet velocity, it is determined that it meets the energy generation efficiency and the work needs [23]. The velocity of the wake can be restored to more than 80% of the inflow velocity, so the computational range of the 3-D model is sufficient.


## 4.2. Turbine Simplification Results

Equivalent roughness is calculated according to (8), then compositive roughness is obtained from (11) by combining with natural roughness of tidal farm which is 0.025, the results shown in Table 3.Table 3
Simulated results in different combinations.
Combinations Pressurep1 (Pa) Pressurep2 (Pa) Energy Losshf Resistance Coefficientλ Chézy Coefficient C Equivalent RoughnessnE Compositive Roughnessn No. v ω (m/s) (rad/s) 1 0.614 0.4 101097 101079 0.015 1.598 7.005 0.274 0.275 2 0.921 0.6 100777 100737 0.036 1.687 6.817 0.282 0.283 3 1.229 0.8 100301 100228 0.068 1.760 6.674 0.288 0.289 4 1.536 1.0 99664 99551 0.109 1.809 6.583 0.292 0.293 5 1.843 1.2 98854 98693 0.160 1.848 6.514 0.295 0.296 6 2.15 1.4 97836 97616 0.222 1.887 6.445 0.298 0.299 7 2.457 1.6 96593 96304 0.296 1.923 6.385 0.301 0.302Figure11 is drawn according to the correspondence of each combination, in which the blue points represent the calculated results of compositive roughness, and the black line represents the trend line which is obtained by regression of scatter points. The function is expressed as (22)n=0.01899·ln⁡v+0.28441Figure 11
Relationship between inflow velocity and compositive roughness.Compositive roughness of turbine increases logarithmically with the increasing of inflow velocity.Equation (22) is applied to that the impeller’s rotating speed is changed according to the design tip speed ratio λr= 5.3 with the flow velocity changing, and the equivalent roughness is calculated by the method mentioned above. Equation (22) can be used for the calculation of 2-D tidal current field but should be used according to the actual situation. Based on the “roughness increment - compositive roughness” relationship expressed by (11) and the “compositive roughness - inflow velocity” relationship expressed by (22), the “bridge” which connects 3-D model and 2-D model is constructed. Thus the local roughness correction, that is, the “roughness increment - compositive roughness” relationship is obtained through the 3-D model to simulate the resistance loss (i.e., CFD simulation results of hf) from upstream to downstream, and the resistance loss is converted into equivalent roughness of the turbine by (8). Through regression of the simulated results, the “compositive roughness - inflow velocity” relationship is obtained and substituted into the grids where turbines are located for calculation. The roughness is changed with the velocity in the 2-D flow field, which reflects the effect of the turbine on the flow field; thus, the connection of 3-D model and 2-D model can be achieved.


## 4.3. 2-D Far-Field Hydrodynamic Characteristics

In this paper, 5 different installation locations of the tidal turbine are considered, and the comparison of the flow field before and after the installation is simulated in the process of spring tide under the 5 conditions. Taking the center of the double impeller as the center of the circle, the radius 100 m, 200 m and 300 m is taken as the exploratory circumference to count and compare the changes of flow field. Take the B4 berthage as an example, and the exploratory circumference as shown in Figure12.Figure 12
Exploratory circumference schematic.Difference of flow velocity is the velocity after installation minus the velocity before installation at the same time, and maximum difference of flow velocity represents maximum value of deceleration in flow velocity at a certain time. Distribution of maximum difference of flow velocity near turbine was shown in Figure13. To reflect deceleration effect of device better, Figure 13 only makes the contour lines of maximum difference of flow velocity in range from -0.015 m/s to -0.05 m/s. Since maximum difference of flow velocity does not appear to be the same at all times, area enclosed by contour lines is an extreme case, which means that maximal affecting range of the turbine on flow field at any time will not exceed area shown in Figure 13.Figure 13
Distribution of maximum difference of flow velocity in (a) B1 location; (b) B2 location; (c) B3 location; (d) B4 location; and (e) B5 location.
(a) (b) (c) (d) (e)It can be seen from Figure13 that turbine location at both end and middle of channel has a certain effect on surrounding flow field, and a certain range of velocity deceleration zones which are “8-shaped” areas appears in facing flow side and reverse flow side of turbine, and flow velocity on two sides of turbine increases slightly. In the process of spring tide, difference of flow velocity in near flow field has largest change, and affecting area is biggest when B4 berthage is installed with turbine, B3 berthage is next to B4 berthage, B5 berthage is followed by B3 berthage, B1 berthage is second to B5 berthage, and B2 berthage is smallest. B4 and B3 berthage are in narrow section of channel, and B2 berthage is located in open section of channel. This indicates that the effect of installation of turbine on flow field in narrow sea area is relatively larger, and the effect of installation of turbine on flow field in open sea area is relatively smaller. In addition, water depth of B2 berthage is largest in the 5 berthages, but the effect on the flow field is the smallest. Affecting range of B4 and B3 berthage in channel on flow field is substantially consistent with waterway direction, indicating that the effect of installation of turbine on flow field is also related to water depth, topography, and channel direction.Hourly changing process of extreme value of deceleration and relative difference of flow velocity in different exploratory circumference with radius of 100 m, 200 m, and 300 m are respectively shown in Figures14, 15, and 16. Wherein maximum of deceleration represents maximum value of decrease of flow velocity at every point on circumference at a certain time, definition of minimum of deceleration is just the opposite. The relative difference of flow velocity represents the ratio of the maximum of deceleration to the maximum of flow velocity before the installation of the turbine.Figure 14
Variations in flow velocity of circumference with radius of 100 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)Figure 15
Variations in flow velocity of circumference with radius of 200 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)Figure 16
Variations in flow velocity of circumference with radius of 300 m: (a) and (b) represent B1 location; (c) and (d) represent B2 location; (e) and (f) represent B3 location; (g) and (h) represent B4 location; and (i) and (j) represent B5 location.
(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)As can be seen from these figures, the extremal curves of deceleration at each berthage show a certain periodicity, which is the result of the periodicity of the tidal current. The curves of B2 berthage are special, and the curve patterns of other berthages are similar. Most of the deceleration minimum curves of all circumferences are located above the horizontal axis and most of the deceleration maximum curves are below the horizontal axis. This shows that there is a certain degree of deceleration and growth of flow velocity appearing near the turbine, which is also consistent with the results of the 3-D near-field simulation. The peak value of the deceleration maximum curves of each turbine generally decreases with the increase of the circumferential radius, which means that the farther away from the turbine, the smaller the effect of the turbine on the flow field. From the curves of relative difference of flow velocity, the effect of the turbine on the far-field of flow field is smaller as a whole, and the maximum change in flow velocity does not exceed 8% of the flow velocity before the installation of the turbine. The degree of deceleration and growth of all berthages were compared as shown in Table4.Table 4
Comparison of deceleration and growth degree.
Circumferential Radius Comparison of Deceleration Degree Comparison of Growth Degree 100 m B4> B3> B1> B5> B2 B3> B2> B1> B5> B4 200 m B4> B3> B1> B5> B2 B3> B2> B1> B4> B5 300 m B4> B3> B1> B5> B2 B4> B2> B3> B1> B5As can be seen from Table4, the results of comparison of the deceleration degree (B4> B3> B1> B5> B2) are approximately the same as the comparison of the maximum value distribution of the flow velocity difference (B4> B3> B5> B1> B2). For the growth degree, B4 berthage changes greatly. There is no growth in its 100 m range, but the degree of growth near 300 m is the largest. In order to better reflect the specific changes in the affected area, the change of the extremum of flow velocity difference is formed by different location under the spring tide, as shown in Table 5.Table 5
Extremum of flow velocity variation.
Location Radius 100 m Radius 200 m Radius 300 m Deceleration Growth Relative Difference Deceleration Growth Relative Difference Deceleration Growth Relative Difference (m/s) (m/s) (%) (m/s) (m/s) (%) (m/s) (m/s) (%) B1 -0.023 0.006 1.758 -0.022 0.007 1.754 -0.016 0.006 1.391 B2 -0.001 0.017 0.419 -0.001 0.011 0.226 -0.003 0.009 0.281 B3 -0.064 0.034 4.754 -0.043 0.024 3.354 -0.034 0.007 3.147 B4 -0.091 0.002 7.461 -0.088 0.01 6.877 -0.053 0.024 3.947 B5 -0.019 0.002 1.503 -0.012 0.003 0.953 -0.009 0.003 0.714It can be seen from Table5 that B4 location has the largest effect on the flow field and the affected area is also the most extensive, and B2 location has the least effect on the flow field. The deceleration maximum of B1, B3, B4, and B5 decreases with the increase of radius, the change of B2 location is smaller, and there is fluctuation around the radius of 300 m. The growth maximum of B2 and B3 decreases with the increase of radius, and the growth changes of B1, B4, and B5 are smaller and the fluctuation also appears. This shows that when the degree of deceleration and growth of flow velocity is greater, the value of deceleration and growth will decrease with the increase of radius. When the deceleration and growth of velocity are less than 0.01 m/s, the value of deceleration and growth will only fluctuate and there is no obvious law.To sum up, for the installation of tidal turbine in tidal farm, the following conclusions can be obtained by combining the consideration of the bathymetric data and the tidal resource in the engineering sea area:( 1 ) Compared to other berthages, the effect of test point B4 on the flow field is the largest after turbine’s layout; this point can be considered to avoid in the layout of berthages.( 2 ) Tidal current power resources of B1, B3, and B5 test points are more abundant; the degree of effect after turbine’s layout is also within the acceptable range.( 3 ) Combining the offshore distance of each point, test points B1 and B3 are the ideal layout points for the five berthages through synthetical consideration.


## 5. Conclusion

( 1 ) 3-D CFD simulation results show that the tidal turbine has a certain degree of hindrance to the upstream inflow, and the faster the rotating speed of the impellers, the more obvious the hindrance. After water flow through the turbine, flow velocity decreases significantly, and a strip of wake flow with the length of about 240 m (about 10DE) is formed downstream of impellers. The distribution of variation gradient of flow velocity in the strip of wake flow is not well-distributed, and the closer the velocity of the impeller is, the more severe the loss is. The wake flow of the double impellers is gradually approaching downstream, converging into one, then narrowing and disappearing in the distance. Flow velocity of inner and outer flow downstream the double impellers is greater than the inflow velocity because the flow section is narrowing.( 2 ) The tidal level which is simulated by 2-D tidal model is in good agreement with the measured values, the flow velocity distribution of the tidal farm is reasonable, the calculated values of the flow velocity and the flow direction are consistent with the measured values, and the variation trend is consistent, which indicated that the mathematical model of the tidal current is reliable and can be used to simulate the flow field of the channel of Putuo-Hulu Islands.( 3 ) The simulation results of the 2-D tidal mathematical model show that turbine’s installation at both ends and the middle of the channel has a certain effect on the surrounding flow field, a certain range of velocity deceleration zones which are “8-shaped” areas appears in the facing flow side and reverse flow side of the turbine, and the flow velocity on two sides of the turbine increases slightly. The water area at both ends of the channel is wide and affected little by the turbine. The area where the flow velocity is affected is larger in the middle of the channel because of the superimposed effect of the narrowed area and the turbine. The size of the affected area is also related to the inflow velocity and the topography; meanwhile, the reciprocating flow of the tidal current is one of the reasons for the expansion of the affected area.( 4 ) From the statistical results of the relative velocity difference, the effect of the turbine on the flow velocity of far-field is relatively smaller. The effect of turbine on the far-field of flow field is smaller as a whole; the maximum change in flow velocity does not exceed 8% of the flow velocity before the installation of the turbine.( 5 ) Compared to other berthages, the effect of B4 berthage on the flow field is the largest after turbine’s layout; this point can be considered to avoid. Tidal current power resources of B1, B3, B5 berthages are more abundant; the degree of effect after turbine’s installation is also within the acceptable range. Combining the offshore distance of each berthage, B1 and B3 berthages are the ideal location for the five berthages through synthetical consideration.



---

*Source: 7983907-2019-01-09.xml*