Abstract
Hydrodynamic effects of mussel farms have attracted increased research attentions in recent years. The understanding of the hydrodynamic impacts is essential for predicting the sustainability of mussel farms. A large mussel farm includes thousands of mussel droppers, and the combined drag on the mussel droppers is sufficient to possibly affect the longevity of the entire long-lines. This article intends to study the drag and wake of an individual long-line mussel dropper using computational fluid dynamics approaches. Two equivalent rough cylinders, namely, Curved-Model and Sharp-Model, have been utilized to simulate the mussel dropper, and each rough cylinder is assigned with surface roughness. The porosity is not considered in this article due to its complexity from inhalant and exhalant of mussels. Two-dimensional laminar simulations are conducted at Reynolds number from 10 to 200, and three-dimensional large eddy simulations are conducted at subcritical Reynolds number ranging from 3900 to
Introduction
According to the report from Food and Agriculture Organization (FAO), almost 109 molluscs species have been cultivated in aquaculture farms all over the world. 1 Among all these molluscs, mussel farming is the most popular aquaculture activity with a very rapid growth due to the characteristics of carbon sequestration and environmental friendliness. The suspended long-line culture is the most common way for growing mussel. Long-line structures are approximately 100 to 200 m in length and consist of two parallel backbone ropes. The crop ropes (mussel droppers) are hung vertically from the backbone ropes down to 10 to 15 m of the water column. The commonly used surface mussel long-line is depicted in Figure 1.

Schematics of the long-line: (a) elevation view, (b) plan view, (c) mussel dropper line arrangement, and (d) close-up view of the mussel dropper.
Large-scale, open-water mussel farming is a sustainable means of providing food to satisfy human demand. However, offshore harsh environment leads to challenges for structure implementation, namely, the movement of the long-lines and the supporting structures induced by the impacts of waves and current, and unpredictable crop losses. The interactions between high energy waves, tidal current, and the long-lines significantly affect the stability of the structures. 2
Previous works investigated the suspended mussel droppers as rough cylinders, and the empirical values of hydrodynamic coefficients (these are drag and inertia coefficients),
Most of the previous research works concentrated on the dynamic response and shielding effects of the mussel farm. The understanding of hydrodynamic coefficients of the mussel dropper still has not been completely settled down. Also, the effects of mussel shell on the flow pattern of the mussel dropper are barely documented. With these two concerns, this study intends to examine the impact of surface roughness of the mussel dropper on the drag and wake of an individual dropper. Two different equivalent cylinders with the same roughness of the mussel dropper have been tested under two CFD schemes.
It is noted that the inhalant and exhalant siphons allow mussels to exchange flow to acquire nutrients from the surroundings, causing the mussel dropper more likely to be porous instead of solid. Previous investigations showed that a porous cylinder had a drag coefficient 20% smaller than the solid cylinder, and the vortex shedding was also weaker behind the porous cylinder. 14 Therefore, to study the drag and wake of the mussel dropper, both roughness and porosity should be considered. However, we focus mostly on the effects of roughness on the drag and wake of the mussel dropper in this article, and thus, the assessment was conducted without considering the complexity of the porosity.
The remainder of the article is organized as follows. In the “Methodology” section, we first introduce the mathematical formulations and numerical discretization of the CFD approaches, which include a two-dimensional (2D) laminar simulation and a three-dimensional (3D) large eddy simulation (LES). The results of the simulations are presented in section “Results and discussion.” The conclusions are given in section “Conclusion.”
Methodology
Laminar simulation
Two-dimensional simulations are conducted for laminar flow around smooth and rough cylinders at Reynolds numbers
We assume the flow is incompressible. The governing Navier–Stokes equations for an incompressible 2D flow regarding the conservation of mass, momentum, and energy are defined as
where
where
The governing equations in this study are a set of partial differential equations (PDEs). We utilize the finite volume method (FVM) with structured grids to discretize these PDEs. The computational domain is represented by numerical grids at which the variables can be transferred and calculated. Structured grids include three basic topologies: H-, O-, and C-grids. All three topologies have been employed for investigating smooth cylinder. The names of the grid topology are derived from the shapes of the grid lines. 16
The discretization is to convert the PDEs into non-linear algebraic equations. For unsteady flows, an elliptic problem has to be solved at each time step. Steady flows are solved by an equivalent iteration scheme, and then problems turn out to be solutions of linear equation systems. Some convergence criteria need to be checked after all iterations are finished. The repeat of the loop depends on the condition if the convergence is satisfied. We utilized commercial code ANSYS ICEM to create three different computational domains for smooth circular cylinder as depicted in Figure 2.

