Development of a porous media model with application to flow through and around a net panel

ARTICLE IN PRESS Ocean Engineering 37 (2010) 314–324

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Development of a porous media model with application to flow through and around a net panel Øystein Patursson a,b,, M. Robinson Swift c, Igor Tsukrov c, Knud Simonsen b, Kenneth Baldwin a, David W. Fredriksson d, Barbaros Celikkol c a

Jere A. Chase Ocean Engineering Center, University of New Hampshire, Durham, NH 03824, USA Faculty of Science and Technology, University of the Faroe Islands, FO-100 Torshavn, Faroe Islands c Mechanical Engineering Department, Kingsbury Hall, University of New Hampshire, Durham, NH 03824, USA d Department of Naval Architecture and Ocean Engineering, U.S. Naval Academy, 590 Holloway Road, 11D Annapolis, MD 21402, USA b

a r t i c l e in f o

a b s t r a c t

Article history: Received 18 February 2009 Accepted 1 October 2009 Available online 9 October 2009

The flow characteristics through and around a net panel have been investigated through computational fluid dynamics (CFD) and measurements. A finite volume approach was used for solving the Reynolds averaged Navier–Stokes equations combined with a k2e turbulence model for describing the flow. For computational efficiency, the net was modeled as a sheet of porous media rather than a large number of cylinders connected by knots. The model resistance coefficients needed for the porous media equations were found by optimizing the fit between computed lift and drag forces on the net panel and lift and drag measured in tow tank experiments. Lift and drag acting on a flat panel of knotless nylon net (2.8 mm twine thickness and 29 mm mesh size) stretched on a frame were measured at different speeds and angles of attack, and fluid velocity was recorded in the region behind the net. The optimization process used to obtain the best fit porous media coefficients was simplified through the use of an analytical model. Final comparisons between CFD predictions and measurements of lift and drag coefficients and velocity reduction behind the net panel were made for two of the speeds and angles of attack. The agreement between measured and modeled data was good with a mean normalized absolute error of 6.2%. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Computational fluid dynamics Porous media model Net panel Current reduction Tow tank experiments

1. Introduction The flow field inside and around sea-cages used in aquaculture is important both from an engineering perspective, due to the effect of fluid velocity on drag force, and from an ecological perspective with regards to exchange of water in the cages. The flow characteristics can be described by a combination of computational fluid dynamics (CFD) modeling, where the momentum and continuity equations are solved numerically, and measurements. Since cages are mostly net a numerical description of this component is necessary. The netting can be approximated by small cylinders connected in knots, but the number of cylinders required to model the net in a cage can be substantial—more than 10 million in a single salmon farming cage, thus it is computationally expensive to model the exact geometry. Another factor is that when the cages have been in the water for some time the nets are biofouled and the resulting complex geometry can have greatly reduced permeability. A general

 Corresponding author. Present address: Aquaculture Research Station of the Faroes, Faroe Islands. Tel.: +298 474747; fax: + 298 474748. E-mail address: oystein@fiskaaling.fo (Ø. Patursson).

0029-8018/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.oceaneng.2009.10.001

approach which addresses these issues is to model the net as a sheet of porous media. In the present work 3D CFD is used where the net is described as a thin volume with added resistance governed by the equations for modeling flow through porous media. A 2D study using the same technique was reported in Patursson et al. (2006) while Vincent and Marichal (1996) used a similar but different description of the flow through the net for modeling axisymmetric structures made of net. Helsley and Kim (2005) also use CFD for investigating mixing behind an aquaculture cage. Other approaches including analytical models have also been used (e.g. Løland, 1991; Gui et al., 2006; O’Neill, 2006). In this study, tow tank measurements on a net panel are used to ‘‘calibrate’’ the model, and to demonstrate the capabilities of the model. The basic shape used by most researchers is a flat net panel generally stretched on a stiff frame (Løland, 1991; Le Bris and Marichal, 1998; Zhan et al., 2006). The net panel was set up in a tow tank and towed at different speeds and oriented with different angles of attack. Measuring current reduction behind the net panel and drag and lift forces on the net panel as a function of tow speed and angle of attack provided a good basis for evaluating the capabilities of the model to simulate the flow through and around a net.

ARTICLE IN PRESS Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

The ‘‘calibration’’ phase consisted of finding the best porous media resistance coefficients for the specific net used in the tow tank experiments. One objective was to develop a general methodology that can be applied to any net of interest. Because CFD runs are time consuming, a simplified analytical method was derived to determine the porous resistance coefficients of the net from measurements of drag and lift forces as a function of angle of attack and tow speed. The method included an analytical description of the drag and lift forces on the net panel, applying the equations for flow through porous media and appropriate approximations. Using a fitting procedure, the porous resistance coefficients that minimize the error between measured and calculated drag and lift forces were found. To validate the approach, CFD predictions using the best-fit porous resistance coefficients were then compared with measured drag and lift forces and current reduction. In order to check the correctness of the analytical method for finding the porous resistance coefficients, the effect of increasing and decreasing each of the coefficients was investigated.

2. Theory The equations describing the flow in the computational domain are the Reynolds averaged Navier–Stokes (RANS) equation   Dui 1 @P @ @ui @uj 1 þ Si ¼  þ gi þ neff þ ð1Þ Dt r @xi @xj @xi r @xj and continuity equation @ui ¼0 @xi

ð2Þ

ð3Þ

