Modelling Open Ocean Aquaculture Structures Using CFD and a Simulation-Based Screen Force Model

: The numerical framework of the open source CFD solver REEF3D is utilised to study the ﬂuid–structure interaction of an open ocean aquaculture system and waves. The presence of the net is considered in the momentum equations of the ﬂuid using a forcing term based on Lagrangian– Eulerian coupling and the hydrodynamic loads on the net. They are deﬁned semi-empirically using a screen force model. Here, the hydrodynamic force coefﬁcients are calculated from the net geometry and ﬂuid velocity. The necessary force coefﬁcients are predicted from a new simulation-based screen force model. Here, CFD simulations are performed to obtain the hydrodynamic loads on net panels for varying net geometries, angles of attack and velocities. Then, a Kriging metamodel is applied to ﬁt a polynomial to the data. The proposed net model is validated against measurements for waves and current through rigid net panels and applied to simulate the dynamics of an open ocean aquaculture cage in waves. In current, the model predicts the drag forces and velocity reduction within a 10% error band, whereas it tends to under-predict the lift forces by up to 20%. In waves, the model tends to over-predict the crests with increasing wave height, but the deviations are also within a 10% error band.


Introduction
The Open Ocean Aquaculture (OOA) industry has seen strong growth over recent years due to an increasing demand for aquaculture products. Compared to traditional designs, OOA structures are larger in size and encounter larger environmental loads in offshore conditions. Therefore, advanced tools are necessary for the design process to better understand the dominating factors for dangerous load scenarios. Separated studies on the deformation of cages and wave loading, respectively, have been conducted using experimental [1] or simple numerical tools [2][3][4]. However, this approach is not appropriate if strong interaction between the structure and environment is expected. In contrast, Computational Fluid Dynamics (CFD) models can include the two-way coupling between the fluid and the structure. First CFD attempts were based on two-dimensional calculations and a porous medium approach to account for the effect of the net on the fluid [5][6][7]. This kind of approximation is necessary due to the large length scale difference between the twines of a net and the size of the complete structure, which prevents the resolution of the net on the same numerical grid as the fluid domain. The limitations of the porous medium approach were overcome in [8] by including an additional source term in the Navier-Stokes equations. The model was successfully validated for rigid [8] and flexible [9] net sheets and cages, as well as complete OOA structures [10].
In the CFD model of Martin et al. [10], the hydrodynamic loads were approximated using a screen force model which assumes velocity-dependent drag and lift forces in accordance with [11]. Empirical approaches are required to approximate the involved force coefficients as a function of geometrical properties of the net and fluid properties. Multiple researchers [12,13] focused on fitting polynomials to measured data for this purpose. The disadvantage of experimental approaches is the relatively small range of applicability due to the cost and time restrictions of model test facilities. Alternative approaches focused on analytical analyses of the potential flow around multiple cylinders [14]. In [15], the crossflow principle around tandem cylinders was utilised to formulate a Fourier series expansion for the force coefficients. Here, the Reynolds number, the net solidity and the angle of attack between net and fluid were proposed as characteristic parameters. The additional coefficients arising from their model were not finally fixed. Recently, they were fitted to existing experimental data [8], and the resulting screen force model was successfully applied for CFD simulations. However, the authors highlighted multiple shortcomings of this approach. Particularly, the difference of forces on nets with similar solidity obtained with different twine diameters and lengths could not be modelled correctly due to the dependency of the chosen formula on the net solidity only. Furthermore, a difference between twisted and smooth net materials could not be achieved. These issues were resolved in [16] by utilising a CFD model and Kriging metamodelling to derive polynomial expressions for a new simulation-based screen force model. In comparison to experimentaland analytical-based formulations, this approach allows for the analysis of a large number of different nets with arbitrary geometry. It could be shown that the results of numerical simulations of a small portion of a net panel can be used to approximate the forces on the complete panel with high accuracy. The arising semi-empirical formulae for predicting screen force coefficients were successfully validated against measurements. However, only perpendicular inflow was considered for the drag forces, and an angle of attack of 45 • was investigated for the lift forces. In the present paper, the simulation results in [16] are used to derive new polynomials valid for the complete range of angles of attack. Thus, the application of the simulation-based screen force model to floating OOA structures can be demonstrated.
The remaining paper starts with a brief overview of the numerical framework in Section 2.1. Then, the new hydrodynamic coefficients for the screen force model are derived (Section 2.2) and validated against measurements for waves and current through rigid net panels (Section 3). Finally, the validated model is applied to the simulation of a typical OOA structure in irregular waves to emphasise the applicability of the proposed model (Section 4). Concluding remarks can be found in Section 5.

