The Numerical Simulation of Hydrodynamics of Fishing Net Cage

The aquaculture industry is playing an increasingly important role in the fish production industry as the demand for seafood increases. At present, the net cage is widely used in aquaculture industry all over the world. Because of environmental impact concerns and limited near-shore locations, more and more net cages for aquaculture will be located offshore and exposed to more waves and currents. Thus, knowledge of their hydrodynamic behavior of net cage under the action of waves and current is of great importance to the design of gravity cages in the open sea.


Introduction
The aquaculture industry is playing an increasingly important role in the fish production industry as the demand for seafood increases. At present, the net cage is widely used in aquaculture industry all over the world. Because of environmental impact concerns and limited near-shore locations, more and more net cages for aquaculture will be located offshore and exposed to more waves and currents. Thus, knowledge of their hydrodynamic behavior of net cage under the action of waves and current is of great importance to the design of gravity cages in the open sea.
Techniques used to investigate the facilities have typically included the use of scaled physical and numerical models, and, where possible, field measurements. Comparing model tests and field measurement, the numerical simulation method is low in cost, saving time and being easy to manage. In this chapter, the latest progress of this technology by our research group will be introduced in detail. The contents of the progresses include the following parts: The numerical simulation of net cage in irregular waves; The numerical simulation of multi-net cage in waves; The numerical simulation of flow through and around fishing net cage.

Numerical methods
The fishing net cage is mainly composed of two parts: floating collar and fishing net. The floating collar is model as a rigid body and the fishing net is simulated as flexible body. The numerical model of net cage system is described in detail as follow.

Fishing net
By applying the lumped-mass model, the net is assumed to be a connected structure with springs and limited mass points, which are set at each knot and the centre of the mesh bar, as shown in Fig. 1 (1) where D F  and I F  are the drag and inertia forces, respectively, R   is the acceleration of mass point, T  is the tension force in twine, B  is the buoyancy force, and W  is the gravity force. According to the research of Zhao et al. (2008), the inertial force I F  on fishing net in waves is much smaller than drag force, so here it is omitted. ( 2) where T is the tension force in twine, l 0 is the undeformed length of twine, l is the deformed length, d is the diameter of twine, C 1 and C 2 are the elastic constants of material which can be obtained by matching, referring to Gerhard (1983). For polyethylene (PE), C 1 =345.37×10 6 , and C 2 =1.0121; for polyamide (PA), C 1 =784.9×10 6 , and C 2 =1.6988. The units of T and d are N and m, respectively.
To consider the direction of fluid forces acting on net mesh bar, a local coordinate system Oτηξ is defined, as described by Fig. 2. The origin of the local coordinate system is set at the centre of a mesh bar and the η axis lies on the plane including τ and V.
Besides the gravity and buoyancy, the τ component of the drag force on lumped-mass point can be obtained by Morison Equation, as follows.
where C D represents the drag coefficient of the τ component, D is the diameter of twine, and l is its length. The same expression can be applied to the other drag forces (F Dη , F Dξ ) of η, ξ components. The drag coefficients related to Reynolds number are described in following section. www.intechopen.com where Re n = V Rn D/ , s=-0.077215665+ln(8/Re n ), is the viscosity of water, C n and C are the normal and tangential drag coefficients for mesh bar, V Rn is the normal component of the fluid velocity relative to the bar, and is the density of water.
For the knot part, the Fredheim and Faltinsen (2003) suggested it could be reasonable to use a drag coefficient in the range of 1.0-2.0 when modeling the knot part as a sphere. Here C D is set as 1.0 for the knot part.

