A numerical framework for modelling the dynamics of open ocean aquaculture structures in viscous fluids

This paper presents a complete numerical framework for modelling open ocean aquaculture structures in waves and current using Computational Fluid Dynamics (CFD). A structural dynamics model is incorporated to account for the motions and deformations of the net. It is based on the lumped mass method, a non-linear material law and implicit time step advancing. The presence of the porous net is considered in the momentum equations of the fluid using a forcing term based on Lagrangian-Eulerian coupling and the acting forces on the net. The proposed framework is suitable for simulating the interaction of nets of arbitrary geometry and under large motion with fluids including complex free surfaces. This is in contrast to existing models which either neglect important non-linearities, the physical interaction with the fluid or are limited to certain net geometries. In addition, the fluidstructure interaction of floating objects with mooring lines, nets and fluid is accounted for in the model. A new floating algorithm is presented for the interaction of the fluid and the rigid structure. It is based on a continuous direct forcing immersed boundary method and a level set representation of the object in the Eulerian fluid domain. This effectively avoids computationally expensive reconstruction processes of existing approaches and enables the application to large three-dimensional structures. The complete numerical framework is first validated against existing measurements for forces on rigid and flexible nets, net deformations and moored-floating structures with and without a net in waves. Then, a semi-submersible and a mobile floating open ocean aquaculture structure are investigated, and the possibilities of the numerical approach are highlighted.


