A non-linear implicit approach for modelling the dynamics of porous tensile structures interacting with fluids

A new model for the simulation of large motions of porous tensile structures and their interaction with the surrounding fluid is developed in this paper. The discrete structure is represented by several non-linear elastic bars and knots connecting up to four bars. An implicit system of equations is derived from the fundamental relations of dynamics, kinematics and material and solved using an improved Newton’s method. The Navier-Stokes equations are solved in a numerical domain to account for the interaction with the fluid. The presence of the porous structure is respected in these equations through an additional forcing term based on a modified Lagrangian-Eulerian coupling algorithm. Here, the forces on the structure are distributed on multiple Lagrangian points embedded in the fluid domain. Integration over a suitable Kernel function is applied to distribute these forces on the surrounding fluid. The derived numerical model is suitable for simulating the interaction of porous tensile structures of arbitrary geometry, non-linear material 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 rely on explicit time integration. The validation process shows excellent agreement between the numerical simulations and existing experimental data and demonstrates the applicability of the new methodology for a wide range of applications.


Introduction
Offshore aquaculture has seen growing interest recently because of increasing size of the sites and greater concern over traditional aquaculture due to their environmental impact on coastal regions. The change of environment significantly increases the importance of the accurate prediction of the expected loads and the structural response due to an increased fluid-structure interaction (FSI). A major part of the structure is covered by flexible membranes. These are characterised as spatially intrinsically two-dimensional structures with tensile stress resistance but neglectable bending stiffness. In the case of aquaculture cages, the membranes have high porosity and two dominant stress directions. In the past, segregated approaches considered the motion of these nets without incorporating interaction with the fluid (Løland, 1991) or assuming the validity of potential theory (Kristiansen and Faltinsen, 2015). These studies cannot be regarded as appropriate for offshore conditions due to the non-linearly increasing importance of the FSI for the accurate prediction of the structural and fluid dynamics. In contrast, numerical simulations using computational fluid dynamics (CFD) can be applied to understand the structural and environmental challenges by modelling the forces on and the fluid dynamics in and around the cage. Amongst others, Lewandowski and Pichot (2007) simulated the flow around and inside a rigid net using the Reynolds-averaged Navier-Stokes equations without considering the net motion.
The significant length scale difference between the flow around the whole structure and the flow through each of its voids prevents the resolution of the complete porous structure within the discrete fluid domain. Therefore, a more elaborated approach separates the calculation of the structural and fluid dynamics while respecting their interaction by an appropriate coupling algorithm. The most dominant coupling algorithm is based on a porous medium representation of the porous sheet. For this purpose, a porous volume is defined around the membrane and the governing volume-and Reynolds-averaged Navier-Stokes equations are solved using a finite volume method (Bi et al., 2014a;Patursson et al., 2010;Zhao et al., 2014). Recently, Christensen (2016, 2017) improved the idea by utilising a less restricted model for the proper definition of the porous resistance coefficients. Martin et al. (2020) analysed this approach and revealed several limitations which prevent the applicability of the porous volume analogy for arbitrary shapes and large deformations. They overcame this issue by proposing a new coupling model based on Lagrangian-Eulerian considerations which are generally efficient and broadly applicable. The main idea is built on the developments of Ryzhakov and Oñat (2017) for closed membranes and Aquelet and Wang (2007) for airparachute interaction. Here, disturbances from the solid parts on the surrounding fluid are distributed using Lagrangian points. In comparison, the porosity of the structure led to several modifications of the original approach and the usage of the screen force model of Kristiansen and Faltinsen (2012) to approximate the hydrodynamic forces from the fluid on the porous membrane. The coupling model of Martin et al. (2020) is improved and extended for dynamic problems within the presented work.
Several approaches for modelling the dynamics of tensile structures were presented in the past. Tsukrov et al. (2003) established a finite element modelling using perpendicular onedimensional two-node elements with three degrees of freedom. Lader and Fredheim (2006) introduced the lumped mass method which represents the discrete structure as massless bars connected by mass knots. The solution of the dynamics of the knots was found in terms of their acceleration from Newton's second law, and Runge-Kutta time integration was applied to calculate the knot velocities and positions from the accelerations. The constitutive equations are not automatically satisfied as no iterations are performed which can lead to stability issues and very small time steps in comparison to the fluid solver. The time step restriction is also necessary due to the explicit time integration. A minor modified version of the original approach was successfully coupled to a porous medium model to simulate flexible porous sheets (Bi et al., 2014a) and cylinders (Bi et al., 2014b) in steady flow conditions. In order to overcome the issue of small time steps, implicit methods were proposed. The development of an implicit quasi-static net model was presented in (Martin et al., 2018). The missing time step restriction increases the efficiency of the model, but the approach lacks justification for applications including large motions and snap loads. LeBris and Marichal (1998) introduced an implicit dynamic net model based on the satisfaction of the kinematic relation between knot position and bar length. The original approach was based on inelastic material which leads to very high condition numbers due to missing elements on the main diagonal of the system matrix. Vincent (1999) 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). Thus, the need for an implicit non-linear dynamic net model arises to accurately model the motion of these distinct types of porous tensile structures. In this paper, such an approach is presented by taking Newton's second law as the basis. The external forces due to gravity, inertia and drag are calculated by extending the idea of the screen force model. A non-linear system of equations for the unknown tension forces is derived using highorder finite differences. The system is solved using an improved Newton's method leading to convergence within three to ten iterations for the considered validation cases.
The emerging numerical FSI algorithm is coupled to a numerical wave tank modelling the transport of the interface between two phases. Hence, the model is suitable for simulating the interaction of porous tensile structures of arbitrary geometry, non-linear material 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 rely on explicit time integration.
The remaining paper is structured as follows. Section 2 presents the derivation of the new implicit non-linear structural model, whereas section 3 provides details about the numerical fluid solver and the coupling algorithm. In section 4, the model is verified, and sections 5 and 6 are devoted to several validation cases of porous rigid and flexible sheets and cages in varying wave and current conditions. A possible application of the proposed model is presented in section 7. Conclusions arising from the previous sections are given in section 8.