Two-dimensional computational domains of flow around smooth circular cylinder and boundary conditions: H-, O-, and C-grids.
Different grid topologies used in this study enable validation for grid independence study under the same algorithm scheme. The most optimal grid topology can be identified after comparison, to meet the accuracy requirement and the computational efficiency. For these three different grids, the boundary conditions are set as follows:
H-grid: we set left side as flow inlet
O-grid: left and right sides are set as flow inlet
C-grid: left, top, and bottom sides are set as flow inlet
Only if stated, otherwise, all positions and length scales are normalized by the characteristic length
The coupled scheme is implemented for pressure–velocity coupling where convective term discretization is conducted by semi-implicit method for pressure-linked equations (SIMPLE) schemes. 17 The spatial discretization of the pressure and momentum is conducted through second-order and second-order upwind differencing schemes, respectively. The temporal discretization is conducted using second-order implicit differencing scheme. The algebraic equations are solved by the Gauss–Seidel iterative method in conjunction with the Algebraic Multigrid (AMG) solver. The AMG method greatly reduces the number of required iterations (and thus, CPU time) to obtain a converged solution, particularly when the model contains large number of control volumes. The time step varies from 0.01 to 0.05 for steady and unsteady flow conditions, respectively. The convergence criteria for the inner iterations (time steps) are set as 10−9 for the discretized continuity, momentum, and energy equations.
Mussels form dense aggregations which grow on suspended aquaculture long-lines. The diameter of the mussel droppers varies as mussels mature. For example, the typical averaged diameter of the mussel dropper approaches approximately 0.15–0.2 m, with the typical peak tidal velocities of the mussel farm ranging from 0.05 to 0.2 m/s, Reynolds numbers of a mussel dropper are between

Two-dimensional equivalent rough cylinder models: (a) real mussel dropper, (b) geometries of equivalent rough cylinders, and (c) grids of equivalent rough cylinders (picture of real mussel dropper is reproduced from Plew 2 ).
Based on the geometries of the equivalent rough cylinders, we define them as Curved-Model and Sharp-Model. By adding the crowns in the geometries of the rough cylinders, mussel shells can be represented without losing the geometrical characteristics. Surface roughness,
LES
Three-dimensional LES are conducted at Reynolds number
where
To ensure the grid orthogonality,
18
the computational domain for 3D simulations utilizes an O-grid, which is extruded by a 2D O-grid in ANSYS ICEM. The diameter of the 3D O-grid is 40D, where

Three-dimensional computational domains of flow around smooth and rough cylinders.
We describe the grids of the 3D computational domain by polar coordinates. Cell numbers are presented by
The first layer of the cell is set as
Inflow boundary: the left portion of the cylindrical domain is set as velocity inlet, where the flow is uniformly distributed with the turbulence intensity of 1%.
Outflow boundary: the pressure outlet is imposed at the right portion of the cylindrical domain. A Neumann boundary condition
Surface boundary: for the surface of the cylinder, the Dirichlet boundary condition is used for the velocities
Span-wise boundary: the periodic condition is implemented for the span-wise direction of the cylinder.
The pressure implicit splitting of operators (PISO) scheme is utilized for the velocity–pressure coupling. Pressure and momentum are discretized by second-order and bounded central differencing schemes, respectively. For the temporal discretization, the bounded second-order implicit differencing scheme is used. Time step is set as
Numerical accuracy
To evaluate the grid sensitivity and numerical accuracy, different computational domains and cell numbers have been examined for 2D laminar and 3D LESs, respectively. The descriptions of different computational domains (H-, O-, and C-grids) for smooth cylinders are illustrated in Table 1. The comparison of different grids is for one thing to conduct the validation, and for another, to select the most optimal grid topology for further rough cylinder simulations.
Different computational domains for 2D laminar flow around smooth cylinder.
Previous numerical simulation
15
found that the obtained Strouhal number

Comparisons of drag coefficient