Introduction
The demand for aquatic food products is expected to increase by thirty million tonnes until 2050 (Ferreira et al., 2014), and aquaculture production has to play a significant role in satisfying the demand to avoid overfishing. In recent years, the traditional aquaculture industry has faced increasing criticism for the environmental impacts of their near-coast structures on the surrounding marine habitat (Grigorakis and Rigos, 2011). This raises the need for alternative concepts in the aquaculture industry. One of the most promising and economically valuable ideas involves the increase of the dimensions of the structure in combination with the relocation of the structures offshore. However, this implies larger loads on the structure and higher risk of fish loss due to the influence of strong current and larger waves. Open ocean aquaculture (OOA) is concerned with the adaptation of fish cages and the risks during operations to an environment with significant exposure to wave action and severe sea conditions (Ferreira et al., 2014). As a consequence, alternative design choices adapted from offshore related fields of engineering arose during recent years (Chu et al., 2020). This includes semi-submerged and submerged fish cages (Xu et al., 2013) to reduce the loads as well as floating rigid structures with one or multiple nets attached. These new types of structures require more advanced tools in the design process to better understand the dominating factors for eventual structural failure and fish loss. Numerical modelling is a relatively inexpensive and flexible way of assessing these factors and reduce the risks if appropriate approaches are chosen during establishing the frameworks for OOA structures.
Historically, computational methods developed for investigating aquaculture systems relied on linear potential theory for load calculation and empirical formulae for estimating the velocity reduction through the net. The most comprehensive study applying these methods was performed by Kristiansen and Faltinsen (2015). They validated a lumped mass net model coupled to a dynamic beam equation for the floater, linear wave theory to approximate the excitation forces and the formula of Løland (1991) to account for the shading effect of the net against experiments of a traditional aquaculture structure. A reasonable agreement could be achieved for small wave steepness and low current which is in accordance with linear theory. For the same type of structure, Shen et al. (2018) validated their numerical model, which is based on the same assumptions, against measurements in both regular and irregular waves. They concluded that in severe sea states, the deformation of the net is of a more limiting factor in the design choices than the stresses in the floater. First numerical investigations of OOA structures were presented in (Fredriksson et al., 2005(Fredriksson et al., , 2003 using a finite element method for the net and linear potential theory for the excitation forces due to current and waves. The considered structure was a spar buoy with an octagonal rim held together by tensioned stays woven into the net. No shading effect was taken into account. The authors stressed the importance of non-linear effects of extreme waves and fluid-structure interaction for the correct prediction of the structural dynamics. Xu et al. (2013) compared a numerical model based on potential theory and a rigid floater for submersible and moored floating net cages with physical model tests. The authors highlighted that the lowering of the structure resulted in smaller deformations of the net and reduced mooring line tension forces compared with the floating configuration. Also, varying wave steepness have minor effects on the structural loads in this condition. More recently, Li et al. (2018) presented a numerical analysis of the concept of a vessel-shaped cage system using a panel method and linear wave theory as the basis. No validation against measurements was provided. However, it was concluded that a deformable net model is necessary and non-linear effects have to be taken into account if offshore environmental conditions are to be investigated. These conditions are characterised by large loads on the structure and hence, strong non-linear fluid-structure interaction. Thus, the existing numerical tools cannot be regarded as appropriate for investigating OOA structures. In contrast, computational fluid dynamics (CFD) offers a two-way coupled numerical approach that can be applied to understand the environmental loading and structural response by modelling the hydrodynamic forces affecting the dynamics of the floating structure, the flexible net and mooring system, and their effect on the surrounding fluid. To the best of the authors' knowledge, the only numerical approach so far was the two-dimensional model presented by Chen and Christensen (2018). They solved the incompressible Navier-Stokes equations in a two-dimensional numerical wave tank and coupled it to a lumped mass net model using a porous medium approach. They validated their model against experiments for a net attached to a moored-floating cylinder in waves and current and showed promising results. However, the approach was not tested for three-dimensional structures such as OOA structures.
The complex problem of modelling the whole system can be split into the sub-problems of solving the interaction between the fluid and a rigid moored-floating structure and the flexible porous net, respectively. Several approaches exist to solve the fluid-structure interaction (FSI) involving rigid bodies only. Moving mesh methods fit a Eulerian grid to the body which is distorted by the changing position of the rigid body. Accurate results in the vicinity of the body can be achieved, but the method faces problems for large body motion and complicated structures. The grid distortion can be avoided effectively by applying the concept of dynamic overset grids. Here, the grid around the body moves with the structure, and the fluid information is transferred to a fixed background grid at a certain distance from the body. Amongst others, Carrica et al. (2007) and Chen et al. (2019) successfully developed numerical towing and wave tanks using dynamic overset grids. The drawback of this method is the increase of computational costs due to the interpolations necessary in the overlapping region of the grid, the parallelisation of the complete solver, and eventual instabilities due to incomplete interpolation stencils. Alternative methods based on an implicit representation of the structure in the Eulerian fluid grid have been developed to provide the computational efficiency necessary for simulating large three-dimensional structures. The most important implicit methods are based on the immersed boundary method by Peskin (1977) for simulating elastic membranes in fluid. Later, Fadlun et al. (2000) introduced this idea to rigid bodies as the direct forcing method. Here, the boundary conditions between fluid and solid are respected by incorporating an additional source term in the momentum equations. The term is calculated on fluid grid points near the surface using a reconstructed solution from the fluid domain and the known velocity at the nearby solid surface. It was pointed out by Uhlmann (2005) that this reconstruction procedure can lead to spikes in the time series of the forces because the reconstruction stencil changes as the body moves. Yang and Stern (2015) presented a solution to this problem by including an overset grid so that the reconstruction stencils remain constant in time. This, however, brings along the aforementioned challenges of these type of grids. Therefore, Uhlmann (2005) proposed a continuous version of the direct forcing method by calculating the forcing term on Lagrangian markers which discretise the surface of the structure. The terms are then distributed on the Eulerian domain using the interpolation procedure of the original immersed boundary method (Peskin, 1977). The introduced smearing effectively removes force spikes but still keeps the nominal order of accuracy for the FSI problem. Additional computational costs arise from the back-and forthtransformations between the Eulerian grid and the Lagrangian markers. These costs can be avoided by using a completely Eulerian calculation as developed by e.g. Yang (2018). Here, the Eulerian grid is split into a fluid and solid domain. The momentum equation is solved in the whole domain including a forcing term. This additional term is calculated from the rigid body velocities at each grid point and smeared over the fluid-solid interface using a smoothed Heaviside step function. The author validated the solver against two-dimensional benchmark cases and showed high stability and accuracy. Here, the author used markers moving with the body for reconstructing the Heaviside step function in each time step. In the present paper, the efficiency of this approach is further enhanced by utilising a ray-casting method and a level-set function to distinguish between fluid and solid.
The second FSI problem is concerned with the interaction of the fluid and the net. The net is an elastic structure with non-linear material properties (Lader and Fredheim, 2006) undergoing potentially large motions and deformations. The lumped mass method was established as the most common solution for solving net dynamics. Originally presented by Lader and Fredheim (2006), the lumped mass method discretises the net into massless bars with mass knots in between. The dynamics of the knots are found by solving Newton's second law and a Runge-Kutta time integration. In (Bi et al., 2014a,b), the method was validated for flexible net sheets and cylinders in steady current flow. However, the explicit time integration and missing fulfilment of the constitutive equations within each time step lead to severe time step restrictions. Implicit methods are therefore more suitable for FSI problems involving a coupled fluid solution. In , an implicit quasi-static net model was proposed where the force equilibrium is solved for each knot with additional constraints on the connectivity of the bars. The missing time-stepping reduces the cost and simplifies the coupling to the fluid solver. However, the lack of dynamic effects prevents the application to large deformation problems including snap loads. A dynamic implicit net model based on the satisfaction of the kinematic relation between knot position and bar length was introduced in (LeBris and Marichal, 1998). The original approach assumed inelastic material which can lead to ill-posed problems due to high condition numbers. Marichal (2003) successfully overcame this drawback by including elastic material into the model. However, their derivation relied on a linear material assumption and linearised equations. This is a severe drawback considering the non-linearity of net material (Lader and Fredheim, 2006). Therefore, the authors of this paper included non-linear material properties which led a non-linear system of equations for the unknown tension forces . The method was successfully validated against model tests for top-fixed deforming net sheets and cages in waves and current. In the present paper, this model is extended to dynamic problems involving moving mounting points.
In contrast to conventional membranes, nets have high porosity and consist of multiple individual twines which are passed by the fluid. The length scale of the flow around each twine is significantly smaller than the length scale of the flow around the whole floating structure. This prevents the resolution of the net on the same numerical grid as the fluid domain, and an alternative representation of the FSI between net and fluid has to be introduced. One possible representation is the definition of a porous medium around the net. The fluid momentum equation is solved in the whole domain including an additional resistance coefficient in the porous zone (see (Bi et al., 2014a;Christensen, 2016, 2017;Patursson et al., 2010;Zhao et al., 2014) for various definitions of the coefficient). As shown in , this approach is not appropriate for complex shapes and deformations and lacks a physical basis. The authors of this paper proposed an alternative method based on the same principles as the immersed boundary methods by Peskin (1977) and Uhlmann (2005). Here, a forcing term is calculated from the hydrodynamic loads on the net and distributed on the fluid domain using an interpolation kernel. This method was successfully validated for rigid  and deforming nets  and will be utilised in the present paper.
Based on the review of existing approaches, this paper presents itself as the first attempt to establish a numerical framework using CFD for OOA structures. All parts of the framework are chosen or newly developed aiming to simulate large, three-dimensional structures including fluid-structure, fluid-net and net-structure interaction in an efficient manner. The paper is structured as follows: section 2 presents the different modules of the complete numerical framework. Several validation cases are presented in section 3. The validated numerical model is then applied to two typical OOA structures in section 4. The paper concludes with final remarks in section 5.

Numerical framework
The different parts of the proposed numerical model are introduced in the following. The framework includes a two-phase numerical wave tank, a floating algorithm for modelling the interaction of rigid structures and fluid and an implicit solver for the net dynamics. Details regarding the coupling of the net solver to the fluid solver and the rigid body solver are also given. A flowchart at the end of the section (Fig. 4) provides an overview of the complete framework.

Fluid dynamics solver
The dynamics of any incompressible fluid obeys the conservation of mass and momentum. These conservation laws are described by the three-dimensional Navier-Stokes equations and continuity equation which are written in the convective form as Here, u is the velocity vector, ρ is the density, p is the pressure, ν represents the kinematic and turbulent viscosity and g is the gravitational acceleration vector. The effect of turbulence is incorporated by adding turbulent viscosity to the diffusion term using the Boussinesq approximation and a modified k-ω turbulence model (Bihs et al., 2016). Following a one-fluid approach, the two phases, air and water, are covered by a single set of equations but space and time dependent material distributions. The transition between the phases is implicitly represented by the zero level set of the smooth signed distance function Φ (Osher and Sethian, 1988). The linear advection equation is solved for propagating Φ in space and time. After each propagation, the reinitialisation equation of Sussman et al. (1994) is solved to keep the signed distance properties of Φ. The density and viscosity is then calculated using with w indicating water and a air properties. Further, the smoothed Heaviside step function H is given by with = 2.1∆x and ∆x the characteristic length of the discrete domain in the vicinity of each evaluation point. The set of equations (1) -(3) is solved on a staggered rectilinear grid using finite differences. Fifth-order accurate weighted essentially non-oscillatory (WENO) schemes (Jiang and Peng, 2000;Jiang and Shu, 1996) adapted to non-uniform point distances are applied for convection terms, and diffusion terms are discretised with second-order accurate central differences. Convection and source terms are treated explicitly with the third-order accurate TVD Runge-Kutta scheme of Shu and Osher (1988), and an implicit Euler method is applied for the temporal discretisation of the viscous term. This effectively removes a quadratic reciprocal dependency on the cell size from the CFL condition (Bihs et al., 2016). An incremental pressure-correction algorithm (Timmermans et al., 1996) is used for solving system (1)-(2). In each k-th Runge-Kutta sub-step, a velocity field is predicted using the pressure gradients of the previous step: with α k = 1.0, 1/4, 2/3, β k = 0.0, 3/4, 1/3 and k = 1, 2, 3. Afterwards, the Poisson equation is solved for the pressure correction term p corr utilising a fully parallelized BiCGStab algorithm with geometric multigrid preconditioning from the HYPRE library (van der Vorst, 1992). An n-halo decomposition strategy and the message passing interface (MPI) handles inter-processor communication. The pressure and divergence free velocity fields are finally calculated from

Including rigid body motion in the fluid solver
A new continuous direct forcing approach is implemented to account for the presence of rigid floating objects in the fluid solver. The same level-set routines as for the free surface are utilised. Hence, high efficiency through parallelised routines is ensured and simulations of large three-dimensional structures are possible. This is in contrast to previous research on continuous methods which mostly concentrated on two-dimensional benchmark cases. The floating object is transferred to the solver as an STL geometry consisting of multiple non-connected triangles. This information is sufficient to create a signed distance field Φ s representing the geometry in the Eulerian fluid domain by applying a ray casting algorithm (Bihs et al., 2017) to receive inside-outside information near the body and the reinitialisation algorithm of Sussman et al. (1994). An example of a complex STL geometry and its level set representation is shown in Fig. 1. The generated level set function Φ s is used for distinguishing between fluid and solid domain by extending (4) and (5): with s indicating solid.  Following the derivation of Yang (2018), the conservation laws with the forcing term hold in the solid domain. Here, P(u) represents the operator for projecting the velocity field into a divergence free rigid body velocity field. Comparing (13) with (1) and (2) reveals that the only difference between these two sets of equations is the term f and the diffusion term, respectively. Therefore, a single set of equations can be solved in the whole domain using H(Φ s ) for representing the transition of fluid to solid and the additional term f for incorporating the correct boundary conditions at the interface. At a discrete level, f reads at the new time step n + 1 The velocity at the new time step is unknown a priori. To overcome this issue and avoid expensive implicit calculations, the valid approximation P(u (n) ) = u (n) is made and the pressure is taken from the previous time step as a good approximation. Then, and by comparing with (7), it can be identified that A good approximation of the updated velocity field is u ( * ) itself. Therefore, the predictor step (7) is first executed without the forcing term. Then, f ( * ) is calculated from and added to the predicted velocity field before solving the Poisson equation (8).
For the calculation of the rigid body velocity field, the translational motion of the rigid body is described by Newton's second law, and the rotational motion is described in a bodyfixed coordinate system using the Euler parameter e = (e 0 , e 1 , e 2 , e 3 ) T with the property e T e = 1. Their relation to Euler angles is provided in A. Following Shivarama and Fahrenthold (2004), a first-order Hamiltonian system can be derived. The translational equations are converted into a system of first-order differential equations as well. Hence, the rigid body system can be integrated with the same explicit Runge-Kutta method as the fluid solver. The body forces and momenta are calculated on the triangulated surface using a trilinear interpolation of the fluid properties and an integration over all N triangles: Here, n is the surface normal vector, τ is the viscous stress tensor and r represents the distance vector to the centre of gravity. The calculated moments are transferred into the body-fixed coordinate system before solving the system. Once the body velocities are calculated, the projection can be calculated as withẋ s the translational and ω s the rotational rigid body velocity vector in the inertial reference frame.

Modelling the net dynamics
The main concept of the previously proposed net model  is presented in the following. The net is discretised in a finite number of mass points (knots) connected by non-linear elastic bars. Due to reasons to be mentioned in the next section, a structural element is defined as the combination of four points and their four connecting bars and covers multiple physical meshes of the net (see Fig. 2).
Figure 2: Illustration of the definitions for the net model: structural elements consist of four knots (thick dots) connected with bars (thick black lines). Each element (hatched area) covers multiple physical meshes (thin grey lines). The contribution of the structural elements to each knot is shown in matching colours.
A system of equations is formulated for the dynamics of the knots by distributing the external forces F ex from each structural element to the attached knots. This leads to the dynamic equilibrium for each knot x i and its N i neighbouring knots. Here, dots indicate temporal derivatives and T ij represents the tension force vector of each bar: with T ij the tension force magnitude and b ij the unit vector of the bar. The mass matrix m i is calculated considering the surrounding N S,i structural elements using with m air,s the mass of the partial element, n s the unit normal vector of the element pointing in relative velocity direction and m a,s the added mass contribution. The external force vector consists of gravity and buoyancy forces as well as hydrodynamic forces. The latter consist of inertia forces I due to the fluid acceleration a f , and velocity related forces D which are calculated using the screen force model of Kristiansen and Faltinsen (2012). In the inertia system of the fluid, D can be split into drag and lift force components in normal (n d ) and tangential (n l ) direction of the local relative velocity vector u rel,s = u f,s −ẋ s : The coefficients c d and c l are calculated from a truncated Fourier series expanded for their dependency on the angle of attack α between fluid velocity vector and structural element vector: The definition of the constants c d,0 and c l, π 4 can be taken from (Kristiansen and Faltinsen, 2012). The determination of the Fourier coefficients is based on non-linear fitting to experimental data. This raises issues as most experimental data is normalised by the undisturbed inflow velocity. However, the velocity at the location where the forces are evaluated in the numerical simulation is disturbed by the presence of the structure. In , an equation is derived to approximate the inflow velocity without the structure. Based on momentum equilibria considerations, the undisturbed inflow velocity u ∞ is approximated by solving the intrinsic equation An implicit solution for the structural dynamics problem is found by starting from the trivial relation between the position of the knots x i and x j and the length of the bar l ij between them: According to Lader and Fredheim (2006), the material of nets follow the non-linear constitutive law with l 0,ij the unstretched bar length and (C 1 , C 2 ) = (1160 N, 37300 N ) for squared meshes. Inserting (30) in the right-hand side of (29) yields The left-hand side of (29) is replaced by the dynamic equilibrium (21). This is done by expressing the position vectors by its accelerations using high-order backward finite differences and replacing them with forces (see B for the derivation). Thus, a non-linear expression f is found for each bar b ij : with T the unknown global vector of tension force magnitudes. A system of equations can be formulated and solved using the improved Newton's method (Chun, 2005) with F the vector of the expressions (32) and J its Jacobian matrix (see  for the expression). In practice, (33) converges well if the initial condition is chosen properly. More details can be found in . The converged result of (33) can then be used for calculating the acceleration, velocity and position of the knots in a straightforward manner.

Coupling the net solver to the fluid solver
The solidity Sn of a net is defined as the ratio of solid front area to the total area and is calculated as (Fredheim, 2005) with d t the twine diameter and l t the twine length. Sn typically varies between 0.1 and 0.3 for aquaculture nets. The resulting difference between the length scale of each mesh and the total structure prevents the resolution of the whole net in the discrete space. This raises the need for an alternative fulfilment of the boundary conditions at the fluid-structure interface because the physical fluid-structure interface is not present in the simulation. Following the Lagrangian approach developed in , these boundary conditions are replaced by a source term S which expresses the physical loss of fluid momentum due to the presence of the net, which leads to a pressure jump. The term is determined from the known external forces on the net and is added to the Navier-Stokes equations before the pressure correction step. The momentum loss is assumed to be uniformly distributed over the net surface regardless of the difference between twines and voids. This implication requires uniformly distributed markers holding S and that the area which is covered by each marker is nearly equal to the cell size of the encircling fluid domain. However, it is not practical to evaluate the structural dynamics on a scale where the knots fulfil these requirements. Therefore, additional Lagrangian markers are introduced to distribute S. They are defined by splitting the structural elements in triangles according to the Eulerian grid size in their vicinity. The Lagrangian markers are then defined in the geometrical centres of each triangle. An example of the distribution of the markers on the discrete net structure is shown in Fig. 3. The forces are distributed on the fluid grid points x e using the interpolation with L e the number of markers within a defined kernel around x e which is taken from (Peskin, 1977): The forces vectors s(x L ) at the marker with position x L = (x L , y L , z L ) are calculated by integrating the external forces in (21) over the triangle area A L :

Coupling net and rigid body dynamics
Additional remarks are given regarding the coupling of the fluid-net and fluid-rigid structure solver. At first, it is decided to explicitly couple the net dynamic solver to the remaining algorithm due to the different time advancement procedures. Therefore, the fluid velocity field of the previous time step is chosen for calculating the loading on the nets. The field is interpolated to the points where the external forces are evaluated using the same kernel function as provided in section 2.4. After the external forces are evaluated, the new tension force distribution is calculated and used for updating the dynamic properties of the knots. Furthermore, the topmost knots are assigned to the rigid body motion to ensure a tight coupling. For this purpose, it is sufficient to replace the third or fourth term of the dynamic equilibria (32) with the known acceleration of the rigid body at these points. Two-way coupling is enabled by adding the tension forces acting on the top knots to the rigid body motion solver as external loading. It is also to be noticed that the later applied mooring system is incorporated in the same manner. It is referred to  for more details.

Validation of the net model in current
The structural net model and its coupling to the fluid solver are validated for steady current flows in the following. A validation involving wave loads can be found in . At first, the forces on rigid net sheets, i.e. without considering net deformation, are investigated. Afterwards, the deformation of cylindrical cages and net sheets are compared to experiments in steady current flow. The main focus in these cases lies on the correct deformation and shading of and forces on nets with different geometrical properties.  (25), and the velocity in the wake of the sheet is measured 0.715 m behind the sheet and 0.9 m below the free surface.

Drag and lift forces on a rigid net sheet
The numerical domain has the dimensions 12 × 7 × 5.6 m and is shown in Fig. 5. The width and length of the tank are reduced due to efficiency reasons, but without affecting the results. The free surface is modelled as a symmetry plane. A uniform discretisation with ∆x = 0.07 m is used. Convergence studies are not presented as the coupling is insignificantly dependent on the fluid grid size and the flow is mostly uniform.   Fig. 6b shows the drag force distribution for varying angles of attack for u ∞ = 1.0 m/s. The forces increase with the net solidity and decreasing α due to an increasing area of attack, and they increase quadratically with the inflow velocity as expected. The maximum deviations between numerics and experiment are 20% for the smallest velocities and solidities. A better agreement is seen for the lift forces at α = 45 • (Fig. 6c) with deviations smaller than 10%. For Sn = 0.16, Føre et al. (2020) report inaccurate measurements because the expected lift forces should be larger than for Sn = 0.15. Therefore, the numerical model seems to predict a more accurate result here. The distributions of the lift forces generally follow the drag forces due to the same reasons. However, the lift forces reach a peak value at around 45 • and decrease again for larger angles because the net acts as a closed surface with flow detachment rather than as a porous sheet in these cases.
The flow decelerates through the net as indicated in Fig. 6d. The shading effect increases with the solidity of the net which is respected in the numerical model through the dependency on the drag forces. Further, the model predicts a lower influence of the inflow velocity on the relative velocity reduction for all solidities. The same effect is reported in the experiments, except for in the case of the largest solidity. Additional measurements should be conducted to clarify the eventual existence of changing shading effects for high solidity nets. Generally, the numerical model over-predicts the momentum loss but the deviations stay within a 10% error bound.

Forces on and deformation of a cylindrical net cage
The complexity of the simulation is increased by considering the net deformation within the FSI solver. The physical model tests of Bi et al. (2014b) in which a top-fixed cylindrical net cage was investigated in steady current flow are considered for validation. The domain size is 5 × 0.45 × 0.4 m which is adopted from the experimental setup. It is discretised using ∆x = 0.01 m. The free surface is modelled as a symmetry plane for computational efficiency. The comparison of simulated and measured results for the sinker weight of 8 g is presented in Fig. 8. The drag forces qualitatively follow the theoretical distribution of forces on a cylinder for different inflow velocities. The velocity reduced by the presence of the net is well captured numerically as seen in Fig. 8b. Here, the model under-predicts the experimental data by up to 10%. Reported pictures from the experiments are used to compare the deformation of the cage though this implicates large uncertainties. The qualitative comparison of the rear vertical centreline is presented in Fig. 8c. As the velocity increases, larger hydrodynamic loads increase the deformation of the net structure. The numerical model presents a good approximation of the physical model test.
Similar effects can be observed if the sinker weight is increased (see Fig. 9a). The additional gravity forces counteract the hydrodynamic forces and, hence, decrease the deformation of the net in the flow direction. In principle, the numerical model agrees well with the experiment. It is however noticed that two different deformation pictures are shown for the same setup of u ∞ = 0.242 m/s and 8 g sinker weight in (Bi et al., 2014b). A possible explanation is the strong vortex shedding for this configuration, leading to oscillatory motion of the structure. The pictures might be taken from two different time instances. This explanation is confirmed numerically because it is possible to find matching distributions for both pictures within the converged solution. In addition, Fig. 9b presents the distribution of the x-velocity through the net. It clearly shows the momentum loss of the fluid through the front and back of the net. The simulation predicts similar velocity reduction for all weight configurations, whereas the experiments report an increasing reduction for more deformed nets with larger angles of attack between fluid and structure. This effect is a new observation presented by Bi et al. (2014b) but stays in contrast to previous measurements of, e.g. (Patursson et al., 2010) which showed a minor dependency of the angle of attack on the velocity reduction. The latter could also be confirmed numerically by the present approach .       Bardestani and Faltinsen (2013) presented model test results for a deforming net sheet in steady current flow. In contrast to the previous case, tension forces acting at the top of the net were measured for different net geometries, sinker weights and inflow velocities. The proper determination of the structural forces is the result of accurate replication of the loads on and deformation of the structure itself. Hence, it is a reliable measure for validating numerical models for the interaction of net and floating structures. In the experiment, a wave tank at the Department of Marine Technology at the Norwegian University of Science and Technology was used. The tank has the dimensions 13.67 × 0.6 × 1.3 m with 1 m water depth. The same tank size is used in the numerical simulations, and a cell size of 0.05 m is considered. All nets have the size 0.51 × 0.76 m and are placed in the middle of the tank with the top of the net fixed close to the free surface, which is replaced by a symmetry plane in the simulation to save computational costs. The considered nets consist of square meshes modelled by 10 × 10 elements and have a solidity of 0.16, 0.19 and 0.23 with a twine diameter of 2.5, 2.5 and 1.8 mm. The corresponding twine lengths are calculated from (34). A cylindrical sinker was attached to the bottom of the net. Its weight varies between 1.2, 1.4 and 1.6 kg. The simulations are executed until a steady state is reached and then compared to the experimental data. The resulting tension forces over changing inflow velocities in Fig. 11 reveal three characteristics for the investigated setup. At first, a non-linear increase of the tension force can be observed for increasing velocities irrespective of the solidity or sinker weight. This is caused by the increase of hydrodynamic forces as shown above. Further, the tension increases with increasing sinker weight even though the deformation of the net decreases as indicated in Fig. 10. This shows that the increase of gravity forces exceeds the decrease in hydrodynamic forces. Hence, the sinker weight significantly influences the behaviour of the system. This is confirmed by the last observation that the tension forces increase with the net solidity, due to the increase of net area, but the effect is negligible here. Qualitatively, the numerical framework agrees with the experiments on these phenomena, and most predicted results are within a 10% error bound. The largest deviations are seen for the highest velocity and solidity with errors of up to 20%. A reasonable explanation for the deviations is the fact that the flow around the sinker is not resolved in the simulation. Hence, the attachment between sinker and net is missing, and the important drag and inertia forces are approximated using Morison's formula.  Figure 11: Comparison of the numerical and experimental tension forces at the top of a net sheet in steady current flow. In (Bardestani and Faltinsen, 2013), the results are normalised by the tension forces without inflow T 0 . The shown data are obtained after multiplication with the numerically calculated T 0 . The numbers in the legend refer to the sinker weight in kilogram.

Moored-floating cylinder in waves
The presented floating algorithm in section 2.2 is validated against measurements of a mooredfloating cylinder in waves. The setup replicates the experiments of Kristiansen (2010) which were performed in the same tank as described in section 3.3. A cylinder with diameter 0.1 m, the same length as the width of the tank and a mass of 3.94 kg/m is placed in the tank. A mooring system, consisting of ropes and springs with stiffness 151.2 N/m 2 and pre-tension 38.1 N/m, is attached to the cylinder. The other end of the lines is coupled to a pulley 2.43 m away from the cylinder at a height of 0.136 m above the free surface. A two-dimensional setup of the experiments is investigated in the numerical wave tank shown in Fig. 12. The waves are generated in a wave generation zone of one wavelength, and a numerical beach damps the waves at the end of the tank. A description of these methods can be found in (Bihs et al., 2016). The mooring system is modelled as two springs mounted at the centre of the cylinder. The cylinder is shown in blue, the mooring lines in green and λ is the wavelength.
First, a convergence study is conducted using decays test of the moored-floating cylinder (see Fig. 13). In the experiments, only a surge decay of the moored-floating cylinder is reported. The comparison of the time series using numerical grids with ∆x = 0.01, 0.0075 and 0.005 m is presented in the Figs. 13a and 13b. For all chosen cell sizes, the numerical model captures the first peak well. Numerical damping results in under-prediction of the subsequent amplitude for coarser grids. Similar observations are indicated for the free heave test (Fig. 13c). Although the medium grid seems to accurately predict the amplitudes, further refinement is needed for a convergence of the phase. As a consequence, the finest cell size is chosen for further testing.  Num ∆x =0.01m Num ∆x =0.0075m Num ∆x =0.005m (c) Numerical heave decay test without mooring system.  In shorter waves, the surge motion oscillates with the encounter frequency plus a sub-harmonic component at half the frequency. This effect is captured well by the numerical model as seen in Fig. 15, which shows the power spectrum of the surge response for the shortest wave. Further, the model follows the linear relation between the wave period and surge amplitude response present for short waves. As shown in Fig. 16, the object thereby damps the waves almost completely as the wave height in the wake is significantly reduced. The heave motion gradually increases as the waves increase in height. Here, linear theory significantly overpredicts the heave response, whereas the numerical solution agrees well with the experiment. In longer wave periods, the heave response amplitude reaches a value slightly above 1 as the cylinder follows the surface of the long waves. As indicated in Fig. 16, highly non-linear fluidstructure interaction occurs in these conditions. Around the surge resonance frequency, the numerical model predicts amplitudes close to the experiments, whereas theoretical formulae typically overestimate the response. Following the work of Kristiansen (2010), further analysis of the wave excitation forces can be conducted by comparing the significant amplitude components of the acceleration signal. This is caused by the direct link of acceleration and hydrodynamic load through Newton's second law. The resulting linear, second and third harmonic components of the acceleration amplitudes are provided in Fig. 17 using discrete Fourier analyses. In surge direction, the excitation is mainly driven by the linear component. High-order components become more relevant for longer waves, in particular around the surge resonance. This increase might be caused by wave overtopping and viscous effects due to flow separation. For the vertical accelerations, second-order harmonics exceed the linear part for wave periods between 0.8 s and 1.0 s. It is noticed by Kristiansen (2010) that this effect occurs because the second component occurs at a frequency close to the natural heave frequency of the system. As a consequence, the disregard of these high-order components can lead to a significant underestimation of the loads on the structure.

Moored-floating cylinder with a net attached in waves
As a final step of the validation process, the complete numerical framework is tested against measurements of a moored-floating cylinder with an attached net in waves (Bardestani and Faltinsen, 2013). The considered experiment is a combination of previous cases. The moored floater is the same as in section 3.4, and the attached net with solidity 0.23 and sinker weight 1.64 kg is taken from section 3.3. Tension forces acting in the topmost twines are reported for regular waves with a wave steepness of 1/14 and wave periods between 0.4 s and 1.3 s. The dimensions of the numerical wave tank are adapted from the previous section, but a third dimension is added due to the presence of the net. A grid size of 0.005 m is defined for the x-and z-direction and a coarser resolution of 0.08 m is chosen in y-direction because no three-dimensional effects are expected. The complete setup is shown in Fig. 18.   19 compares the predicted maximum tension forces with the experimental data. For small waves, the maximum tension forces are of similar magnitude as in a hydrostatic fluid due to the small motion of the cylinder in these waves. As the wave height increases and the wave period approaches the eigenperiod of the system, the maximum tension forces are approximately five times higher than in the hydrostatic condition and snap loads occur. These loads arise from the relative motion between the cylinder and sinker. Typically, when the cylinder is in a wave trough, the maximum elongation of the net reduces and the net becomes slack (see Fig. 20a). When the wave trough passes and the cylinder is accelerated upwards by the following wave crest, the net accelerates but its motion is restricted by the sinker mass. Thus, a large force is observed in the net which might lead to damage to the net in practice. The behaviour of the system changes for very long waves. Here, both the cylinder and the sinker follow the curvature of the waves, and the occurrence of snap loads becomes less likely. The numerical model quantitatively agrees very well with the experimental data in short and long waves. The snap loads tend to be over-predicted by up to 15%.  Additionally, the change of the motion of the moored cylinder due to the net with sinker is presented in Fig. 21. It shows that the amplitude responses in both, heave and surge, reduce. For the surge motion, the qualitative behaviour of the system is less influenced but smaller amplitudes can be expected. This might be caused by the inertia and drag of the net. In comparison, the characteristics of the system change significantly in the heave direction. The reason is the relatively heavy sinker weight which requires large excitation forces to accelerate in the fluid. Hence, it constrains the cylinder motion for small and medium wavelengths resulting in small amplitude responses.

Engineering applications
The functionality of the proposed numerical framework is elaborated for two types of OOA structures in regular and irregular waves as well as current.

Semi-submersible OOA structure
Semi-submersible offshore fish cages are characterised by being attached to a pre-tensioned mooring system which holds the structure in place but also influences the seakeeping properties of the system. Therefore, the accurate prediction of the dynamic responses to varies sea states is of importance and exemplarily investigated for a structure in the style of Ocean Farm 1. The structure was originally developed by SalMar, Norway and later replicated by The structure is assumed to be rigid and consists of 8 pontoons with straight columns attached. The columns are connected at three different heights via additional thinner columns so that a hexadecagon with a diameter of approximately 1 m is formed. An additional pontoon and column is placed in the middle of the structure slightly below the others and connected to the other columns with thin pipes. The complete structure is shown in Fig. 1. In comparison to the model of Zhao et al. (2019), the thinnest pipes are not considered because of their negligible contribution for the loading and minor influence on the fluid. The draught of the model is 0.28 m in the simulations and reached by adjusting the overall mass of the structure and a free heave decay test. Uniform mass distribution is then assumed to calculate the moments of inertia for the rotational motions of the structure.
A characteristic of this semi-submersible OOA design is that the net is fastened more tightly to the structure than in traditional aquaculture cage systems (Chu et al., 2020). As a consequence, the deformation of the net can be neglected. In the numerical model, this assumption is incorporated by making the net as part of the rigid structure instead of calculating its dynamics separately. The external forces on the net are added to the rigid body solver and used for determining the shading effect as before. The net covers the complete structure and is assembled using a cylinder for the side walls and a cone for the bottom. Each part of the net consists of twines with a length of 8 mm and thickness of 0.6 mm resulting in a solidity of 0.145. A mooring system is attached to the structure for the simulations in waves. The experimental setup includes four mooring lines, each consisting of a rope with a linear spring at the end. The stiffness of the lines is calculated as 195 N/m based on the reported relation between force and line elongation. Further, the pretension is set as 1.91 N. The exact position of mounting points on the structure remains unclear. It is therefore decided to place it on the lower connection columns such that the virtual line extensions intersect in the geometrical centre of the structure. The moored-floating structure is placed in a numerical wave tank of the dimensions 10 m × 2 m × 1.8 m which is a shortened version of the physical wave tank (Fig. 22). The water depth is set to 1 m. The waves are generated at the inlet using the wave relaxation method, and a numerical beach reduces wave reflections at the end of the tank. The simulations include regular waves with height H = 0.06 m, 0.1 m and period T w = 1.0 s, 1.2 s, 1.4 s, which are taken from the experiment. Additional simulations without the net are conducted to study the importance of the net for the motion of the structure. Further, the response of the structure in irregular waves is simulated to gain a deeper understanding of the structural response. Several JONSWAP spectra with a significant wave height of 0.1 m and peak periods between 0.5 s and 3.5 s are chosen for this purpose. The resulting power spectrum is shown in Fig. 25a. Each spectrum is generated by superposing multiple linear wave components as described in (Aggarwal et al., 2019). Power and cross power spectra are calculated using an FFT analysis, and the linear transfer functions (RAO), as well as the coherences γ, are subsequently determined using with S motion and S wave the power spectra from auto-correlation analyses and S wave,motion the cross-spectrum.    Fig. 23. After the first peak, the motions are significantly damped due to the bottom net. The convergence of the results is observed. For the heave motion, the differences are generally small, whereas at least a certain grid resolution is required to capture the pitch motion sufficiently. It is therefore decided to use ∆x = 0.008 m for the simulations including waves. The results of these simulations are shown in Fig. 25 and Fig. 26 as the power spectra and response amplitude operators of the structure and maximum front and aft mooring line forces.
As shown in Fig. 26a, the heave amplitude increases with decreasing wave frequency and increasing wave amplitude. The maximum heave response is expected at f = 0.4 Hz. A second peak, which is indicated in the power spectrum in Fig. 25b, might occur at even lower frequencies. Further, the results indicate a highly damped system as the response to highfrequency excitations is small (Fredriksson et al., 2003). Similar observations can be stated for the surge motion in Fig. 26b. The surge motion increases non-linearly with decreasing wave frequency and approaches values closer to one in very long waves with f < 0.6 Hz as the structure increasingly follows the wave envelope. Also, the regular wave tests reveal larger surge motion for steeper waves due to increased wave energy. Here, the net plays a minor role as the motion without the net shows similar amplitudes. However, the horizontal forces on the net account for about 30% of the total horizontal forces on the system. This indicates that the horizontal forces are generally small, amongst others caused by the low solidity of the net, and that the surge response of the system is mainly dominated by the mooring system. This possible explanation is substantiated by observing the free surface travelling through the structure in Fig. 24. The fluid is accelerated along and wakes are developed behind each member of the structure. In contrast, the damping effect of the net is not visible. For the rotational motion of the structure (Fig. 26c), a strong increase of the pitch amplitude indicates a possible resonance close to the lowest investigated wave frequencies. The response in pitch is relatively small for wave frequencies larger than 0.6 Hz compared to the translational motions. This might be caused by a rather horizontally than vertically acting mooring system. The maximum tension forces in the front and aft mooring lines increase naturally with increased structural motion and reach local maxima close to the maxima of the heave and surge responses. This strengthens the argument that mooring reaction forces are the driving excitation forces for the dynamics of the OOA structure. Generally, the front line forces are larger than the forces in the aft due to the undisturbed impact of the wave loads. The difference between the front and aft forces tends to increase with larger encountered wave periods, whereas the wave steepness mostly affects the aft mooring line as steeper waves travel less disturbed through the upper part of the structure. Further, the simulations without the net reveal that the aft forces tend to decrease if the net is present, as already observed experimentally. This might be caused by the shielding and damping effect of the net. Generally, it is noticed that the results from the regular wave tests mostly coincide with the irregular test results, which indicates that both wave inputs are valid approaches to determine the response of OOA structures.  The obtained transfer functions are based on the assumption that the considered system is linear. The coherences for the motions and tension forces are presented in Fig. 27 to investigate the validity of this assumption. The shown distributions hint at a linear system for wave frequencies between 0.3 Hz and 0.8 Hz because γ is close to unity. The translational motions tend to become non-linear at smaller frequencies than the pitch motion which has a coherence close to one up to f = 1.1 Hz. The strongest non-linear effects are expected for the tension forces which is caused by the coupling to the wave loads and all degrees of freedom.
As a final remark, it is reported that the simulation of the fluid-structure interaction in irregular waves takes around 185 h for 300 s of simulation time on 64 cores (Intel Sandy Bridge) with 2.6 Ghz and 2 GB memory per core.

Mobile floating OOA structure
The offshore aquaculture facility Havfarm 2 (Fig. 28) is developed by Nordlaks and NSK Ship Design (Nordlaks). It can be considered as an example for a mobile floating OOA structure. The main structure is represented by a large, slender ship-shaped hull with several net cages attached. The design process faces challenges due to the complex interaction of multiple nets with the fluid and the resulting water quality change in the cages. The quality is expected to improve with increasing discharge through each net. For this purpose, Havfarm 2 is equipped with a dynamic positioning (DP) system which can change the heading angle between incident flow and structure. However, this increases the external loads which have to be withstood by all components involved. The DP system will be further applied to vary the location of the farm, dependent on the sea state and weather forecast. Hence, the manoeuvrability of the farm and thus, the prediction of global forces is important for the operation of Havfarm 2.
Model tests were performed in the ocean basin of SINTEF Ocean, Trondheim, Norway, in a 1 : 40 model scale to investigate the fluid-structure interaction experimentally. Amongst others, towing tests with different heading angles between the structure and towing direction were conducted and are taken as a reference here.
The prototype of Havfarm 2 consists of multiple rectangular beams forming four equally sized box-shaped spaces. In each of these, a cylindrical net with solidity 0.22 is tightly fastened to the rigid structure. They are simulated as non-deforming nets moving with the rigid structure as explained above. Further, a flexible conical net with the same solidity is attached to the bottom of each cylinder, which requires dynamic modelling. A sinker weight of 2 kg is pre-tensioning this part of the net during the towing tests. Prior to this, free decay tests in heave and pitch are conducted to choose the grid size for the subsequent simulations. Three different grids with uniform cell sizes of 0.07 m, 0.05 m and 0.03 m are considered. Fig. 29 shows the predicted time series in comparison to the experimental results. For the heave motion, the coarser grids tend to over-damp whereas the finest grid predicts most peaks and the period sufficiently. Also, the refinement of the grid results in increased pitch amplitudes close to the experimental values. The period is captured well for the first two peaks and over-predicted for the remaining peaks. A further grid refinement could be necessary to improve the pitch motion and obtain a grid independent solution. It is further noticed that small differences between the numerical and experimental geometry exist because of the neglect of very thin bracings. Their effect on the motion of the structure is, however, assumed to be of minor importance.
Based on the obtained results, a cell size of 0.03 m is chosen for a box around the structure. The box is placed in the middle of a tank of the dimensions 30 × 20 × 10 m with a smooth growth of the cell size towards the boundaries. The water depth is chosen as 8.0 m to avoid interaction with the bottom of the domain. Constant inflows of 0.83 m/s and 1.0 m/s are predefined at the inlet and the flow freely leaves the tank at the outlet. The structure is rotated relative to the inlet with heading angles of 0 • , 15 • and 45 • . The structure can freely heave, roll and pitch during the simulations. However, the motions are not investigated further because the heave motion remains small and the rotational motions are below 1.0 • for all cases.  The mean velocities inside each net, 0.1 m below the free surface, are computed and compared to the experimental data in Fig. 30. The measurements showed large oscillations for which reason the variations are included in terms of one standard deviation. Additionally, slices of the x-y plane around the structure at this height are shown in Fig. 31 to reach a further understanding of the results. It is at first noticed that the computed x-velocities are mostly within the chosen interval of the experimental results. The values also coincide well with the theoretical formula for the velocity reduction through net panels by Løland (1991). This formula is however limited to the case of α = 0 • . For the front cage and small heading angles, the simulations predict an accelerated flow resulting in velocities higher than the inflow velocities. This is also visible in the Figs. 31a and 31b, where the two vertical beams in the front form a narrow channel passed by the flow. At these angles, the shading effects of the nets result in decreasing velocities in the cages behind. The flow separation at each beam becomes increasingly important for the flow field in each cage with increasing angles (see Fig. 31c). At large heading angles, the interaction between the cages becomes less significant. Thus, the differences between the predicted velocities in the different cages as well as the x-velocities itself become small as shown both experimentally and numerically. The mean y-velocities are generally smaller than the mean x-velocities, and the fluid oscillates more in this direction particularly for the front cage and large heading angles. A possible physical explanation is the development of an oscillating wake behind each beam. This causes the  recirculation zones passing the probe points in the centre of each cage periodically. The analysis of the velocity inside the upper part of the cages reveals that the intended improvement of water quality through increased discharge cannot be achieved by increasing the heading angle. However, this changes for the lower, flexible part of the cage. Fig. 32 shows the velocity distribution in the centre x-z plane along the longitudinal axis of the structure for different angles of attack α. As expected, the shading effect of the nets causes the increasing deceleration of the flow along the structure for α = 0 • (see Fig. 32a). Thus, a lower discharge and less deformation are predicted for the cages in the back. By increasing the heading angle (Figs. 32b and 32c), the flow in front of each cage is less disturbed by the wake of cages placed in the front. As a result, similar deformation and discharge are predicted for all cages which consequentially indicates improved water quality in the back cages.
The water quality control through the rotation of the structure has the drawback of increased loads. In order to quantify this, Fig. 33 presents the loads on the net and the structure in x-and y-direction for the different cases. In general, the forces increase with increasing inflow velocity. Also, the loads on the nets are generally more crucial to consider in x-than in y-direction and are even dominant at small heading angles in x-direction. The increase of the loads on the nets is further less dependent on the heading angle than the  rigid structure forces. This is caused by the symmetry of the cage geometry in comparison to the changing structural area exposed to the undisturbed inflow. Thus, the structural forces become the dominant factor for large heading angles. The same conclusion can be drawn in y-direction (Fig. 33b). The increase for larger heading angles is caused by the increased area but might also be influenced by intensified vortex shedding.

Conclusions
A new numerical framework for modelling the motion of OOA structures in waves and current was proposed and applied in this paper. It enables the study of the effects of waves and current on the motion of the system taking into account the fluid-structure interaction around and inside the cages, the motion of the rigid structure as well as the deformation of the net. The interactions of floating structures, mooring, nets and fluid are incorporated as two-way coupling problems. Efficient numerical approaches were chosen and newly developed to solve the governing equations of fluid and structural dynamics. In particular, the coupled solution of the rigid structural motion and the fluid flow was a necessary step to meet the requirements arising from the transition from traditional fish cages towards structures suitable for offshore environments. Several validation cases including rigid and deforming nets, moored-floating objects and a combination of these were presented. Reasonable agreements with available experimental data could be presented for all considered tests and deviations were justified on a physical level. The response of a semi-submersible OOA structure in regular and irregular waves was investigated thoroughly. These structures typically contain rigid nets which implies no volume reduction of the net during operation. For the considered design concept, relatively small vertical motion and strong motion reduction for high-frequency excitation could be observed. This is a characteristic of offshore semi-submersible platforms due to their low centre of gravity, a relatively large mass and the mooring system. Further, the importance of incorporating the net into the investigation increases with the wave height and period due to the increased wave energy and the non-linear growth of the drag forces on the net. The shading effect of the net seems to play a minor role in wave-only cases.
Another application concentrated on the flow around a mobile floating OOA structure in steady current flow. This type of structure can freely rotate and move to different sight locations, which are advantages over semi-submersible structures. The numerical study reveals that the considered structural design results in complex flow patterns with separation and recirculation zones interacting with the upper part of the cages. This complicates the proper adjustment of the discharge through the cages by changing the heading angle. Further, the lower flexible parts of the cages show partly large deformations whose influence on the biomass in the cage has to be taken into account.
It is finally noticed that the predicted small y-loads on the nets for the floating structure might be caused by neglecting the effect of them on the turbulence in the fluid. Both, the potential vortex shedding behind each twine and the change of the fluid vorticity while passing the net can eventually result in increased cross-flow and thus cross-flow forces. This might also change the importance of the net shading for the investigated wave only cases. The quantification of these effects and their inclusion in the numerical framework are left for further research. In future studies the presented framework will be applied to alternative OOA concepts with focus on structural responses under extreme loads.
B Derivation of the left hand side in (31) to the implicit system (32) The dynamic equilibria (21) have to be fulfilled at each knot x at any time. This can be ensured by replacing x (n+1) in (31) with its accelerations using high-order backward finite differences. The weights of each time instance included in the difference are found from (Fornberg, 1998) because of variable time steps in the coupled simulations. Thus, the velocity of the knot is expressed as with c p the weights of the P points of the interpolation. The unknown velocity vectors v (n+1) are approximated by repeating the derivation: Inserting (48) in (47), the left hand side in (31) can be explicitly calculated as: with the definitions After rearranging this result, the implicit function (32) arises in a straightforward manner.