where m is dynamic viscosity and mt is the turbulent viscosity or eddy viscosity. Summation over repeated indices is assumed throughout this paper. When using one of the k2e models, mt is a function of k and e, which are the turbulent kinetic energy describing the energy in the turbulent eddies and the turbulent dissipation rate, describing the dissipation rate of k. Equations describing k and e are needed for closure of the equation system. These equations are usually generated by forming scalar transport equations for the wanted variables as described by e.g. Launder and Spalding (1972) and Shih et al. (1995). For this study the model proposed by Shih et al. (1995) was used. The flow through the net is described using the equation for flow through a porous media developed from the equations proposed by D’arcy (1856) and Forchheimer (1901) Si ¼  Dij muj  12Cij rumag uj

and 2

Cn 6 0 Cij ¼ 4 0

0 Ct 0

3 0 07 5 Ct

ð6Þ

If the local axes of the porous media are not aligned with the global coordinate axes, a tensor rotation approach has to be employed. These general equations were solved using the CFD software FLUENT 6.3 from ANSYS Inc. The volume integrated steady state versions of Eqs. (1) and (2) were discretized using a finite volume approach and solved using the segregated iterative solver employing an algebraic multigrid method with under relaxation. The momentum equation was discretized using a second order upwind scheme, derivatives were evaluated using the node based formulation, and the SIMPLE (e.g. Ferziger and Peric´, 2002) method was employed for the pressure–velocity coupling. The grid was generated by the meshing software Gambit 2.3 from ANSYS Inc., and two grid generation methods were employed—TGrid that generates meshes from tetrahedral cells and Hex Core that generates mixed meshes using tetrahedral and hexahedral cells. The velocity at the inlet was constant across the inlet face and values of k and e were specified in the inlet water. Wall boundaries were modeled using the standard wall functions. More detailed information about the software is provided by Fluent (2006, 2007).

3. Measurements

where ui are the fluid velocity components, r is the fluid density, xi are the spatial coordinates, gi is the acceleration due to gravity, Si is a source term used to describe the added resistance of the 02 0 is the root mean net, P ¼ p þ 23 rk, p is pressure, k ¼ 32 urms , urms square of the turbulent velocity fluctuations, neff ¼ meff =r and

meff ¼ m þ mt

315

ð4Þ

where umag is the magnitude of the fluid velocity and Dij and Cij are prescribed material matrices consisting of the porous media resistance coefficients in local x1 , x2 and x3 direction. The local axes are the principal axes of the porous media. If x1 is normal to the net panel and x2 and x3 are normal to each other, but parallel to the plane of the net panel, the Dij and Cij matrices are of the following form: 2 3 0 Dn 0 6 0 D 0 7 Dij ¼ 4 ð5Þ t 5 0 0 Dt

A series of measurements of drag and lift forces on a net panel and velocity reduction behind the net panel were conducted in the University of New Hampshire (UNH) 37 m long, 3.66 m wide and 2.44 m deep tow/wave tank as described by Patursson (2007). The test setup shown in Fig. 1 was constructed and used in the tank to support the net panel and measure the forces. Drag force (parallel to velocity) and lift force (perpendicular to velocity) were measured using a loadcell arrangement. The water speed relative to the carriage and net was measured using an acoustic doppler velocimeter (ADV) in different positions behind the net, mostly in the wake region, but also outside the wake. The velocity of the carriage was measured using a set of light gates measuring travel time over a set distance. 3.1. Measurement setup The net panel was positioned well below the water surface (0.73 m) to minimize drag forces from surface wave creation due to the pressure difference between the front and the back of the net panel. Based on preliminary CFD investigations, it was decided to position the net panel in the center of the cross-section of the tank. The CFD investigations also indicated that the wave height of the surface wave created by the net in this position would be on the order of a few millimeters and can be considered negligible. The setup was designed such that the structure deformation was small and the angle of attack of the net panel could be set within 7 23 ( 7 13 from setting the angle of attack and 713 from bending of the structure). The basic components were a frame that holds the net panel, which was attached to a strut, and a loadcell arrangement that was attached below the carriage close to the water surface (Fig. 1). The strut passed through the loadcell arrangement, and at the upper end it was connected to a universal joint that allowed bending but no swiveling. The two loadcells were positioned perpendicular to each other, one in the drag direction and the

ARTICLE IN PRESS 316

Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

Fig. 1. Overview of the setup for measuring drag and lift forces on the net panel. The dimensions are in meters. Two loadcells were used—one loadcell in the tow direction for measuring the drag force and one perpendicular to the tow direction for measuring the lift force. The net panel was attached to a frame of 26.7 mm OD stainless steel pipe. The frame was attached to a strut passing through the loadcell arrangement and terminated in a universal joint allowing bending but no swiveling. The loadcells were attached to the strut using pinned joints to minimize bending of the loadcells.

Table 1 Accuracy of ADV measurements according to data from the supplier. sADV ðcm s1 Þ

Vel. range (cm s1 )

edat ðcm s1 Þ

edat (%)

12.5 25 50 75

10 10 30 100

0.35 0.35 0.55 1.25

2.80 1.40 1.10 1.66

sADV is the measured speed, Vel. range is the velocity range setting of the instrument and edat is the accuracy.

The net used in the experiments was a 1 m by 1 m knotless nylon net panel with twine diameter d ¼ 2:8 mm (210d/108) and mesh bar length (length of one mesh bar between twine intersections) l ¼ 29 mm when it was stretched onto the frame. The net was attached to the frame using the square mesh orientation. Net solidity was defined as twine projected area divided by outline area and calculated by the equation used by Fredheim (2005)  2 2d d  : ð7Þ S¼

l