Comparisons of Strouhal number
Figure 5 shows that the overall dependences of drag coefficients on Reynolds number are as expected, and the effect of domain dimension on drag coefficient is insignificant. Since all the 2D simulations are conducted in the same algorithm scheme, if a great discrepancy occurs, there must be some unexpected errors need to be checked from the beginning of the loop. For the H-grid under steady flow conditions
There are obvious differences in the dependence of Strouhal number on Reynolds number (Figure 6). Small domain obtained large Strouhal number, which verified the previous numerical investigation. 15 The Strouhal numbers obtained from the medium and the large domains are very close to each other. We then select the medium domain with respect to the requirement of computational efficiency and numerical accuracy. With this medium domain, the drag coefficient and Strouhal number obtained in different grids are compared (Figure 7).

Comparisons of drag coefficient
Figure 7 shows that for the steady flow conditions, the drag coefficients obtained from H-grid are larger than that from others, while the drag coefficients obtained from C-grid are smaller than that from others. For the unsteady flow conditions, the results are almost same for these three different grids. The same trend of Strouhal number can be achieved from different grids under medium domain dimensions. The vorticities behind smooth cylinders in different grids are demonstrated in Figure 8.

Comparisons of vorticities behind smooth cylinders in different 2D grids.
The vorticities generated in C-grid appear quite sparse in far field. H- and O-grids predict similar vortex patterns. The comparisons of drag coefficients and Strouhal numbers show that the medium domain is compatible to meet the requirement of computational efficiency and numerical accuracy. Therefore, we utilize O-grid to conduct the simulations for rough cylinder without any further verification, and the domain dimension is set as 30D based on smooth cylinder cases. The cell numbers and 3D LES grid independence studies for smooth and rough cylinders at Reynolds number
Numerical accuracy studies of different cell numbers
The effect of the cell numbers on drag coefficient and Strouhal number is insignificant for the smooth cylinders. The present simulation results for smooth cylinder are in good agreement with previous work conducted by LES.
18
Therefore,
Results and discussion
Drag coefficient and Strouhal number
Two equivalent rough cylinders (Curved-Model and Sharp-Model) are utilized to model the mussel droppers. Each equivalent rough cylinder has been assigned surface roughness

Effects of surface roughness on drag coefficient and Strouhal number.

Effects of crown shapes on drag coefficient and Strouhal number.
Similar with the trend of smooth cylinder, drag coefficients of the equivalent rough cylinders decrease with the increase in Reynolds numbers, while Strouhal number increases with the increase in Reynolds numbers. In terms of effect of surface roughness, drag coefficient and Strouhal number increase with the surface roughness. In particular, Strouhal number ranges from 0.2 to 0.24 for rough cylinder in laminar flow. The considered range for smooth cylinder is from 0.1 to 0.2 when Reynolds number is below 200, demonstrating that the surface roughness has significant effect on vertex shedding of the rough cylinders.
Figure 10 shows that sharp crown has more effects on drag coefficient and Strouhal number. Drag coefficients and Strouhal numbers of Sharp-Model are greater than that of Curved-Model. As the surface roughness increases, there are oscillations in the Strouhal number, which indicates that sharp crowns induced more violent flow oscillations.
Differing from the 2D simulations, the flow oscillations become even more violent at subcritical Reynolds number in 3D LESs. To study the vertex shedding of the rough cylinders, Figure 11 exhibits the dependence of lift coefficient on non-dimensional time at

Dependence of drag and lift coefficients on non-dimensional time at
Figure 11 shows that the oscillation of lift coefficient is periodic because of the periodic boundary condition on the span-wise direction. Particularly, this periodic oscillation shows more regular pattern after 20 s. However, the wake oscillations of the Sharp-Model exhibit more chaotically at

Dependence of lift coefficients on non-dimensional time at
The oscillation of lift coefficient is partially induced by periodic boundary conditions. Besides, due to the blocking effects of the sharp crowns around the rough cylinder, more active oscillation occurs compared with the situations for Curved-Model. Lift coefficient at the initial stage has larger amplitude around 0.05, as the wake development, boundary layer shear flow transports from 2D to 3D, the averaged amplitude of lift coefficient reduces approximately 30% compared with that of the primary stage.
To evaluate the effect of Reynolds number on drag coefficient at all considered subcritical Reynolds number, Figure 13 presents the dependence of drag coefficients on Reynolds numbers for smooth and equivalent rough cylinders. Given the surface roughness