Fluid Dynamics Solver
The open source CFD code REEF3D [17] is adopted in this paper to simulate the fluidstructure interaction of OOA structures and waves. A one-fluid approach is combined with a direct forcing method in this solver. Here, the two phases, water and air, are distinguished through a signed distance function Φ, and the rigid part of the floating body is described implicitly using the signed distance function Φ s . The conservation laws for mass and momentum in the whole domain are then given as [10] ∇ · u = 0, with u the velocity vector, p the pressure, g the gravitational acceleration vector, S the forcing term for the net representing the momentum loss due to the presence of the net in the fluid, and f the direct forcing term for the rigid body immersed in the fluid domain. In semi-discrete formulation, it can be defined as Here, the smoothed Heaviside step function with = 2.1∆x and ∆x the characteristic cell length is used to ensure a smooth transition between the different phases. This smoothening efficiently avoids pressure jumps which are typically observed for forcing methods. Further, u b represents the rigid body velocity at the grid points in the Eulerian fluid domain: withẋ s the translational, ! s the rotational rigid body velocities and r the vector between the grid point and the bodies centre of gravity. The rigid body velocities are determined using the procedure presented in [18]. In particular, the fluid forces and momenta acting on the floating structure are calculated from integrating the normal and tangential forces at the exact surface of the floating body using interpolated velocity and pressure values. External forces from mooring lines and attached nets are added to the dynamic equations of the rigid body to ensure two way coupled simulations. The free surface between water and air is propagated in space and time using the linear advection equation After each propagation, an equation-based reinitialisation algorithm [19] is applied to keep its signed distance properties. An incremental pressure-correction algorithm [20] is combined with a third-order accurate TVD Runge-Kutta scheme [21] to solve the system of equations explicitly. For this purpose, a rectilinear grid is defined in the numerical domain, and the pressure and velocity fields are defined on the cell centres and faces, respectively. Convection terms are discretised using the finite difference fifth-order accurate weighted essentially nonoscillatory (WENO) scheme [22]. The diffusion and gradient terms are discretised with a second-order accurate central difference scheme. In addition, the diffusion term in the Navier-Stokes equations is propagated in time using the implicit Euler method to remove its influence on the CFL number. The fully parallelized BiCGStab algorithm with geometric multigrid preconditioning from the HYPRE library [23] is applied for solving the implicit equations. A ghost cell approach and the message passing interface (MPI) are implemented for the inter-processor communication.