other in the lift direction. The strut could be oriented arbitrarily around the vertical axis, and the angle of attack of the net panel could be set using a protractor on top of the universal joint. The loadcell arrangement was calibrated in place by pulling from the center of the net panel and measuring the force with a calibrated loadcell. A slightly nonlinear calibration curve was found. The accuracy of the force measurement was 7 0:5 N, but the accuracy of the angle setting, 723 , was of greater importance. The attachment of the ADV to the tow carriage was made using I-beams, C-beams and clamps. The ADV could be positioned with an accuracy of 72 cm using this setup. The ADV used is a downlooking 10 MHz ADVField Probe with optional sensors, connected to an ADVField processor, all from SonTek/YSI, Inc., with datalogging on a computer on the tow carriage. The accuracy for the ADV at the carriage speeds used, obtained from the supplier’s specs, is given in Table 1. The speed of the tow carriage was measured using a set of light gates that measured the travel time of the carriage over 1:395 m 7 0:002 m. The time measurement was made with 1 ms (millisecond) accuracy. The speed measured by the ADV was calibrated against the carriage speed measured by the light gates without the net panel in place. After calibration the error was well within the accuracy from the supplier (Table 1).

l

For the net used, solidity S was estimated to be 0.184. The solidity was also found using a digital photographic method where the pixels of the background color were sorted from the pixels of the net color. Using this technique, the solidity (S) was estimated to be 0.20 when the net was stretched on the frame. The measurements were performed at the tow speeds 12.5, 25, 50 and 75 cm s1 and at the angles of attack 03 , 153 , 303 , 453 , 603 and 903 . In addition, one test was done at 50 cm s1 and 753 . 3.2. Data processing Drag and lift forces acting on the empty frame and the frame with the net were found separately using loadcell data. One value for each of the forces was found for each run as an average of the steady state portion of the force time series. All drag and lift force measurements were the average of at least 120 measurements in each run and at least two tank runs were done for the empty frame. The drag force (D) and lift force (L) on the net were found by subtracting the averaged measurements for each speed setting and angle of attack on the frame from the averaged measurements for the corresponding speed setting and angle of attack on the

ARTICLE IN PRESS Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

frame and net. The drag coefficient (Cd ) and lift coefficient (Cl ) of the net panel at each speed and angle of attack were calculated as D 1 2

ð8Þ

rA/u0 S2

and L

Cl ¼

1 2

ð9Þ

rA/u0 S2

where A is the outline area of the net panel (1 m2 ) and /u0 S is the average tow speed measured by the light gates for all measurements at the same speed setting and angle of attack. The drag and lift coefficients were also calculated using projected thread area (Aproj ¼ AS) instead of outline area of the net panel. They were named Cd;proj and Cl;proj . The current reduction factor (Ur ) in chosen locations was found by Ur ¼

u0  /uS u0

ð10Þ

where u0 is the carriage speed measured using the light gates and /uS is the time series average of the water speed relative to the carriage measured by the ADV with the frame and net panel in place. 3.3. Results of experiments and discussion A plot of drag and lift coefficients (Cd and Cl ) as a function of angle of attack (a) for the different velocities used is shown in Fig. 2. The plot shows a clear trend for both Cd and Cl with a and that, for the most part, Cd and Cl decreased with increasing speed, which indicated a slight dependence on Reynolds number (Re). The variation of Cd and Cl with Re (decreasing Cd and Cl with increasing Re) was small, however, compared to the variation with a. The lift coefficient was non-zero (negative) at a ¼ 903 , and although this seemed strange, it could have been caused by some asymmetric feature of the net itself. During the tows, a vibration of the frame was observed at certain speeds. Namely at the lowest speed and the highest speed. The vibrations at the low speed were induced by some unsteadiness in the motion of the carriage while the vibrations at the highest speed were induced by the vortex shedding from the vertical frame posts and the frame supports. The vibrations resulted in fluctuating force measurements, but since the data were presented as average values of long time series, the accuracy of the data was assumed appropriate for the present study. This was verified by inspecting the standard error

1.5 Cd (12.5cm/s) Cd (25cm/s) Cd (50cm/s) Cd (75cm/s) Cl (12.5cm/s) Cl (25cm/s) Cl (50cm/s) Cl (75cm/s)

0.25

Cd, Cl

0.2 0.15

1 0.75

0.1

0.5

0.05

0.25 0

0 0

15

s

SEM ¼ pffiffiffi n

ð11Þ

where s is the standard deviation and n is the number of samples in the averaged section of the time series used for the force measurement. The SEM, for most of the measurements, was less than 10% and 15% of the corresponding D and L on the net. However, there was a small portion of the measurements where the SEM was larger. These were the lift measurements at 75 cm s1 for the empty frame, where the SEM was around 25% of the corresponding L on the net, and also the measurements where the force on the net was very small; D and L at a ¼ 03 and L at a ¼ 903 . The SEM was calculated for the separate measurements, but for each of the data points some replicate measurements were included for representing the D and L on the net. For the measurements on the empty frame, the average of two replicates was used for each data point. This improved the accuracy of the measurements. An inspection of the replicated measurements agreed with the quantification of the error by the SEM. Overall the measurements at u0 ¼ 50 cm s1 seemed to have the highest accuracy (SEM less than 2% for D and less than 8% for L), and the drag data seemed to be more accurate than the lift data. A plot of Ur as a function of a for the different tow speeds is shown in Fig. 3. The plot shows that for most of the tow speeds, the velocity reduction factor, Ur , behind the center of the net panel was of similar magnitude from a ¼ 603 to 903 , and was larger for lower a, but for a ¼ 03, Ur was smaller than for a ¼ 153. This could be due to the fact that for a ¼ 03 the main contribution to Ur was the front vertical pipe of the frame, while the rest (net panel and rear vertical pipe) was hidden in the wake of the front pipe and affecting the flow less than at a small, but larger than zero, angle of attack. The variation with speed was quite small compared to the associated measurement error, but it is interesting to see that Ur was generally smaller at 12:5 cm s1 than at 25 cm s1 , while in Fig. 2 it can be seen that the drag coefficient was generally larger at 12:5 cm s1 than at 25 cm s1 . This might have to do with increased mixing in the wake at slower speeds. It is also interesting to note that the low drag coefficient at u ¼ 12:5 cm s1 and a ¼ 903 is reflected in a low Ur . Fig. 3 also shows a few measurements of Ur as a function of angle of attack and distance across the tank from the centerline 2.5 m behind the net panel. Here it can be seen that the wake from the frame posts reached at least 2.5 m behind the frame. It can also be seen that there was a slight increase in velocity (negative Ur ) outside of the wake.