Non-linear dynamic numerical model for tensile structures
The considered tensile structure is assumed to consists of a large number of square or rhombic meshes forming a porous cylinder or sheet with two distinct stress directions. Thus, an appropriate discrete representation of the structure consists of multiple knots and connecting elastic bars in the principal directions of the meshes. The definition of so-called macro elements is needed for the coupling to the fluid dynamic solution. Each macro element contains four knots and four bars and, depending on the porosity of the structure, represents multiple physical meshes. Instead of the porosity, the solidity S is considered below. The solidity of a porous sheet is defined as the ratio of solid front area to the total area and therefore, equals one minus the porosity. Using the information about the physical twine diameter d t and twine length l t , S is calculated from (Fredheim, 2005) Each structural element is further split into four screens as shown in Fig. 1. Thus, each knot is related to up to four screens. As the screens are not explicitly resolved in the fluid grid, an approximation of the forces on the structure is introduced. Following the assumptions of Morison et al. (1950) for hydrodynamic transparent structures, the mass of the S i screens is lumped at knot x i so that the total mass at that knot m i can be approximated from with m air,s the mass of the screen in air and n s the unit normal vector of the screen pointing in relative velocity direction. The added mass m a,s is approximated as the mass of the water volume occupied by the screen under the assumption that the net is a mesh of multiple cylinders only. Here, the added mass is only applied in the normal direction of the structure.
Figure 1: Discrete structure representation: the thin lines show the physical structure; the structural (macro) elements are represented by black dots connected with thick black lines; the dashed lines present the discretisation into several screens. The red area is related to the red dot for the force calculations.
The dynamic equilibrium equations are formulated for each x i with N i neighbouring knots under consideration of Fig. 2: Here, G i represents the sum of the static gravity and buoyancy force and H i the external hydrodynamic forces consisting of inertia forces I i due to the fluid acceleration a f , and drag and lift forces D i whose calculation are described below. The tension force vectors T ij are defined as the unknown tension force magnitude T ij times the unit vector b ij of the corresponding bar Assuming non-linear material, a constitutive equation is formulated as with l ij the current length of the bar and l 0,ij its original length. Lader and Fredheim (2006) found this relation to be valid with C 1 = 1160 N and C 2 = 37300 N for nylon nets with squared meshes, which is used in the validation section below. Eq. (6) is reformulated for l ij by eliminating the non-physical solution: x z y x j At each new time step (n+1), the unknown position of the knots is related to the unknown tension forces using the kinematic compatibility equation Inserting the constitutive relation (7) in (8) yields The fulfilment of the dynamic equilibria (3) is ensured by replacing x (n+1) in (9) with accelerations using high-order backward finite differences in time. The weight of each point included in the approximation is found iteratively using the algorithm of Fornberg (1998) because of variable time steps in the coupled simulations. In doing so, the first derivative of x (n+1) is expressed as and then reformulated for the position vector: The now arising unknown velocity vectors v (n+1) are approximated using the same procedure, so that In the scope of the paper, third-order accurate polynomials are chosen by setting P = 3. By inserting (12) in (11), the left-hand side in (9) is approximated as with the definitions Under consideration of definition (3), the substitution of (13) in (9) yields a non-linear function F for each bar b ij : with T the global vector of tension force magnitudes and The solution for T is found from (16) using a Newton's method. The improved iterative procedure is implemented for this purpose. Here, F represents the vector of the functions (16) and J is its Jacobian matrix. It is shown by Chun (2005) that (18) converges with third-order. The derivatives of (16) are found separately for the joint bar b ij and its adjoint bars (see A for the expressions). An approximation of the solution has to be given as initial condition. It is proposed to perform the first time step with the linearised version of this model. The solution of the arising linear system only depends on an approximation of the initial position of the structure. The derivation of this model follows the developments of Marichal (2003) and can be found in B. The pre-processing consists of preparing connectivity matrices and follows thereby the concept described previously (Martin et al., 2018). Once a converged result has been found for the tension forces, acceleration, velocity and position of the knots are found from (3), (11) and (12). The description of the structural model is completed by providing details about the calculation of the velocity related forces D i from the fluid using the screen force model (Kristiansen and Faltinsen, 2012). Considering the inertia system of the fluid, this force can be split into a drag and a lift force component in normal (n d ) and tangential (n l ) direction of the local relative velocity vector u rel,s = u s − v s : The surface integral is approximated by a second-order accurate quadrature rule using the geometrical centre as integration point. The drag and lift force directions are defined as with n s the unit normal vector of the screen pointing in the relative velocity direction. 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 normal vector of the screen The definition of the constants c d,0 and c l, π 4 are given in (Kristiansen and Faltinsen, 2012). The determination of the Fourier coefficients is based on non-linear fitting to experimental data as described by the authors in .
3 Direct coupling of the structural response to the fluid solution 3.1 Numerical model for solving the fluid dynamics The conservation of mass and momentum for incompressible fluids arises in convective form as with u the velocity vector, g the gravity acceleration vector, p the pressure and ν the kinematic viscosity. The system is solved numerically using finite differences on rectilinear grids. The convection term is discretised with the fifth-order accurate weighted essentially non-oscillatory (WENO) scheme of Jiang and Shu (1996) adapted to non-uniform point distances, and the diffusion term is discretised using second-order accurate central differences. A staggered grid approach is chosen to ensure a tight connection of the pressure and velocity field. The solution process follows the incremental pressure-correction algorithm for incompressible flows as proposed by Timmermans et al. (1996). In the predictor step, continuity violating velocities u ( * ) are calculated using pressure gradients from the previous time step: The implicit handling of the diffusion term removes decisive restrictions from the CFL condition (Bihs et al., 2016). A third-order accurate total variation diminishing Runge-Kutta scheme (Shu and Osher, 1988) is applied to advance (26) in time. Here, the time step is determined from the CFL condition. At each sub-step, the projection step solves a Poisson equation for the pressure correction term p corr The total pressure is then found from the obtained pressure increment (Brown et al., 2001) using Subsequently, the predicted velocities are projected onto the space of divergence-free fields so that The algorithm is implemented in the open-source CFD solver REEF3D (Bihs and Kamath, 2017;Bihs et al., 2016). Full parallelisation of the model is provided through an n-halo domain decomposition strategy and the message passing interface (MPI) for inter-processor communication. The solver utilises the fully parallelized BiCGStab algorithm with geometric multigrid preconditioning of the HYPRE library (van der Vorst, 1992) for solving the Poisson equation.
If a free surface is present, defined as the interface between water and air phase, the same set of equations is solved. However, the material properties become space and time dependent which is implicitly described by the zero level set of a smooth signed distance function Φ (Osher and Sethian, 1988). The linear advection equation is solved for propagating Φ in space and time. The fifth-order accurate HJ-WENO scheme of Jiang and Peng (2000) is applied for the spatial discretisation and the temporal discretisation is performed by a Runge-Kutta method. The level set function has to be reinitialized regularly to keep its signed distance property. For this purpose, the reinitialisation equation (Sussman et al., 1994) is solved in pseudo time τ so that Φ converges to a valid solution of the Eikonal equation |∇φ| = 1.