Simulation-Based Screen Force Model
Offshore aquaculture structures are typically equipped with relatively stiff nets for which reason the deformation of the nets can be neglected. Instead, the location of the net is updated in each time step using the translational and rotational velocities of the rigid floating frame. The forces acting on the net are thereby added as external forces to the rigid body dynamics equations. A two-way coupled approach between fluid and net is enabled by calculating the forces on the net using the solution of the fluid dynamics and including the influence of the net on the fluid in the fluid equations. The forcing strategy developed in [8] is utilised to enforce the boundary conditions at the fluid-net interface. Here, the forcing term S in (2) is defined such that it represents the momentum loss of the fluid while passing the net. The term is first calculated locally on Lagrangian markers, which discretise and follow the net in a Lagrangian manner, and then distributed on the Eulerian fluid grid using a Kernel interpolation. Its magnitude is determined using an extended screen force model including inertia and hydrodynamic forces as well as gravity and buoyancy effects. The gravitational force is approximated from the weight of the net surface A L occupied by the corresponding marker L, and the inertia force F I,L is defined as with m a,L the added mass, n L the unit normal vector of A L , a f ,L the fluid acceleration at x L and a rel,L the relative acceleration vector between fluid and structure. The added mass is taken as the mass of the water volume displaced by the twines within A L . Following the screen force model proposed in [15], the hydrodynamic drag and lift force vector at the Lagrangian location L, F H,L is determined from with n d the normal and n l the tangential direction of the relative vector u rel,L between fluid and solid velocity. The drag and lift coefficients, c d and c l , are generally dependent on the geometrical configuration of the net and the fluid properties. The geometry of a net is typically described by the solidity ratio S r which is defined as the ratio of the solid front area to the total area including the voids between the twines: with l the length and d the diameter of the twines. In [8], the drag and lift coefficients are approximated from the analytically derived Fourier series in [15] and fitting of the free coefficients to experimental data. This approach could be validated successfully but at the same time, shortcomings were highlighted. On the one hand, the limited number of experimental data reduces the reliability of the fitting procedure and on the other hand, the analytical approach neglected the effects of very narrow meshes with thin twines because it only includes the dependency on the solidity ratio and Reynolds number. The authors of this paper overcame both these problems by proposing a simulationbased screen force model [16]. Here, high-fidelity URANS-LES simulations of a small portion of the nets and Kriging metamodelling were utilised to derive polynomials for the hydrodynamic coefficients. In contrast to previous attempts, the dependency of the coefficients on both the diameter and length of the twines was included. Furthermore, the twisted configuration of many aquaculture nets was respected by the model. In the original work, the author concentrated mostly on the drag forces on a net panel perpendicular to the inflow. The approach is extended to arbitrary angles of attack in this paper using the polynomials in [16] for α = 90 • and α = 45 • and quadratic interpolations. The coefficients for the latter are calculated under the assumption that the maximum lift forces occur at α = 45 • and are zero at the two extreme situations of α = 0 • and 90 • . Thus, the drag and lift coefficients are calculated as a function of the angle of attack α using with the coefficients of the quadratic interpolation which are calculated using the polynomials fittings (. indicating a scalar product) c d,90 = p d,90 . u rel , d, l, u rel d, u rel l, d l, u 2 rel , d 2 , l 2 , 1 , [u rel , d, l, u rel d, u rel l, d l, u 2 rel , d 2 , l 2 , u 3 rel , l 3 , u rel d l, u rel d 2 , u rel l 2 , d u 2 rel , d l 2 , l u 2 rel , l d 2 , d 3 , 1], c l,45 = p l,45 . u rel , d, l, u rel d, u rel l, d l, u 2 rel , d 2 , l 2 , 1 .  (14) The obtained polynomials are valid in the interval (0.1 m/s < u rel < 1 m/s, 0.5 mm < d < 3 mm, 0.01 m < l < 0.03 m) due to the simulated nets. The extension of these intervals would be straightforward to achieve as the validated CFD model can be applied to arbitrary geometrical configurations.

Current Cases
Three different experiments are considered to validate the proposed model for current flow through rigid net panels. In [7], the drag and lift forces are measured for a rigid net The width and length of the numerical domain are reduced due to efficiency reasons, but without affecting the results. The free surface is modelled as a symmetry plane. A third experiment on a rigid net panel in current is presented in [25]. The net panel has a width of 1.3 m and a height of 0.7 m. Three angles of attack, i.e., 30 • , 60 • , 90 • , and an inflow velocity of 0.5 m/s is taken into account. The net has a solidity of 0.128 here. The same computational domain and grid sizes are used for the simulations in current. A uniform discretisation of ∆x = 0.07 m is used without presenting a convergence study due to its insignificance for the coupling strategy and the uniform flow pattern (see [8] for a detailed study on this). The time step is controlled by a CFL number of 0.5. Figures 1 and 2 present the numerical and measured drag and lift force coefficients for the described nets as a function of the angle of attack. For comparison reasons, the numerical results obtained with the same CFD code but the screen force model fitted from measurements is included as "old" results (see [8,10]). The computed drag coefficients increase with increasing angle of attack which qualitatively agrees well with the measurements. The predicted lift coefficients indicate that the assumption that the maximum lift forces are obtained for α = 45 • is correct. At larger angles of attack, the flow separates at the frame, and the lift forces reduce. In general, the lift forces are significantly smaller than the drag forces. A slightly better agreement is achieved with the proposed model in comparison to the old approach. The same is true for the lift forces in [24] which were over-predicted by the old model. On the other hand, the lift forces reported in [7] are under-predicted by up to 20%. In general, most deviations are within the 10% error bound which is a sufficient result. For additional validation, the velocity of the damped flow behind the net is compared to the experimental results in Figure 3. Here, only cases with perpendicular inflow but different inflow velocities are considered. The numerical model is capable of predicting the right amount of damping for all considered inflow velocities as the deviations are below 10%. This correct behaviour is particularly important for inflow velocities larger than 0.5 m/s because here, the net has a significant effect on the fluid flow.