1.25 Cd,proj, Cl,proj

0.3

of the mean (SEM) of the measurements

30

45 α (°)

60

75

90

Fig. 2. Drag and lift coefficients as a function of angle of attack based on outline area (Cd and Cl ) and based on projected thread area (Cd;proj and Cl;proj ). The Reynolds number ðRe ¼ rud=mÞ based on thread diameter for the different speeds is: Reð25 cm s1 Þ ¼ 609, Reð50 cm s1 Þ ¼ 1218 and Reð12:5 cm s1 Þ ¼ 304, Reð75 cm s1 Þ ¼ 1826.

0.25

a

0.2 Ur

Cd ¼

317

u = 12.5cm/s u = 25cm/s u = 50cm/s u = 75cm/s

0.2 0.15

0.15

0.1

0.1

0.05

0.05

0

0

−0.05 0

15 30 45 60 75 90 α (°)

b

0

α = 90° α = 60° α = 30°

0.2 0.4 0.6 0.8 Distance (m)

1

Fig. 3. Ur 2.5 m behind the net panel; as a function of angle of attack and carriage speed, centered behind the net panel (a), and as a function of angle of attack and distance across the tank from the centerline (b). The dashed lines show the position of the vertical frame post for the different angles of attack. All measurements in (b) were performed at 0:5 m s1 .

ARTICLE IN PRESS 318

Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

4. Sensitivity study for CFD simulations Preliminary CFD runs were made to decide which settings, parameter choices and mesh to use when applying the model. The tank, net panel and frame were always modeled using their exact outer dimensions. In order for the problem to be steady state, the net panel was fixed, and relative velocity was achieved by specifying incident flow equal to the corresponding carriage speed. In general, it was found that, as long as guidelines based on physical consideration were adhered to, the model was remarkably insensitive to specific choices. The thickness of the porous media representing the net, for example, should be comparable to or larger than the thickness of the net, but much smaller than the size of the net panel. Trials indicated that, as long as the magnitudes of the resistance coefficients were inversely proportional to thickness, the choice of thickness did not affect the simulation results significantly, and a thickness of 50 mm was chosen. Similarly, trials indicated no significant dependence on options for near wall treatments, and a no shear stress condition was applied to the free surface, while the bottom and side-walls moved with the incident velocity and a no slip boundary condition was applied to the frame. The differences between the standard k2e model and the realizable k2e model for turbulence closures were small, and the realizable k2e model was selected for its slightly better stability and convergence characteristics for this specific problem. With respect to inlet turbulence, considerations of possible residual turbulence as well as spreading of the wake region in back of the frame members were taken into account, and compromise values of k ¼ 3:75  105 m2 s2 and e ¼ 2:5  107 m2 s3 were chosen. These values correspond to the RMS 0 ¼ 0:005 m s1 and value of the turbulent velocity fluctuations urms a turbulent length scale, describing the size of the energy carrying turbulent eddies, l ¼ 0:15 m. Decisions regarding the mesh were governed by limiting the total number of cells while still being able to resolve the frame member wakes. This was achieved utilizing FLUENT’s Hex Core meshing scheme in such a way that small tetrahedral cells described the region next to the frame around the net panel, small hexahedral cells described the region of the wake behind the frame, larger hexahedral cells described the rest of the interior of the domain and larger tetrahedral cells described the region next to the outer boundaries of the domain. Details of these and related preliminary test runs are provided by Patursson (2008).

Defining the xi0 system as the global coordinate system rotated around the vertical x3 - axis with x10 in the flow direction and x20 horizontal and perpendicular to the flow direction, Eq. (12) can be expressed as @p 1 ¼ Si0 ¼  Dij0 muj0  Cij0 rumag uj0 @xi0 2

ð13Þ

where Dij0 ¼ Rip Rjq Dpq

ð14Þ

and Cij0 ¼ Rip Rjq Cpq where 2

sina 6 R ¼ 4 cosa

ð15Þ

cosa sina

0

0

3 0 7 05

ð16Þ

1

and a is the angle of attack. Within the porous media, the pressure gradient is approximated by @p ðDpÞi  Dxi0 @xi0

ð17Þ

where ðDpÞi =Dxi0 (i not summed) is the pressure gradient across the thickness of the porous media Dxi0 in the xi0 direction. Pressure gradients through the porous media were assumed constant through the thickness though pressure on each face could vary. Velocity through the porous media was also assumed evenly distributed spatially across the net panel. Using a porous media with a thickness much smaller than the other dimensions, the drag, D, and lift force, L, were calculated using Eq. (13) as D ¼ ðDpÞ1 A10 ¼ S10 Dx10 A10 ;