Dependence of drag coefficients on Reynolds numbers for smooth and equivalent rough cylinders
The overall trends of the drag coefficients for smooth and rough cylinders in this study are similar to Hoerner’s report.
22
The drag coefficient of the rough cylinder ranges from 1.1 to 1.2 at subcritical Reynolds numbers. Even with the same surface roughness, the drag coefficients show different values due to the crown shapes on the surface of the rough cylinders. The drag coefficients for the Sharp-Model are larger than that of others. When Reynolds number is less than O(
Wake and flow pattern
To observe the effects of surface roughness on the flow pattern, Figure 14 presents the 2D streamlines of Curved-Model and Sharp-Model at

Streamlines of flow around different equivalent rough cylinders at
When Reynolds number approaches 3900, the wake behind the rough cylinder is completely turbulent flow, but the boundary layer is still not fully developed as turbulent flow. The separation of the flow occurs when flow passes the rough cylinder, and then the flow rolls back to the wake to be absorbed onto the cylinder again, as shown in Figure 15(a).

Flow around Curved-Model at
Turbulent kinetic energy (TKE) varies at different wake positions, in particular, TKE approaches its peak in the wake roll-back position (position B) as shown in Figure 15(b). Position A in Figure 15(a) indicates that the wake is rolling back to form a secondary wake. TKE is defined in the following equation
where

Turbulent vorticity of Sharp-Model at
To study the effects of crowns around the rough cylinder on the wake in the near field, Figure 17 exhibits the turbulent oscillation velocities along the centerline of the cylinder. The comparison between smooth and rough cylinders is conducted for further analysis. In Figure 17,

Time-averaged stream-wise velocity along the wake centerline at
The overall trend of the velocity distribution is similar to the demonstration in Cheng et al.
18
The turbulent fluctuation velocity decreases along the wake centerline in the near-wake region to reach a minimum value and then increases in the far-wake region. The velocities for the Curved-Model and Sharp-Model are smaller compared to the case of smooth cylinder, due to the blockage of the crowns on the equivalent rough cylinders. Cell numbers utilized in Cheng et al.
18
are about 210 million with
The cell numbers could be the reason that makes the present simulation results slightly different with the results from Cheng et al.
18
For example, the valley value obtained from the present LES is smaller than previous simulation at

Turbulence kinetic energy (TKE) along the wake centerline at
Figure 18 demonstrates that the overall trends for smooth and rough cylinders are similar. The velocity is zero in the wake roll-back region from 0.1D to 0.4D. TKE increases as
Q-criterion
23
is utilized to present the iso-surfaces of vorticities in this study. The basic principle is to decompose the velocity gradient

Iso-surface of Q-criterion vorticities at

Iso-surface of Q-criterion vorticities at
Figure 19 shows the 3D vorticities behind the cylinders. When the cylinder varies from the smooth to the rough, the vorticity formulations change differently. The vorticity tubes behind the Sharp-Model indicate that the surface roughness, especially the crown shapes, significantly alters the vorticity formulations.
When Reynolds number approaches O(
Conclusion
In this article, CFD approaches have been utilized to investigate the drag and wake of a long-line mussel dropper. 2D laminar and 3D LESs were conducted at considered Reynolds numbers. The grid independence study and numerical accuracy were examined for convergence analysis. Two equivalent rough cylinders were used to model the mussel dropper, through which the effect of surface roughness on the drag and wake can be assessed. The conclusions are as follows:
For the smooth cylinder in 2D simulations, the comparisons between three different grids showed that O-grid with domain diameter of 30D could feasibly predict the drag and wake without losing accuracy. For the equivalent rough cylinders in 2D simulations, the drag coefficients increased with the surface roughness. The crowns around the rough cylinder significantly affected the drag coefficient. Sharp crowns resulted in larger drag coefficients than that of others. The streamlines for equivalent rough cylinders were very different from those of smooth cylinders.
The drag coefficient of the rough cylinder ranges from 1.1 to 1.2 at subcritical Reynolds numbers. Sharp crowns induced more energy loss and the wake is more fluctuated. As a result, the Strouhal number was in a range from 0.25 to 0.26. This number was much larger than that of smooth cylinders (around 0.21). The equivalent rough cylinder is an efficient way by which the drag and wake of a mussel dropper can be examined using CFD scheme. The porous characteristics induced by inhalant and exhalant of mussels should be emerged in further studies.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was financially supported by the National Natural Science Foundation of China (Grant Nos 51679046 and 51909040).