Lagrangian-Eulerian coupling algorithm
The development of a new coupling algorithm for the simulation of fluid dynamics around static porous structures was the subject of previous research . Therefore, only the main concepts and modifications for dynamic calculations are provided in this section.
The main idea is built on a Lagrangian-Eulerian approach where the hydrodynamic forces are calculated from interpolated fluid velocities, and the disturbances of the fluid through the solid parts of the structure are incorporated using forcing terms in the Navier-Stokes equations. A flowchart of the complete FSI algorithm is shown in Fig. 3. An example of the distribution of the Lagrangian points on the discrete structure is shown in Fig. 4. They are defined so that the disturbances are nearly equally distributed over the area with distances similar to the surrounding cell sizes. This is achieved by splitting the macro elements into triangles and refine these according to the stationary grid size in their vicinity. The Lagrangian points are then defined in the geometrical centres of each triangle. A momentum loss vector F is defined at each fluid point x e = (x e , y e , z e ) of the Eulerian grid and added to the predicted velocity before solving the Poisson equation. In comparison to , the vector is calculated using with f (x L ) the hydrodynamic force vector on the screen corresponding to the Lagrangian point x L = (x L , y L , z L ), L e the number of Lagrangian points within a defined Kernel around x e chosen as the interpolation kernel of Peskin (Peskin, 1977) Advance and reinitialise free surface using (30) and (31) Predictor step (26) (32) Solve (27) for pressure correction Correct pressure and velocity field with (28) and (29) n = n +1 The forces f (x L ) are represented by the integral of the external forces from (3) over the The chosen definition (32) removes effectively the free parameter which arose previously due to the usage of an inverse distance weighting. A final remark applies to the necessity of correcting the fluid velocity around the structure within the determination of c d and c l in (19). Based on experimental data, these coefficients are determined as a function of the inflow velocity u ∞ , whereas the given velocity information at the structure in the numerical domain is influenced by the pressure jump ∆p caused by the porosity. To overcome this mismatch, a virtual inflow velocity is approximated using Froude's momentum theory (Carlton, 2019), which is based on the momentum balance in front and behind an infinitesimally thin screen. This theory can be utilised to relate the velocity at the screen u s to ∆p: A detailed derivation can be found in . The pressure jump arises from the disturbances in (34) normal to the screen and is defined as In combination with (35), u ∞ can be approximated from u s by solving the intrinsic equation with a Newton-Raphson method and u ∞ = u s as initial value.