Forces on floating collar
Although the float collar mainly consists of two concentric pipe rings and a hand rail. Float collar system is usually at the water surface and double floating pipes are the main components to withstand the wave and current induced loads. For the sake of simplicity, the float collar system is reduced to a double-column pipe structure.
When calculating the forces on float collar, the float collar is divided into many minisegments. The forces on the whole collar can be obtained by summing the forces on each mini-segment. Fig. 3 is a sketch of a mini-segment of the float collar with a local coordinate www.intechopen.com Hydrodynamics -Theory and Model 290 system defined in each mini-segment. As to the coordinate system, n and τ are in the normal and tangential directions of the mini-segment, respectively, and then v is normal to the miniplane.
Because the tube diameter of float collar is relatively small compared to the wave length. As Brebbia and Walker (1979) suggested, the n component of wave-induced forces on a minisegment can be calculated using the modified Morison equation to include relative motion between the structural element and the surrounding fluid, which is shown as follows: where n u  and n R   are the velocity vectors for water particles and mini-segments of the n component, respectively; n a  and n R   are the acceleration vectors for water particles and minisegments of the n component, respectively; is the density of water; V 0 is the water displaced volume of a mini-segment; A n is the effective projected area of a mini-segment in the direction of the n component; and C Dn and C mn are the drag and added mass coefficients of the n component, respectively. The same expression can be applied to other waveinduced forces (F , F ). The relationship between the mooring-line forces and elongation is determined experimentally and to be used in the numerical simulation, as follows: where ΔS is the elongation of a mooring line (m), S is the original length of the mooring line (m), and T is the tension in the mooring line (N).
Apart from the above-mentioned external force, the float collar is also subjected to gravity and buoyancy forces. Because it is easy to calculate these forces, the formulation for these forces is not given here.

www.intechopen.com
The Numerical Simulation of Hydrodynamics of Fishing Net Cage 291

Motion of floating collar
The three-dimensional motions of the float collar include surge-sway-heave translation and roll-pitch-yaw rotation. To obtain the motions of the float collar, two sets of coordinate systems are defined, which are the fixed-coordinate system Oxyz, and the body-coordinate system G123, as shown in the Fig. 4. Initially, axes x, y and z are parallel to axes 1, 2 and 3.
where F xi , F yi and F zi are the components of the external forces on mini-segment along fixedcoordinate axes x, y and z, m G is the mass of the float collar, x G , y G and z G are the acceleration of the mass centre of the float collar, and N is the number of mini-segments.
Axes 1, 2, 3 are principal axes with origin at the centre of the mass G, and thus the Euler equations of motion of a rigid body (Bhatt and Dukkipati, 2001) are applied. In the bodycoordinate system, the three rotational equations of motion are given by: where subscripts 1, 2, 3 represent the body-coordinate axes 1, 2, 3; I 1 , I 2 and I 3 are the components of the moments of inertia I along the three principal axes; ω 1 , ω 2 and ω 3 are the components of angular velocity vector ω along the three principal axes; M 1 , M 2 and M 3 are the components of the moment vector M along the three principal axes.

Flow around fishing net
The numerical simulation of the flow field around a plane net is based on the FLUENT software platform. The porous media model is introduced to model the plane net (see Fig.  5), and the finite volume method (FVM) is used to solve the governing equations of the numerical model. In this way, the numerical simulation of the flow field around the plane net is available.

Porous media resistance coefficients
The porous media model employs empirically determined flow resistance in the region of the porous media (Fluent, 2006). For flow through the porous media, the hydrodynamic forces acting on the porous media can be expressed as follows: where S i is the source term for the momentum equation in the i direction, is the thickness of the porous media, A is the area of the porous media, and F is the hydrodynamic force in the i direction.
When the region is outside the porous media model, S i = 0. While, inside the porous media, S i is calculated by the following equation: where D ij and C ij are prescribed material matrices consisting of the porous media resistance coefficients, D n is the normal viscous resistance coefficient, D t is the tangential viscous resistance coefficient, C n is the normal inertial resistance coefficient, C t is the tangential inertial resistance coefficient.
Substituting Eq. (11) into Eq. (10) provides formulas for calculating the drag force (F d ) and the lift force (F l ) of the plane net. The drag force is parallel to the flow direction, and the lift force is perpendicular to the flow direction. www.intechopen.com The porous coefficients in Eq. (11) can be calculated by the data for the drag and lift forces. In a general way, the drag and lift forces of the plane net are obtained from the laboratory experiments. In addition, the forces can be calculated from the Morison equation: where C d and C l are coefficients that can be calculated using empirical formulas proposed by When the plane net is oriented normal to the flow, the porous coefficients (D n and C n ) are chosen from a curve fit between drag force data of the plane net and corresponding current velocities using the least squares method. The other two coefficients (D t and C t ) can be ignored because the lift force is equal to 0 when =90°.
When the plane net is oriented with different attack angles (see Fig. 6), the porous coefficients should be transformed into formula (16) (Bear, 2006). Then, the four porous coefficients can be obtained by minimizing the error between the theoretical values and existing data for the drag and lift forces using the least squares method, which is a common analytical method in error minimization.