L ¼ S20 Dx20 A20 :

ð18Þ

In Eq. (18), Ai0 is the outline projected area of the net panel in the xi0 direction (Fig. 4), A10 ¼ Asina;

A20 ¼ Acosa;

ð19Þ

Dxi0 is the distance through the net panel in the xi0 direction t ; sina

Dx10 ¼

Dx20 ¼

t ; cosa

ð20Þ

t is thickness, and a is the angle of attack. The drag and lift forces were calculated as D ¼ S10 tA;

L ¼ S20 tA

ð21Þ

and the drag and lift coefficients, Cd and Cl , were calculated as 5. Porous media resistance coefficients

Cd ¼

The porous resistance coefficients corresponding to the net used in Section 3 (Dn , Dt and Cn , Ct ) were obtained by minimizing the error between porous resistance model predictions and net panel drag and lift measurements. Because it would have been very time-consuming to use the full CFD computer model to generate multiple predictions, a simplified analytical model was used. In the analytical model, drag and lift were linearly related to the porous resistance coefficients, so standard analytical error minimization techniques could be employed straight forwardly.

2D Au20

r

¼

2S10 t ; ru20

Cl ¼

2S20 t : ru20

Note that D and L are linearly related to the porous media resistance coefficients, as well as velocity and velocity squared, through Eqs. (13) and (21).

5.1. Analytical approach for calculating forces on the net panel Combining Eqs. (1) and (4) and neglecting acceleration, gravitation and fluid stress while passing through the porous media, @p 1 ¼ Si ¼  Dij muj  Cij rumag uj @xi 2

ð22Þ

ð12Þ Fig. 4. Porous media sheet cross-section geometry.

ARTICLE IN PRESS Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

An approximation for ui0 , the velocity through the porous media in the xi0 direction, was needed. It was assumed that 0 1 uo ð1  rn Cd Þ B C 0 u~0 ¼ @ ð23Þ A 0 where Cd was found from measurements and rn was an estimate based on CFD simulations of flow through the porous media for a net panel at chosen angles of attack and different speeds. The expression for the velocity reduction was based on the formulation by Løland (1991). It was assumed that the across flow effect from the lift force was negligible. The parameter rn was estimated from previous CFD results and Eq. (23). The parameter rn changed quite a lot with angle of attack, but did not seem to depend much on speed, and the same value was used for all speeds. The parameter rn did not seem to be highly affected by the porous media resistance coefficients used, although the range of coefficients tested was limited since all were aimed at representing the net panel tested in Section 3. Fig. 5 provides calculated values of rn . Note that the assumptions underlying the analytical model development are best satisfied when a is not small. 5.2. Finding the best porous media resistance coefficients Three different error functions were minimized in order to identify the best way to find the porous media resistance coefficients Dn , Dt , Cn and Ct . The error functions used were a least squared normalized error (LSNE)   N  M  1X D  Dmeasured 2 1 X L  Lmeasured 2 þ ð24Þ LSNE ¼ D L N M a least absolute error (LAE) similar to that used by Zhan et al. (2006) P P jL  Lmeasured j jD  Dmeasured j þ P P LAE ¼ ð25Þ jDj þ jLj and a least absolute normalized error (LANE)   N  M  X   1X D  Dmeasured  þ 1 L  Lmeasured  LANE ¼     D L N M

ð26Þ

In Eqs. (24)–(26), D and L are evaluated from Eqs. (13) and (21), therefore proportional to porous media resistance coefficients, while subscript ‘‘measured’’ indicates values from the tow experiments.

Fig. 5. The variation of rn with angle of attack. The figure is based on CFD simulations using porous coefficients Dn ¼ 80 760 m2 , Dt ¼ 56 230 m2 , Cn ¼ 5:356 m1 and Ct ¼ 1:616 m1 performed at two speeds, u0 ¼ 0:125 and 0:5 m s1 .

319

A script was written in the Matlab software from The MathWorks, Inc. and used to solve the problem of minimizing the error functions for the variables Dn , Dt , Cn and Ct utilizing a minimizing function, Eq. (21) and the above mentioned assumptions. The data used for the fitting were tow tank data explained in Section 3, where drag, Dmeasured , and lift force, Lmeasured , were measured as a function of angle of attack, a (0, 15, 30, 45, 60, 75 and 903 ) and tow speed, u0 (12.5, 25, 50 and 75 cm s1 ). A least squares fitting function was used to solve for LSNE and a function was written to solve for the absolute error functions. The resulting porous media resistance coefficients are given in Table 2. Cd , Cl and Ur were calculated for comparing to the measured values. Cd and Cl were calculated using Eq. (22) and the fitted coefficients in Table 2. Ur was calculated as Ur ¼ rC d

ð27Þ

where r was found in a similar manner as rn and was calculated as r¼

u0  u2:5 Cd u0

ð28Þ

where u2:5 is the current velocity 2.5 m behind the center of the net panel found from CFD simulations. Fig. 6 shows the r-values used. The resulting fit between calculated and measured Cd , Cl and Ur is shown in Fig. 7. The coefficients found using the LAE error function seemed to overpredict the velocity dependence of the drag coefficient and, hence, also of the velocity reduction, and the lift coefficient seemed to have a reversed velocity dependence. The coefficients found using the LSNE and the LANE error functions showed quite similar results, but since the LSNE error function depended on the squared error, this error function had a higher response to questionable data points. This was the reason for the large velocity dependence of the lift coefficient observed in Fig. 7(a). If the lift force measured at 12:5 cm s1 and 603 and at 75 cm s1 and 453 were omitted, the fit using LSNE and LANE would be almost the same. Therefore, it was decided to use the coefficients from the LANE error function for further analysis. Table 2 Porous media resistance coefficients from the different error functions. Error function