Verification of the FSI solver
First, the proper reproduction of the geometrical properties of the structural material is verified using the experimental setup proposed in (Lader and Fredheim, 2006). A 68 meshes long and 4 meshes wide sheet with d t = 1.8 mm and l t = 16 mm is stretched in longitudinal direction using attached weights. The steady state elongation is measured over increasing forces, and the theoretical result is expected to follow (6). Fig. 5a exemplifies the convergence of the elongation of the numerical structural model presented in section 2. Under-relaxation techniques are used to accelerate the convergence and damp the initial shock of the applied force. As can be seen in Fig. 5b, the numerical model follows the theoretical results as expected. Further, the chosen material coefficients show a good approximation of the experimental results. Additionally, a linear material is tested and shown in the same figure. The results indicate proper representation of the elongation for ε < 0.02 but large deviations from physics for ε > 0.08. In a second step, the convergence of the numerical simulation for the validation case of section 6.1 is shown in Fig. 6  The deformed centre line reaches a steady state if at least 9 × 9 elements are used. The reason for the observed fast convergence is the coupling algorithm's mechanism to distribute the screen forces to the fluid domain by splitting the structural elements according to the cell size of the fluid grid. Hence, the chosen number of elements is mostly related to the accuracy of the structural deformation but less influential for load and momentum loss calculations.
Similar studies were conducted to find appropriate grid and structural element sizes for the following validation, showing only the converged results.    5 Validation of the fluid-structure coupling algorithm As the first step of the validation process, the coupling algorithm of section 3.2 is validated against measurements of rigid porous sheets in steady current and waves. Hence, the structural solver is not used. The deviation ε between the numerical result Φ num and the experimental result Φ exp is calculated using The results are shown in Fig. 8. The drag forces increase with increasing inflow velocity and decreasing angle of attack, whereas the lift forces increase with increasing inflow velocity and show a peak for α = 45 • . The velocity behind the sheet decreases due to the presence of the structure. This effect increases with higher inflow velocity but is not influenced by the angle of attack. The distributions of wake velocities in Fig. 8c show the correct behaviour for the numerical model. The deviations between the simulation and the measurements are presented in Tab. 1. Most errors are well below 10%, and the largest deviations are given for the predicted drag and lift forces at large inflow velocities and angles of attack.