Mesh grids and boundary conditions
An example of computational grids of a plane net at an attack angle =90° is shown in Fig. 7. The mesh consists of unstructured tetrahedral elements that are refined in the vicinity of the porous media. The coordinate system for the numerical model is a right-handed, 3D Cartesian coordinate system, where x is positive toward the flow direction. The left boundary of the flume tank is described by the velocity-inlet boundary condition, while the right boundary is described by the outflow boundary condition. The solid surfaces of the flume tank and the free surface are modeled using the wall boundary condition (with zero shear force).

Hydrodynamic simulation of net cage in irregular waves
The real sea is mainly composed of irregular waves. In this section, the hydrodynamic behavior of net cage in irregular waves is analyzed by numerical simulation and physical model test.

Irregular wave model
The linear wave theory (refer to Zhao et al. 2007) and random phase method are employed in the computation of irregular wave field. Based on linear wave theory, the velocity potential of the wave is given by where A is the wave amplitude, g is gravitational acceleration, f is the wave frequency, k is the wave number (equal to 2 /L ), L is the wave length, d is the water depth, z is the vertical distance (positive upward) from mean-water level.
The irregular waves are generated using the random phase method as described by Charkrabarti (1994) which, in the time domain, is the superposition of multiple regular waves.
where A j is given by where A j and K j are the wave amplitude and number of the individual wave components, respectively. ε j is the uniform distributive random phase between 0 and 2 . The surface elevation time series will repeat after a given amount of time, and it is correlated with the two aspects: the number of regular waves used in the superposition and the sample length of the surface elevation. In the present study it was decided to discretize to 60 waves in the spectrum and have a length of each time series giving approximately 100 wave cycles. The input wave spectrum S(f) is determined by the modified JONSWAP spectrum given by Goda ( where T P is the spectral dominant period, T p =T H1/3 /(1-0.132( +0.2) -0.559 ); H 1/3 and T H1/3 are the significant wave height and period, respectively; f is the wave frequency; f P is the spectral peak frequency; is the peak shape factor, =0.07 (f ≤ f P ), =0.09 (f > f P ), is the peak enhancement factor and equal to 3.3. The parameter J is determined by the following expression: 1 0.06238 1.094 0.01915ln 0.230 0.0336 0.185(1.9 )

Physical model set up
A series of experiments of the motion response of float collar and the tension response in the mooring line were conducted in a wave-current flume at the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, China. The wave-current flume is 69m long, 2m wide and 1.8m high, equipped with a random wave-maker and a currentproducing system. The test setup is shown in Fig. 8. As shown in the Fig. 8, the physical model was installed in the centre of the flume. The mooring-line forces were measured by two transducers attached to the bottom of the mooring lines. The anchors on the left side are labeled for load identification. An optical measurement system for determining the float collar motions was developed. Two diodes (front and back points), numbered 1 and 2, were fixed on the float collar for motion analysis. The movement of diodes was recorded by a CCD camera. The camera captures a series of images and transfers each frame to the computer and temporary storage. Later, specially-written software is used to search the position of diodes, which is used to calculate the float collar horizontal and vertical movements as a function of time. The geometric and mechanical properties of gravity cage model are provided in Table 1

Comparison in time domain
The numerical and physical gravity cage models were subjected to regular and irregular waves. The motion response (in surge and heave) of the float collar and the tension response in the mooring line are analyzed and compared in the time domain.

Fig. 9. Comparisons of the average of top one-tenth mooring line forces between numerical and experimental data
To examine the numerical model, quantitative comparisons are carried out between calculated and experimental results. The average of one-tenth highest mooring line tensions and cage motions are selected for comparisons. Fig. 9 shows the comparisons of the average of one-tenth highest mooring line tensions between calculated and experimental results. Fig.  10 and 11 demonstrate comparisons of the top one-tenth cage surge and heave motions between simulated and experimental results, respectively. The mean relative error of the mooring-line forces is 6.9%, and the mean relative errors of the cage surge and heave motions are 9.9% and -4.2%, respectively.

Comparison in frequency domain
For regular wave cases, response amplitude operators (RAOs), determined by dividing the amplitude of the response by the amplitude of the wave, are used to characterize the surge and heave motion response and the mooring-line tension response. For irregular wave cases, the transfer functions of the motion response of the float collar and the tension response in the mooring line were calculated using auto-spectral density method, and compared with the regular wave RAOs. The corresponding transfer functions are calculated as a function of frequency. Combining with statistical wave descriptions, these transfer functions can be used to yield statistics of cage dynamic response in a realistic sea conditions.
The time series of wave elevation and heave motion of float collar from the numerical simulation and physical model tests are shown in Fig. 12. In Fig. 12, the heave motion of float collar is synchronized with the wave elevation. Based on these time series, the corresponding auto-spectral density function can be obtained by using Fast Fourier Transform (FFT) method. The numerical and physical model wave spectra are shown in Fig.  13. The Fig. 13 shows that the wave spectra obtained from the numerical simulation and physical model test are similar to each other.
According to Bendat and Peirsol (1986), the linear transfer functions of cage motion and tension responses can be calculated using the auto-spectral technique.
Fig. 14 shows the heave motion transfer function. In general, the motion response transfer function results for heave were characteristic of a highly damped system. At high frequencies, the cage system has no significant heave motion response to wave forcing. With the decreasing of wave frequency, the heave motion transfer function increases approaching one, indicating wave contouring behavior. This is reasonable, and it agrees with our qualitative analysis according to the common knowledge. In all cases the heave motion transfer function did not exhibit a resonant peak but tend to a value of one as the wave frequency decreases. If a resonant condition exists, it most likely occurs at the waves with frequency lower than 0.6 Hz. The surge transfer function results are presented in Fig. 15. It shows that the responses of the gravity cage in surge, obtained from the numerical simulation and physical model test, are similar with little response in high frequencies. As with heave motion response, the surge motion tends to increase slightly with decreasing wave frequency.
As shown in Fig. 16, the mooring line tension results, obtained from numerical simulation and physical model test, matched up well throughout the frequencies. The transfer functions of the mooring line tension and the cage motion in heave and surge show the same trend.
As the wave frequency increases, the tension response in mooring line decreases. The numerical and physical models both predicted a mooring line tension response less than 0.02 N/mm throughout the wave frequencies.

Hydrodynamic behavior of multiple net cages
In the open sea, the large fish farm often includes multiple net cages assembly together. In this section, a numerical model of four-cage mooring structure in regular waves is presented. The effect of wave directions on mooring line tensions is analyzed.

Cage arrangement description
Mooring system is used to hold cage structures against the forces caused by waves and currents. Its design is very important for the operating and performance of cage structures. A change in the mooring system will change the internal loads on the cage systems. In this section, a four-cage system with grid mooring is analyzed. As shown in Fig. 17, four cages are arranged in two columns. Each net cage is connected to the submerged grid mooring by four bridle lines. The submerged grid mooring is attached to the sea floor by anchor lines. The direction of anchor line is perpendicular to grid line, which forms a 16 degree angle with the sea floor. The grid mooring consists of four submerged, pre-tensioned square grid (1 m × 1 m). The horizontal grid is 0.1 m below the still water surface. The pretension on the mooring system is maintained by the submerged floater at each grid corner. The diameter of floating pipe is 6.25 mm, and its density is 7.1 g/m. The diameter of net twine is 0.72 mm, and the mesh size is 11.7 mm.

www.intechopen.com
The Numerical Simulation of Hydrodynamics of Fishing Net Cage 301 Fig. 17. Offshore grid mooring cage structure

The effect of wave direction on mooring line tension
Note that the waves propagated along various directions in the open sea, and it will affect the hydrodynamic behavior of grid mooring cage structure. Therefore, three different wave directions are considered in this section. Concerning the symmetrical characteristics of grid mooring cage structures, wave directions of 0°, 30° and 45° are chosen. Since the tension force on anchor lines is largest among the three types of mooring lines, only the maximum tension force on anchor lines is given here for the sake of clarity. Fig. 18 shows the maximum tension force on anchor lines of net cage system under the action of waves from different directions. For three different wave directions, the maximum tension forces on anchor lines are 0.375 N, 0.315 N and 0.282 N, respectively. When the wave incident angle is 0°, the left three anchor lines are main components to withstand external force on net cages. If the incident angle becomes 45°, the left three anchor lines and the bottom three anchor lines are main components to withstand the external force. Considering three kinds of wave directions, the uniformity of tension force on anchor lines is different.
When wave incident angle becomes 45°, there are more anchor lines to withstand external force on net cage. Therefore, the maximum tension force on anchor lines is smallest when wave incident angle is 45°.
For different wave incident angles, the transfer load path is different, and hence the uniformity of tension force on anchor lines is different. In the design of mooring cage structure, wave incident angle of 45° may be a good choice.

Numerical simulation of flow through and around fishing net cage
In this section, the flow field inside and around a gravity cage is simulated. The diameter of the cage (D) is 16 m, and the height (H) is 10 m. The cage net is knotless nylon net. The twine diameter was 2.8 mm, and the net mesh size was 29 mm. For the numerical model only the fishing net is modeled, and the effect of the float collar and sinker system on the flow field is ignored. The cylindrical cage is divided into 16 plane nets around the circumference and a bottom net (see Fig. 19). These plane nets can be described as the porous media with different attack angles. Therefore, the simulation of the flow field inside and around a gravity cage is straight forward. Fig. 19. Sketch of the gravity cage model.

Numerical model description
In the numerical model, the flume is 240 m long, 80 m wide and the depth of the water is 20 m. The boundary is considered large enough that the influence of the boundary condition on the flow field around the cage is negligible. The thickness of the porous media is 50 mm, and the porous coefficients are D n =75083 m -2 , D t =38307 m -2 , C n =4.985 m -1 , and C t =1.660 m -1 . The numerical simulations are performed with the incoming velocity u 0 =0.5 m/s and the inlet turbulence quantities: k=1.34×10 -4 m 2 /s 2 and ε=4.76×10 -9 m 2 /s. Four cases are modeled with different cage numbers, n=1, 2, 3 and 4, and the spacing distance between two adjacent cages is 24 m. The cages, which are located in a row along the flow direction, are centered on the width of the flume. The front cage is positioned 40 m downstream from the velocity-inlet boundary, while the others are positioned downstream successively with different cage numbers. Fig. 20 shows the sketch of one cage model.

Simulated results
As shown in the contour plots, the velocity is approximately uniformly distributed inside the cage. There is a small region of velocity reduction upstream of the cage, while there is a rather large velocity reduction region downstream from the cage. Obvious flow velocity reduction exists along the flow direction, and the velocity reduction increases with increasing cage numbers. The width of the wake becomes narrower toward the centerline with increasing distance from the cage. The flow velocities outside the velocity reduction region are approximately 1%-10% greater than the incoming flow u 0 , and more cages lead to greater flow velocity around the cages (see Fig. 21). As the cage number increases from 1 to 4, the maximum velocity reduction downstream from the cages increases from 25% to 69%. The cages influence the flow velocity distribution inside and around them. The rear cage has no influence on the flow velocity distribution inside the former cage (see Fig. 22).
where r i is the flow velocity reduction factor, n c is the number of upstream crossings of the cage net, and C d is the drag coefficient of the net.
The flow velocity reduction factors u/u 0 show a good agreement and the velocity reduction trend of the numerical simulations is consistent with the results of Aarsnes's formula at different cage numbers. Fig. 23 shows the comparison results when n=4.

Conclusions
Three numerical models for calculating net cage in irregular waves, multi-net cage in waves and flow through fishing net cage are developed in the paper. The simulated results are in good agreement with the experimental values. Some main conclusions can be drawn as follows: 1. The gravity cage system shows a characteristic of a highly damped system with small response to high frequencies. The gravity cage surge and heave motions generally decrease as the wave frequencies increase. The heave motion is synchronized with the wave elevation. At low frequencies, the heave motion transfer function is approaching to one indicating wave contouring behavior. The transfer functions of cage motion response (heave and surge) and mooring line tension response have the same trend. 2. When the wave incident angle is 45°, the maximum tension force on anchor lines is smallest. Therefore, in the design of grid mooring cage structure, wave incident angle of 45° may be a good choice. 3. The velocity is approximately uniformly distributed inside the cage. As the cage number increases from 1 to 4, obvious flow velocity reduction exists along the flow direction, and the velocity reduction downstream from the cages increases from 25% to 69%. However, the rear cage has no influence on the flow velocity distribution inside the former cage. www.intechopen.com