Wave Cases
Next, the drag forces on a rigid net panel perpendicular to incident regular waves are simulated and compared to experimental results. In [1], the drag forces on a rigid sheet in regular waves are investigated. Here, a rigid net panel of 1.0 m × 0.5 m with S r = 0.2 is placed in the centre of a wave flume. Three different waves with the wave heights H = 0.1 m, 0.15 m and 0.2 m and a wave period of T = 1.4 s are investigated in this study. They are numerically modelled using Stokes' second-order theory. The numerical domain is adapted from the experiments with a uniform discretisation of ∆x = 0.07 m. The time step is controlled by a CFL number of 0.3. The time series of the numerically predicted drag forces on the panel are compared to the experiments in Figure 4. Generally, a good agreement is observed for both, the frequency and peak values of the force signal. The simulation tends to over-predict the crests with increasing wave height, but the deviations are still below 10%. It is further interesting to notice, that the forces tend to double between each of the wave inputs which highlights the non-linearity of the expected forces on nets in waves.

Application to the Motion of a Semi-Submersible OOA Structure in Irregular Waves
The applicability of the proposed model to the simulation of OOA structures shall be provided in the following. The study is based on the model test description in [26] for a 1:120 model of the Ocean Farm 1. The rigid structure consists of 9 columns with pontoons attached at their bottom. As illustrated in Figure 5, they form a hexadecagon with a diameter of approximately 1 m and one column placed centric. The simulation setup is taken from [10] where the same structure was investigated using the screen force model fitted to experimental instead of simulated data. Four taut mooring lines and a rigid net with a solidity of 0.145 are attached to the structure. The wave input is a JONSWAP spectrum with a significant wave height of 0.1 m and a peak period of 1.4 s (see also Figure 6a). As results, the linear transfer functions (RAO) are calculated using with S motion and S wave the power spectra from auto-correlation analyses. As shown in Figure 6b, the heave response increases with decreasing wave frequency and decreases rapidly for high-frequency excitations. This minor response in high-frequency excitation indicates that the overall system is highly-damped. Furthermore, the surge motion increases with decreasing encounter frequency and follows the wave envelope for frequencies around 0.55 Hz (see Figure 6c). The structure is expected to stay close to the equilibrium position in shorter waves similar to the heave motion. In contrast, the pitch motion in Figure 6d is relatively large for the whole range of investigated wave frequencies. This might be caused by a rather horizontally acting mooring system with less resistance forces in the rotational direction. A possible resonance is indicated at a frequency of 0.55 Hz where the largest amount of wave energy is located. Finally, the maximum tension forces in the front mooring lines increase with amplified structural motion. Thus, a maximum is predicted at a frequency close to the heave and surge maxima. On the one hand, this again indicates the horizontal direction of the mooring resistance forces. On the other hand, the mooring reaction forces seem to be the driving forces for the expected dynamics of the complete semi-submersible OOA structure.

Conclusions
In this paper, a new CFD model for the simulation of OOA structures in waves is proposed. It is based on the two-way coupled solution of the Navier-Stokes equations with additional source terms to account for the rigid floating body dynamics and the damping effect of the attached nets. A new simulation-based screen force model is presented to calculate the force coefficients in the net model. It bases upon the work in [16] and extends the original approach to arbitrary angles of attack using quadratic interpolations.
The validation process for current and wave cases shows that this approach is capable of predicting the forces on and velocity reduction behind rigid net panels of various geometry. In current, the model predicts the drag forces and velocity reduction within a 10% error band, whereas it tends to under-predict the lift forces by up to 20%. In waves, the model tends to over-predict the crests with increasing wave height, but the deviations are also within a 10% error band. Overall, the accuracy is in accordance with previous results obtained with the model. The advantage of this model is the applicability to an infinite range of geometrical net properties and flow conditions, which shall be demonstrated within future research.
The application to a semi-submersible OOA structure further highlights the practicability. Within future research, the complete CFD-based model is applied to different types of OOA structures and the simulation-based screen force model is further improved using additional simulation results.

Conflicts of Interest:
The authors declare no conflict of interest.