Rigid porous sheet in regular waves
A similar experiment was conducted in (Zhao et al., 2008)  The time series of the numerically calculated global drag forces on the two sheets are presented in Fig. 9. For comparison, the measured forces are indicated in red. Spectral analysis is conducted to quantify the period and amplitude of the forces, and the results are shown in Tab. 2 in comparison to the experimental data. Here, the positive and negative force amplitudes are considered separately due to the asymmetry of Stokes waves. In general, the forces increase with increasing wave height due to larger particle velocities below the wave and increasing solidity due to a larger surface passed by the fluid. A good agreement between simulation and experiment can be stated as the deviations are mostly well below 10%. As the measured forces are just provided over two wave periods, increased uncertainties for the experimental data have to be kept in mind.  (Zhao et al., 2008).

Validation of the complete FSI model
The previous section indicates the proper working of the fluid solver and the coupling algorithm for both, constant inflow conditions and waves. In this section, the complexity of the fluid-structure problem is increased by adding the structural solver and hence, the deformation of the structure. First, a single sheet and a cylindrical structure are investigated in steady inflow conditions. Afterwards, the deformation of a porous sheet in waves is presented.

Deformation of a porous sheet in steady current flow
The deformation of a porous sheet in steady current flow is presented following the experimental setup by Bi et al. (2014a). The sheet has a size of 0.3 m × 0.3 m with solidity S = 0.243. It is numerically represented by 9 × 9 elements. The top is fixed during the experiments and a steel bar with a mass of 73 g in air is attached to the bottom of the structure. In the numerical model, the mass and inertia effects of the bar is added to the lowest row of structural elements. The inertia effects are approximates as drag and added mass forces for a cylinder using Morison's formula (Morison et al., 1950). This procedure also holds for the validation cases below. A slice of the computational domain is shown in Fig. 10. The experimental results for the deformation are extracted from pictures of the whole structure presented in their paper and therefore, prone to larger uncertainties. Based on that, the qualitative comparison in Fig. 11c shows a satisfying performance of the numerical model. For larger velocities, the model tends to predict a larger curvature in the middle part of the sheet. This also affects the calculated global drag forces shown in Fig. 11a. Here, deviations below 10% are shown for the lower range of investigated inflow velocities, but a 20% under-prediction is given for largest velocity. This is probably linked to the slightly different deformation causing larger lift forces but smaller drag forces. Additionally, Fig. 11b shows the distribution of the velocity through and behind the sheet for two inflow velocities. As the flow passes the structure, a velocity drop is visible due to a loss of fluid momentum. The magnitude and position of the velocity reduction are well presented by the proposed model irrespective of the inflow velocity. Exp