Dn ðm2 Þ

Dt ðm2 Þ

Cn ðm1 Þ

Ct ðm1 Þ

LSNE LAE LANE

75 854 76 486 51730

35 409 84 741 26 379

4.8419 5.0727 5.0980

1.4439 1.1300 1.6984

Fig. 6. The variation of r with angle of attack. The figure is based on the same data as Fig. 5.

ARTICLE IN PRESS 320

Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

0.3

0.5

0.1

12.5 cms−1 25 cms−1 50 cms−1 75 cms−1

0.4 0.05

UR

Cl

Cd

0.2

0.3 0.2

0.1

0.1

0 0

0 0

30

60

90

0

30

α (°)

60

90

0

30

α (°)

12.5 cms−1 25 cms−1 50 cms−1 75 cms−1

0.4 0.05

UR

Cl

Cd

0.2

90

0.5

0.1

0.3

60 α (°)

0.3 0.2

0.1

0.1

0

0

0 0

30

60

90

0

30

α (°)

60

90

0

30

α (°)

90

60

90

0.5

0.1

0.3

60 α (°)

0.4 0.05

UR

Cl

Cd

0.2

0.3 0.2

0.1

0.1

0

0

0 0

30

60

90

0

30

α (°)

60

90

0

30

α (°)

α (°)

Fig. 7. Data from the fits using LSNE (a), LAE (b) and LANE (c). The lines correspond to the optimized analytical models; the symbols to measurements.

S = 0.243

S = 0.317 1

0.2

0.3

0.75

0.15 0.1

Cd

0.4

Cd

Cd

S = 0.13 0.25

0.2 0.1

0.05 30

60

90

0.5 0.25

0 0

0.159 ms−1 0.316 ms−1 0.966 ms−1

0 0

30

α (°)

60

90

0

0.159 ms−1 0.316 ms−1 0.966 ms−1

0.3 Cl

Cl

Cl

0.1 0.05

0.02 0 60 α (°)

90

0.2 0.1

0 30

90

0.4

0.15

0.04

60 α (°)

0.06

0

30

α (°)

0 0

30

60 α (°)

90

0

30

60

90

α (°)

Fig. 8. Data from the fit of the data from Rudi et al. (1988) using LANE. The lines refer to the optimized analytical model; the symbols show measurements.

ARTICLE IN PRESS Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

321

Table 3 Porous coefficients for the different data sets from Rudi et al. (1988). S

Dn ðm2 Þ

Dt ðm2 Þ

Cn ðm1 Þ

Ct ðm1 Þ

0.130 0.243 0.317

112 250 152 760 845 370

34 352 37 283 105 900

3.3015 6.1493 12.401

1.9928 2.8146 6.5737

To examine the robustness of the routine for finding porous coefficients, it was tested on data measured by Rudi et al. (1988). The reported data included drag and lift coefficients for several net panels at different angles of attack. The measured and calculated drag and lift coefficients for three net panels of different solidities are shown in Fig. 8, and the associated porous coefficients are given in Table 3. The measured value for drag at a ¼ 03 seemed to be overpredicted, at least for some of the measurements, so it was omitted.

6. Results from computational fluid dynamics CFD simulations were run modeling the net panel and frame at six angles of attack and two different speeds. The results were then compared to the tow tank measurements. The speeds 12.5 and 50 cm s1 , and angles of attack, 0, 15, 30, 45 and 903 , were chosen to include the most reliable tow tank measurements. The mesh was created using the Hex Core meshing scheme with 10 mm cell size next to the frame and in the wake of the frame. The growth factors are 1.2 from the frame and 1.15 in the wake. The standard wall functions were used. All walls were moving walls except for the surface that was modeled as a wall with zero shear force. The turbulence model used was the realizable k2e model, and the incoming turbulence k ¼ 3:75  105 m2 s2 and e ¼ 2:5  107 m2 s3 was based on the RMS value 0 ¼ 0:005 m s1 and a of the turbulent velocity fluctuations urms turbulent length scale, describing the size of the energy carrying turbulent eddies, l ¼ 0:15 m. The porous coefficients used were obtained from the LANE error function technique, and are provided in Table 2. Contour plots of velocity on a horizontal cut through the center of the net panel at various angles of attack using a current speed of 50 cm s1 are shown in Fig. 9. It is interesting to note the structure of the wake behind the net. The wake continues for a long distance, many times the width of the net panel, behind the net panel in the flow direction. In the across flow direction the dimensions of the wake are approximately the same as those of the net panel. The reduction of the current speed is a continuous process that starts a short distance in front of the net and continues through the net panel and some distance behind the net panel until it reaches a level about 11% less than the incoming flow. From this point the wake is very steady. The distance from the net panel to the start point of this steady level is about one to two times the width of the net panel. This means that the chosen point of measuring current in Section 3 for finding the velocity reduction is within this region of steady reduced flow. The variables Cd , Cl and Ur were compared as shown in Fig. 10, and the difference between measured and modeled data was generally small with a few exceptions.

 The difference between measured and modeled Ur was quite large at small angles of attack. This was mainly due to the fact that the wake from the frame around the net panel was narrower in the model than what was measured, and the high Ur at a ¼ 153 in the measured data was due to the wake from the frame, while the effect from the frame was smaller or

Fig. 9. The plots show velocity contours obtained from the CFD simulations on a horizontal cut through the center of the net panel at different angles of attack. The angles of attack used are from top to bottom 903 , 453 , 303 and 03 . The current speed used is 50 cm s1 .

 