Deformation of a porous sheet in regular waves
In the previously considered experiment of Zhao et al. (2008), additional measurements of the deformation of a porous sheet in regular waves were presented. The sheet is 0.78 m wide, 0.6 m high and consists of squared meshes with d t = 0.0018 m and l t = 0.06 m. The solidity ratio is approximated as 0.06. The numerical equivalent is modelled using 8 × 8 elements.
The top is fixed during the experiments and an iron bar with a mass of 82 g in air is attached to the bottom of the structure. The same wave tank of 30.0 m × 2.0 m × 2.0 m with a water depth of 1 m is used (see Fig. 12). The centre of the sheet is at (x, y, z) = (10, 1.0, 0.6) m. The same waves as given above are investigated. Additionally, results for waves with wave height H = 0.15 m and wave periods between 1.1 s and 1.5 s are reported. All waves are numerically modelled using second-order Stokes theory. The computational domain is discretised using a uniform grid with ∆x = 0.04 m.   Fig. 13a, the flow pushes the structure in the positive x-direction. Here, the centre of the sheet deforms first because of the increased inertia of the lower part of the sheet due to the additional weight and higher velocities near the free surface. After the wave crest passes the structure, reduced velocities lead to a flattened profile rotated counter-clockwise. When the wave trough approaches as shown in Fig. 13c, the structure is pulled back and starts to rotate clockwise.
The amplitude and speed of this cycling motion depend on both, wave height and period, as it can also be seen in Fig. 14. It shows the maximum deformation of two probe points P 1 = (10, 1, 0.3) m (bottom position) and P 2 = (10, 1, 0.6) m (centre position) following the structural deformation. The deformation increases with the wave height. For both locations, doubling of the amplitudes is given for almost doubling the wave height. Similarly, the increase of the wave period increases the amplitude of the motion due to longer periods of almost unidirectional incident flow velocities. The centre point generally moves less than the bottom point in the numerical model, whereas the measurements indicate larger motion of the centre point for higher and longer waves. Tab. 3 provides the deviations between numerics and experiments for all considered cases. The motion of the lower part of the sheet is represented well as the error is below 10%. In contrast, the centre part deforms insufficient in the simulations for large waves showing under-predictions of up to 18%. These deviations can be linked to the wave representation or the calculated forces on the structure. Unfortunately, both are not reported in (Zhao et al., 2008) so that further analysis of the errors cannot be given.    6.3 Deformation of a porous cylindrical structure in steady current flow As a final validation case, the deformation of a porous cylinder in steady current flow is considered. The setup and measurements are taken from (Lader and Enerhaug, 2005).  Fig. 15. The domain is discretised using ∆x = 0.08 m. The inflow velocity u ∞ varies between 0.13 m/s and 0.56 m/s. As the results, Fig. 16 shows the global drag and lift forces and volume and area reduction coefficients for the different inflow velocities and attached weights. The reduction coefficients represent the ratio of the volume and area of the deformed structure to the initial structure and are calculated as proposed in (Lader and Enerhaug, 2005) to be consistent with the experiments. Their accurate prediction demands a correct force calculation and proper velocity reduction through the front part of the cylinder since both influence the deformed shape. The global forces on the structure, shown in Fig. 16a and Fig. 16b, increase with increasing inflow velocity, and the influence of the changing weights is only of importance for velocities larger than 0.33 m/s. For smaller velocities, the numerical model agrees well with the experiments due to consistent deviations below 10%. For the largest inflow velocity, the lift force is under-predicted. Here, the simulations show consistent results as the lift forces with the largest additional weight is generally the lowest due to the smallest deformation. In contrast, the experimental data shows the largest lift force for this configuration without providing a physical explanation for this phenomenon. The volume and area reduction (Figs. 16c and 16d) is negligible for inflow velocities smaller than 0.2 m/s. For larger velocities, the numerical model accords well with the experiments by predicting increasing volume and area reduction with increasing velocity. The largest deviation is observed for the predicted area reduction coefficient for velocities between 0.23 m/s and 0.27 m/s. A possible explanation is a slightly different deformation process of the numerical cage in comparison to the physical one. In the experiment, the deformation at this velocity seems to be related to a bend of the cylinder, whereas the cage also deforms through a reduced diameter numerically. However, a large uncertainty associated with accessing the area reduction by tracking only three points has to be considered. As expected, the largest deformation of the cylinder is predicted for the lowest additional weight.  The time series of the forces on the left array of nets are shown in Fig. 18. The global forces in x-and z-direction decrease for nets in the wake of another structure due to the momentum loss through the porous sheets in front. The simulations indicate a slight increase of these forces of the nets in the back if d is reduced. However, they are less influenced by d than the forces in the y-direction. For the bigger distance, no influence between the two arrays can be stated as the forces in the y-direction are close to zero. The decrease of the distance results in an increase of lateral forces due to the direct influence of the accelerated flow around the neighbouring structure.  Figure 19: Deformation of the backside centre lines of the nets. The numbering of the nets is indicated in Fig. 17. The results for d = 3 m (Orig) are shown in black, the results for d = 1.5 m (Narrow) in red.
The changing loading conditions also influence the deformation of the structures as can be seen from Fig. 19. As expected, the reduction of velocity in the wake of the first structure results in less rotation of the nets behind. Further, the deformation increases with decreasing distance between the rows due to the increase of forces. In accordance with these observations, the distribution of the tension forces (Fig. 20) can be explained. The largest tension forces are to expected near the clamping on the top, and the strain reduces for the structures in the back of the array. Generally, the strain is larger on the frontside than on the backside due to the momentum loss of the fluid while passing the net. x