nonexistent in the CFD prediction. There was larger diffusion of the wake behind the model at slower speeds, and this was also reflected in the comparison with a smaller difference between measured and modeled data at u0 ¼ 0:125 m s1 . The measured data point of Cl at u0 ¼ 0:125 m s1 and a ¼ 603 looks suspicious and might be a bad data point. The measured Cl was less than zero at a ¼ 903 . This might have been a correct measurement due to some anti-symmetries in the net panel, but the model did not have the capability to predict a non-zero Cl at a ¼ 903 .

6.1. Accuracy of the porous resistance coefficient procedure To test whether the procedure for finding the porous coefficients was actually minimizing the error of the CFD simulations at the chosen values, the effect of varying the porous resistance coefficients about the value found using the method in Section 5.2 was explored. This was done by calculating the same error function as used for the fit in Section 5.2, LANE, between the

ARTICLE IN PRESS 322

Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

u0 = 0.125ms−1 CFD Tank

0.08 0.06 Ur

0.2

0.04

Cl

Cd

0.2 0.1

0.02

0

−0.02

0.1

0 0

30

60

90

0 0

30

α (°)

60

90

0

30

α (°)

60

90

60

90

α (°)

u0 = 0.5ms−1 0.08 0.06 Ur

0.2

0.04

Cl

Cd

0.2

0.02

0.1

0.1

0 0

0

−0.02 0

30

60

90

0

30

α (°)

60

90

0

30

α (°)

α (°)

Fig. 10. Comparing the results from the CFD simulations to measurement data.

All data points

Limited number of data points

0.15

0.1

0.14

0.08 LANE

LANE

0.13 0.12 0.11

0.06 0.04 0.02

0.1 0.09 −10

Dn Dt Cn Ct

−5

0 Offset (%)

5

10

0 −10

−5

0 Offset (%)

5

10

Fig. 11. Sensitivity of analytical method to offset in porous resistance coefficients. Each line represents the variation of the error function, LANE, as a function of the offset of the respective coefficient, while the other coefficients are not changed. The figure on the left was made using all data points in the measurements, while the figure on the right only included four data points, u0 ¼ 0:125 and 0:5 m s1 at a ¼ 453 and 903 .

CFD results and the measured data. Due to the long simulation time required, it was chosen not to run CFD simulations for all the data points in the measurement series. A reduced number of data points were chosen that showed a similar error function as the whole data set. The choice of data points was examined using the analytical model (Fig. 11) where the porous resistance coefficients were offset a percentage of the fitted coefficients, found using the LANE error function, given in Table 2. The offset was done for one coefficient at a time while the other three coefficients were held fixed at the fitted value. Fig. 11 shows that the analytical error function using only four data points, u0 ¼ 0:125 and 0:5 m s1 at a ¼ 453 and 903 , had similar behavior to the error function using all the data points in the measurement series. The error itself was smaller due to the fact that these measurements generally had smaller errors than many of the other data points. The accuracy of the porous resistance coefficient procedure was tested by running the model with the porous resistance coefficients varied around the predicted value and calculating the

error function similar to what was done for Fig. 11, but for the limited number of data points. The model was run for the four scenarios u0 ¼ 0:125 and 0:5 m s1 at a ¼ 453 and 903 for each of the sets of porous resistance coefficients. For each of the sets of porous resistance coefficients, the LANE error function was calculated in addition to another error function (LANEu ) where Ur was included LANEu ¼    K  M  N  X X    1X D  Dmeasured  þ 1 L  Lmeasured  þ 1 Ur  Ur;measured   M  N     D L Ur K ð29Þ The results are shown in Fig. 12. The results show that the coefficients with the largest effect on the error were chosen well. Cn which had the largest effect on the calculated error was chosen within a few percent, while Ct which

ARTICLE IN PRESS Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

Original error function

323

Including current reduction

0.08 Dn Dt Cn Ct

0.08 LANEu

LANE

0.06

0.1

0.04

0.06 0.04

0.02

0.02

0

0 −20

0 20 Offset (%)

40

−20

0 20 Offset (%)

40

Fig. 12. Sensitivity of CFD method to offset in porous resistance coefficients. Each line represents the variation of the error function, LANE, as a function of the offset of the respective coefficient while the other coefficients were not changed. The figures were made using the four data points, u0 ¼ 0:125 and 0:5 m s1 at a ¼ 453 and 903 . The figure on the left shows the original LANE error function (Eq. (26)), and the figure on the right plots the LANEu error function where the current reduction was included (Eq. (29)).

also had a fairly large effect on the error was underpredicted by around 10%. The other two coefficients were not predicted very accurately, but they had a very limited effect on the error, especially Dt .

7. Discussion The porous media model was applied to the 1 m by 1 m square net within the pipe frame subject to incident flow at various speeds and angles of attack. The porous media resistance coefficients (Cn , Ct , Dn , Dt ) were obtained such that a best fit was achieved between calculated and measured lift and drag coefficients. To expedite this process, an analytical model was used which approximated local flow conditions through the net at the incident angle. Three error minimization criteria were applied, compared and evaluated. This process indicated that the least absolute normalized error (LANE) method provided the most desirable weighting of error contributions, and the coefficients obtained using LANE should be the ones used in subsequent CFD applications. The analytical model with LANE minimization procedure was also applied to Rudi et al. (1988) net panel lift and drag measurements thereby generating resistance coefficients for three additional net types. Employing the porous resistance coefficients obtained using LANE for the net used in this study (Table 2), the full, threedimensional, CFD porous media model was applied to the net panel (net plus solid pipe frame) used in the tow tank experiments. Computed and measured lift and drag coefficients, as well as velocity reduction in the wake (not used to obtain the resistance coefficients), were compared to provide an initial validation of the porous media CFD methodology. The sensitivity of the results to the resistance coefficients chosen was also explored by offsetting resistance coefficient values and plotting LANE as a function of resistance coefficient variation. The porous media method of modeling the flow through a net panel gave CFD predictions which agreed well with measured data for both drag and lift forces on the net panel (LANE ¼ 2:4%) and current reduction behind the net panel (LANEu ¼ 6:2%). Also the proposed analytical method of finding the porous resistance coefficients of the net gave coefficients for which model results agreed fairly well with the data (LANE ¼ 10%). The proposed mesh and setup of the CFD model gave good results overall. The wake of the frame around the net panel was not modeled well, however, and this resulted in a rather large discrepancy between predicted and measured current reduction behind the net panel at small

angles of attack when the contribution of the frame becomes significant. This was not considered a problem since the focus was on modeling the net and not the structure holding the net panel. Since the measured current reduction data were not used in the method for finding the porous resistance coefficients, the discrepancy between measured and modeled velocity reduction at small angle of attack did not affect the coefficients found, but the comparison between measured and modeled velocity reduction served as a first validation of the approach. This validation supports future application to full net enclosures since the difference was very small except for the mentioned discrepancy due to the net’s rigid frame.

8. Conclusion This method should be very useful for modeling the flow through and around structures of net. Application of the method, however, requires the resistance coefficients for each individual net that is modeled. Finding the coefficients is not straight forward, and the method used here requires values for measured drag and lift forces as a function of speed and angle of attack. Future research should include validating the approach for large structures of net such as fish farming cages and working towards easier methods for finding the porous media resistance coefficients.

Acknowledgments The present work is funded in part by the Statoil Group, which is a consortium of oil companies, including Statoil, Phillips, Enterprise, and Veba who are involved in oil exploration on the shelf of the Faroe Islands, and in part by the Cooperative Institute for New England Mariculture and Fisheries (CINEMAR) and the National Oceanic and Atmospheric Administration (NOAA). References D’arcy, H., 1856. Les fontaines publiques de la ville de dijon. Dalmont, Paris. Ferziger, J.H., Peric´, M., 2002. Computational Methods for Fluid Dynamics, third ed. Springer, Berlin. Fluent, 2006. FLUENT 6.3 User’s Guide. Fluent Inc., Lebanon, NH, USA, September. Fluent, 2007. GAMBIT 2.4 Modeling Guide. Fluent Inc., Lebanon, NH, USA, May. Forchheimer, P., 1901. Wasserbewegung durch boden. Zeitschrift des Vereines Deutscher Ingenieure 45, 1782–1788. Fredheim, A., 2005. Current forces on net structures. Ph.D. Thesis, NTNU Department of Marine Technology, Faculty of Engineering Science and

ARTICLE IN PRESS 324

Ø. Patursson et al. / Ocean Engineering 37 (2010) 314–324

Technology. Trondheim, Norwegian University of Science and Technology, Doctoral Theses at NTNU, 2005:64. Gui, F.-k., Li, Y.-c., Zhao, Y.-p., Dong, G.-h., 2006. A model for the calculation of velocity reduction behind a plane fishing net. China Ocean Engineering 20 (4), 615–622. Helsley, C.E., Kim, J.W., 2005. Mixing downstream of a submerged fish cage: a numerical study. IEEE Journal of Oceanic Engineering 30 (1), 12–19. Launder, B.E., Spalding, D.B., 1972. Lectures in Mathematical Models of Turbulence. Academic Press, London, England. Le Bris, F., Marichal, D., 1998. Numerical and experimental study of submerged supply nets: applications to fish farms. Journal of Marine Science and Technology 3, 161–170. Løland, G., 1991. Current forces on and flow through fish farms. Ph.D. Thesis, University of Trondheim. O’Neill, F.G., 2006. Source models of flow through and around screens and gauzes. Ocean Engineering 33, 1884–1895. Patursson, Ø., 2007. Tow tank measurements of drag- and lift force on a net panel and current reduction behind the net panel. Technical Report NVDRit2007:10, University of the Faroe Islands, To´rshavn, Faroe Islands.

Patursson, Ø., 2008. Flow through and around fish farming nets. Ph.D. Dissertation, Ocean Engineering, University of New Hampshire, Durham, NH, USA. Patursson, Ø., Swift, M.R., Tsukrov, I., Baldwin, K., Simonsen, K., 2006. Modeling flow through and around a net panel using computational fluid dynamics. In: OCEANS’06 MTS/IEEE - Boston, Boston, MA, USA. Rudi, H., Løland, G., Furunes, I., 1988. Experiments with nets; forces on and flow trough net panels and cage systems. Technical Report MT51 F88-0215, MARINTEK, Trondheim, Norway (in Norwegian). Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k–e eddy viscosity model for high Reynolds number turbulent flows. Computers and Fluids 24 (3), 227–238 (Introduction of the realizeable k–epsilon model). Vincent, B., Marichal, D., 1996. Computation of the flow in the codend. In: Fishing Technology and Fish Behavior Working Group. International Council for the Exploration of the Sea. Zhan, J.M., Jia, X.P., Li, Y.S., Sun, M.G., Guo, G.X., Hu, Y.Z., 2006. Analytical and experimental investigation of drag on nets of fish cages. Aquacultural Engineering 35, 91–101.