Conclusions
In this paper, a new methodology for modelling the non-linear dynamics of porous tensile structures and their interaction with the surrounding fluid was proposed. An efficient structural model was derived for arbitrary deformations and non-linear material. It bases on solving Newton's second law for the unknown tension forces. A new approach for calculating the fluid loading on the structure was proposed as the structure is not directly resolved in the fluid. Here, fluid properties were interpolated on the structural domain using a Kernel function, and the hydrodynamic forces were approximated using a Morison-type approach. High-order backward finite differences were included to approximate the structural motion. Finally, a single matrix-vector problem arose which was solved using an accelerated Newton's method. In contrast to existing explicit algorithms, the implicit time and deformation handling increase stability and effectively remove strong time step restrictions from the fluid solver.
Two-way coupling was provided by including the loss of fluid momentum, due to passing a porous sheet, in the Navier-Stokes equations as a source term. It was determined from a Kernel integration of the hydrodynamic and body forces over multiple Lagrangian points which follow the structural deformation. This represents an innovative extension of the classical forcing approach for porous and hydrodynamic transparent structures. Hence, its application to couple the fluid with slender elastic structures, such as mooring lines and floating beams, and deforming porous media like vegetation is straightforward.
An extended verification and validation study was presented. The proposed model was validated against existing experiments for rigid and elastic porous sheets and cylinders with varying geometries and solidities in current and regular waves. Deviation bands of less than 10% were regularly achieved which indicate a proper calculation of the loads, the wake velocity field and the structural response. Hence, the chosen kernel for the distribution of the momentum loss vector on the fluid domain is a valid alternative to the previous choice and removes effectively additional parameters from the model. The application to the simulation of the flow around multiple elastic net cages delivers insight into possible applications. Here, the proposed numerical model provides the possibility to study fluid dynamics around and inside the cages as well as the effects of waves and currents on the cage array deformations. The resulting fluid disturbances have not been numerically resolved yet but can now be modelled using the proposed approach. Within future work, the framework will also be extended to floating cage systems with mooring systems attached.

A Analytical expressions for the derivatives of F ij
The derivatives of (16) have to be calculated separately for the bar b ij between knot x i and x j and the adjoint bars b ip and b jp : • Derivative for the joint bar b ij : • Derivative for the bars b ip with i = j: • Derivative for the bars b jp with j = i: The resulting Jacobian is inverted using the partially pivoted LU decomposition of the C++ library Eigen (Guennebaud et al., 2010).

B Derivation of a linear system for solving the structural problem
The derivation of a linearised version of the proposed model is based on the kinematic relation Under the assumption of linear elastic material, the length of the bar between knot x i and x j is written as with l 0,ij the original length and κ the elasticity constant. In each time step, the tension force is subject to an incremental increase so that T (n+1) = T (n) + ∆T . Thus, the right-hand side of (42) can be linearised with the argument of small elasticity (κ 1) and small tension fluctuations (∆T 1): The location and velocity of each knot is approximated using first-order backward finite differences in time. Inserting them into (42) and linearising the left-hand side by neglecting terms of higher-order, yields under consideration of (44) Thus, a linear system of equations arises for the tension forces at the new time step using (5) and proper rearrangement: The resulting system matrix is inverted using the partially pivoted LU decomposition of the C++ library Eigen (Guennebaud et al., 2010).