FSI Modeling of Fish Nets For Offshore Aquaculture A two-way coupled fluid-structure interaction model for indus- trial scale problems Master’s thesis in Applied Mechanics OSKAR TYLÉN Department of Mechanics and Maritime Sciences CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2023 www.chalmers.se www.chalmers.se Master’s Thesis In Applied Mechanics FSI Modeling of Fish Nets For Offshore Aquaculture A two-way coupled fluid-structure interaction model for industrial scale problems OSKAR TYLÉN Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology Gothenburg, Sweden 2023 FSI Modeling of Fish Nets For Offshore Aquaculture A two-way coupled fluid-structure interaction model for industrial scale problems OSKAR TYLÉN © OSKAR TYLÉN, 2023. Supervisor: Joakim Hägglund, Aker Solutions AB Examiner: Lars Davidson, Department of Mechanics and Maritime Sciences Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology SE-412 96 Göteborg Sweden Telephone +46 (0)31 772 1000 Cover: Visualization of FSI model. Colorscale on continuous net surface represents mag- nitude of forces acting on net. Grey scale represents velocity reduction of flow as water passes through net. Rendered in Star-CCM+. Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Department of Mechanics and Maritime Sciences Gothenburg, Sweden 2023 iv FSI Modeling of Fish Nets For Offshore Aquaculture A two-way coupled fluid-structure interaction model for industrial scale problems Master’s Thesis in Applied Mechanics OSKAR TYLÉN Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology Abstract In the development of offshore aquaculture for salmon farming, a number of challenges that are not present in the design of traditional nearshore aquaculture systems arise. On offshore locations, the harsher environment that in a sense is what makes it beneficial to move aquaculture offshore is also what increases the complexity of system analysis. Greater current speeds, wave heights and risk of storms all induce a greater demand for accuracy in numerical analysis of loads acting on systems and the resulting dynamic be- havior. To answer this demand, computational fluid dynamics (CFD) is a great tool that can increase the fidelity of system analysis and the understanding of local flow effects affecting the system. Due to the large range of scales involved, it is however not compu- tationally feasible to include high fidelity resolution of net structures in analysis of full scale aquaculture systems. This gives rise to the need of including simplified models used to analyze net structures in otherwise high fidelity CFD models. In the present master’s thesis work, a two-way coupled fluid-structure interaction model based in commercial software Star-CCM+ and Abaqus is developed. The model is based on analytic force models developed for fish nets, and includes the net motion and large scale flow effects resulting from said force models. The forces and deflections calculated by the model are validated towards experiments on a circular net cage exposed to steady currents. In the validation it is shown that the use of simplified net models within a CFD framework is a viable method for consistent prediction of the forces acting on a net. For large levels of net deflection however, the structural model used fails to accurately predict the deformed shape of the net cage, causing force measurements to deviate from experimental data. In addition to model validation, the importance of certain model parameters and components is tested. In this process it is shown that for the force model to remain accurate, it is important to include dependence on net solidity and inflow angle as well as an appropriate reduction of velocity downstream of the net. Keywords: Fluid Structure Interaction, FSI, Offshore Aquaculture, Fish Net Modeling v Acknowledgements I would like to express my gratitude to anyone who, in any way, have afforded their time to aid the progress of the present work. At Aker Solutions, I would first and foremost like to thank supervisor Joakim Hägglund for taking part in plenty of productive discussions and showing a tireless interest in the project. The expertise, feedback and guidance provided has been instrumental for the execution of the project. Secondly, I would like to express my gratitude towards Toon de Jonge for much needed help in understanding the complexity of, and the available options for finite element modeling of the problem at hand. At the Department of Mechanics and Maritime Sciences at Chalmers University of Technology, I would like to thank Lars Davidson for taking on the examination of the project. Lastly, I would like to thank those at Aker Solutions AB who provided me with the opportunity to write this thesis and the resources required to do so. Oskar Tylén, Gothenburg, May 2023 vii Nomenclature Abbreviations CFD Computational fluid dynamics FE Finite Element FEM Finite Element Method FSI Fluid-structure interaction FVM Finite volume method GHG Green house gas RANS Reynolds Averaged Navier-Stokes Re Reynolds number RRMSE Relative root mean square error SQP Sequential quadratic program- ming Greek Variables α Inflow angle α Mass damping coefficient β Stiffness damping coefficient ϵ Strain vector on voigt notation λ Lagrangian multiplier vector σ Lagrangian multiplier vector σ Stress vector on voigt notation λ Lamé’s first parameter µ Dynamic viscosity ν Kinematic viscosity ν Poisson’s ratio νt Turbulent viscosity ρ Density θ Adjacent inflow angle Indices ∞ Denotes free stream value * Optimizer [0] Initial value circ.cyl Circular cylinder D Denotes drag force direction int Interpolated KF Kristiansen & Faltinsen (2012) L Denotes lift force direction N Denotes normal force direction nn Nearest-neighbor r Denotes quantity measured rela- tive to net T Denotes tangential force direction T Transposed t Target Notation ϕ̄ Denotes time average of quantity ϕ ϕ̈ Denotes 2nd order time-derivative of quantity ϕ ϕ̇ Denotes 1st order time-derivative of quantity ϕ â Denotes unit vector in direction of a Roman Variables C Damping matrix c Centroid coordinates D Constitutive matrix F Net force vector f FEM forces g Inequality constraints h Equality constraints K Stiffness matrix L Lagrangian M Modal mass matrix n Net surface normal vector u Fluid velocity vector x Net surface position vector ∆xsource Threshold distance between tar- get and source in interpolation across domains A Area a Fourier coefficient b Fourier coefficient ix Nomenclature C Force coefficient CM Added mass coefficient CRF Current reduction factor dt Twine diameter E Young’s modulus G Shear modulus k Material shear stiffness factor KRe Re scaling factor L Net mesh width p Pressure P k Production of turbulent kinetic energy Re Reynolds number s Source term vector Sn Net solidity t Time V Volume x Contents Nomenclature ix List of Figures xiv List of Tables xvii 1 Introduction 1 1.1 Aquaculture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Salmon Farming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Offshore Aquaculture . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2.1 Modeling Challenges . . . . . . . . . . . . . . . . . . . . . 3 1.2 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Societal, ethical and ecological aspects . . . . . . . . . . . . . . . . . . . . 6 2 Theoretical Basis 7 2.1 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Reynolds Averaged Navier-Stokes . . . . . . . . . . . . . . . . . . . 7 2.1.1.1 Closure Models . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Two-dimensional finite-element formulation . . . . . . . . . . . . . . 8 2.2.2 Material Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2.1 Material Symmetry . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2.2 Rayleigh Damping . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Fish Nets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1 Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1.1 Inflow angle and force directions . . . . . . . . . . . . . . 10 2.3.2 Force Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.2.1 Reynold’s number dependence . . . . . . . . . . . . . . . . 12 2.3.2.2 Solidity dependence . . . . . . . . . . . . . . . . . . . . . 12 2.3.2.3 Angle dependence . . . . . . . . . . . . . . . . . . . . . . 13 2.3.2.4 Added mass force . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.3 Reduced Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.1 Sequential Quadratic Programming . . . . . . . . . . . . . . . . . . 16 3 Computational Model 17 3.1 Modeling of Net Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 xi Contents 3.1.1 Added Mass Force . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.1.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Computational Domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.1 FSI Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.1.1 Mesh Interpolation . . . . . . . . . . . . . . . . . . . . . . 22 3.3 CFD-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.1 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.2 Momentum Equation Source Term . . . . . . . . . . . . . . . . . . 23 3.4 Structural Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4.1 Finite Element Model . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4.2 Loads and Boundary Conditions . . . . . . . . . . . . . . . . . . . . 26 3.4.3 Geometric Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4.4 Material Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4.4.1 Material Damping . . . . . . . . . . . . . . . . . . . . . . 29 3.4.4.2 Compressive stiffness . . . . . . . . . . . . . . . . . . . . . 29 4 Model Validation 30 4.1 Material Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Model Sensitivity Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.1 Value of Model Components . . . . . . . . . . . . . . . . . . . . . . 32 4.2.1.1 Angle Dependence . . . . . . . . . . . . . . . . . . . . . . 33 4.2.2 Parameter Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5 Results 36 5.1 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.1.1 Full Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.1.1.1 Verification of Source Term Effects . . . . . . . . . . . . . 37 5.1.2 Impact of individual model components . . . . . . . . . . . . . . . . 38 5.1.2.1 Orthotropic Material Model . . . . . . . . . . . . . . . . . 39 5.1.2.2 Angle Dependence . . . . . . . . . . . . . . . . . . . . . . 40 5.2 Model Sensitivity Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2.1 Net Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2.2 Element/Cell Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2.3 Material Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6 Discussion 45 6.1 Structural Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6.2 Force Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.2.1 Model Component Importance . . . . . . . . . . . . . . . . . . . . . 46 6.2.2 Angle dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.2.3 Dynamic Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.4 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 References I Appendices III A Validation Results III xii Contents A.1 Løland (1991), orthotropic, one-way coupling . . . . . . . . . . . . . . . . . III A.2 Løland (1991) w/ Re dependence, orthotropic, one-way coupling . . . . . . IV A.3 Løland (1991) w/ Re dependence, orthotropic, two-way coupling . . . . . . V A.4 Løland (1991) w/ Re dependence, isotropic, two-way coupling . . . . . . . VI B Sensitivity Studies VII B.1 Solidity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII B.2 Twine Diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII B.3 Material Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IX B.4 CFD Mesh Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X B.5 Finite Element Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XI xiii List of Figures 1.1 Equivalent CO2 emissions per 100 grams of protein for common animalistic protein sources. [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Global production of fish between 1950 and 2018, split into production from capture fisheries and aquaculture. [3] . . . . . . . . . . . . . . . . . . 2 1.3 Conceptual solution procedure of the FSI-model, as defined by project aim 4 2.1 Parameters Sn, dt and L. Net of similar type as that used in Moe-Føre et al. (2016). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Nomenclature describing the two-dimensional rotation of a differential sec- tion of a net panel and the forces acting on the section in relation to the local flow velocity vector u . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Solidity dependencies for different models in literature. Assuming Ccirc.cyl D = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Comparison of the angle dependence obtained using different Fourier coeffi- cients in equation 2.20. Normalized by CD(Re, Sn, θ = 0). CD(Re, Sn, θ = π) = 0. Histogram shows the distribution of θ across a circular net cage at u∞ = 0.5 m/s as obtained from the model developed in the present work. . 14 2.5 Current reduction factor as applied in AquaSim [30] . . . . . . . . . . . . . 15 3.1 Conceptual two-dimensional wake effect . . . . . . . . . . . . . . . . . . . . 17 3.2 Drag coefficient of a cylinder with diameter dt = 2.42 mm (equation 2.14) in the Reynolds number range of the data of Moe-Føre et al. . . . . . . . . 18 3.3 Net surface together with two-dimensional slice of three-dimensional CFD- mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.4 Nearest neighbor interpolation. When Ct is the target centroid, either Cnn or Vnn will be used as the data source, depending on whether the data is stored on vertices or centroids in the CFD code. . . . . . . . . . . . . . . . 22 3.5 Paths for data transfer between solid and fluid mesh. Numbers on paths denote line in algorithm 1 where path is used. . . . . . . . . . . . . . . . . 22 3.6 2D section of 3D hexahedral CFD mesh used in fluid domain . . . . . . . . 23 3.7 Source term cells in CFD model at initial state. Cells in which source term si is added marked in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.8 Main finite element mesh. Nodes with attached weights on bottom cir- cumference and nodes with translational degrees of freedom fixed on top circumference marked in red. . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.9 Fluid loads per unit area exerted on continuous net surface . . . . . . . . . 26 3.10 Conceptual images of the geometrical differences between the true net and the finite element shell geometry . . . . . . . . . . . . . . . . . . . . . . . . 27 3.11 In plane shear deformation of net structure . . . . . . . . . . . . . . . . . . 28 xiv List of Figures 3.12 Cylindrical coordinate system used to define material behavior of circular net cage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.1 Circular net cage in its non-deformed state. Image courtesy of Moe-Føre et al. (2016) [15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Deformation of circular net cage of solidity Sn = 0.19 at current velocities [0, 0.26, 0.5, 0.76, 0.93] m/s. Image courtesy of Moe-Føre et al. (2012) [15]. 30 4.3 Force data obtained from tests carried out by Moe-Føre et al. "N19" denotes a net with Sn = 0.19 and so on [15]. . . . . . . . . . . . . . . . . . . . . . 31 4.4 Side by side deformation comparison of net used by Moe-Føre et al. and FE representation of net in computational model . . . . . . . . . . . . . . 32 4.5 Comparison of experimental results of Zhan et al. and analytical results obtained from equations 3.1 and 3.2 using coefficients presented in table 4.3. 34 4.6 Optimized angle dependence of force model . . . . . . . . . . . . . . . . . . 34 5.1 Horizontal forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, nr. 3 in table 4.2. . . . . . . 36 5.2 Vertical forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, nr. 3 in table 4.2. . . . . . . . . 37 5.3 Deformation of net cage in experiments by Moe-Føre et al. compared to the deformation of the net cage in the full computational model, nr. 3 in table 4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.4 Normalized velocity magnitude |u|/u∞ of flow field surrounding net cage at u∞ = 0.5 m/s. Instantaneous snapshot. . . . . . . . . . . . . . . . . . . 38 5.5 Distribution of angle θ across all finite elements at u∞ = 0.50 m/s . . . . . 38 5.6 Horizontal forces calculated by the computational model whne including different amount of model component. Model numbering as defined by table 4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.7 Comparison of the net cage deformation obtained at u∞ = 0 using different material models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.8 Comparison of the net cage deformation obtained at u∞ = 0.50 m/s using different material models. . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.9 Horizontal forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, using different coefficients for angle dependence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.10 Vertical forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, using different coefficients for angle dependence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.11 Vertical and horizontal forces obtained using default values together with range between Sn ± 15%. . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.12 Effect of FE size on forces measured by FSI model . . . . . . . . . . . . . . 43 5.13 Comparison of the net cage deformation obtained at u∞ = 0.76 m/s using different target element sizes in the FE discretization. . . . . . . . . . . . . 43 5.14 Comparison of the net cage deformation obtained at u∞ = 0.50 m/s using different material stiffness. . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 xv List of Figures 6.1 Deformation phenomena of computational model that are likely caused by lack of physicality in modeling assumptions. u∞ = 0.93 m/s. a) Large scale folding of net surface on upstream side of net cage. b) FSI model net in orange overlayed on image from Moe-Føre et al. . . . . . . . . . . . . . . 45 xvi List of Tables 2.1 Fourier coefficents for equations 2.20-2.21 obtained by Martin et al. (2020) [26] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.1 Constant material parameters . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Order of implementation of FSI model components . . . . . . . . . . . . . 33 4.3 Optimum values of Fourier Coefficients and parallel drag coefficient . . . . 34 4.4 Tested value for each individual parameter in the performed parameter sensitivity study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 xvii 1 Introduction With a growing global population and an urgent need to reduce environmental impact, the search for sustainable sources of food is a great challenge of modern societies. Today, food production accounts for about one fourth of the global anthropogenic greenhouse gas (GHG) emissions [1]. When breaking down the GHG emissions of food production, about 60% can be attributed to the production of traditional protein sources. This includes emissions from livestock, animal feed and land use for livestock [1]. When reducing the total food production emissions, a large part of the effort should hence be directed towards emissions relating to protein production. Figure 1.1: Equivalent CO2 emissions per 100 grams of protein for common animalistic protein sources. [1] As seen in figure 1.1, the traditional protein sources yield vastly different amounts of emissions per unit mass of produced protein. A natural step towards reduction of GHG emissions relating to protein production would therefore be to increase the amount of protein coming from pig, fish, and poultry, while reducing the amount coming from beef, lamb, and mutton. 1.1 Aquaculture In any attempt to increase the amount of protein coming from fish, aquaculture will play an important part. Since 2016, it has been the main global source of fish for human consumption [2]. With the stagnation of wildly caught fish taking place from the 1990s to present date, seen in figure 1.2, it can be expected that the share of produced fish coming from aquaculture will increase dramatically if the ongoing increase in total fish production is to be continued [3]. 1 1. Introduction Figure 1.2: Global production of fish between 1950 and 2018, split into production from capture fisheries and aquaculture. [3] 1.1.1 Salmon Farming In Norway, a vast majority of farmed fish is salmon. Unfortunately, this industry is plagued by a number of issues relating to the health and well being of farmed salmon. This is problematic both from an economic perspective, as the total yield of salmon is reduced, and from an ethical perspective, where the conditions of the fish can easily be questioned. As with any biologically based production process, aquacultural farming is heavily affected by disease. In the case of salmon farming this is mainly due to the presence of salmon lice, which are known to cause increased mortality and decreased growth of farmed salmon [4]. They are a naturally occurring parasite in wild salmon populations, however the growth and spread of sea lice are amplified by the conditions related to salmon farming. These conditions include the concentration and size of salmon, concentration of farms, and flow through speeds of water. It is estimated that on average, 9% of Norwegian salmon farm revenue is lost due to sea lice infestation [5]. Another issue is the oxygen concentration in the water residing within the farm. This is affected by multiple environmental factors such as salinity, salmon concentration and current velocity [6]. As a result of this correlation, it is considered beneficial to position densely populated sea cages in places where currents allow for constant flushing of the system. In traditional systems, operated in relatively stationary masses of water, the oxygen level must instead be controlled through injection of oxygen within the farm. The control and efficiency of this system is however not entirely straightforward to measure, and hence oxygen concentration risks drifting below sufficient levels despite considerable efforts to avoid this. The most commonly proposed solution to these issues, both from regulatory authorities and the fish farming industry itself, is to move aquaculture production to offshore loca- tions. This would not only reduce issues of parasites and low oxygenation, but would also offer plenty of space for expansion of an industry that is currently fighting for space in sheltered, shallow nearshore waters [7]. 1.1.2 Offshore Aquaculture While the prospect of moving aquaculture production from nearshore to offshore loca- tions is promising for solving some of the main issues of today’s aquaculture, it does not come without issues of its own. The high levels of hydrodynamic energy improving the flushing of aquaculture systems will also induce a significant increase of wear and tear, on the aquaculture system and the fish itself. As a result of this, greater engineering effort 2 1. Introduction will be required to ensure that offshore aquaculture systems are designed in a safe and sustainable way [8]. Part of this effort is put into modeling of the response of different im- portant parameters at different conditions. In Norway developers are required to perform numerical analysis of wave, wind, and current effects as well as substantial analysis of the loads acting on aquaculture systems [9]. For this purpose, numerical methods such as computational fluid dynamics (CFD) could be a great asset in modeling full scale aqua- culture systems in a vast range of operating conditions. Chu et al. (2020) lists a number of challenges in offshore aquaculture where numerical modeling could be useful [7]: • Current speeds and wave action At offshore locations, current speeds are higher, and wave heights are generally greater than at nearshore locations. This difference induces greater motions and structural responses of aquaculture systems, leading to greater risk of structural fatigue and failure, as well as injuries to fish. For the purpose of understanding and mitigating the increased risk, CFD can be a very useful tool. • Accidental storm incidence As offshore aquaculture systems would be designed to operate for a large number of years before decommissioning, they would have to withstand storms and weather incidents of unusual magnitude. Hence, the design of these systems would have to; ensure that operators are not placed in hazardous situations from severe weather conditions, and ensure that as conditions return to normal the welfare of the fish has not been compromised. Once again CFD can be a useful tool that allows engineers to evaluate the effects of such conditions on full scale systems at a level of detail that is difficult (and expensive) to achieve through model testing. • Conductive environment for fish welfare For the best level of growth, the farmed fish also demand certain levels of environmental parameters such as dis- solved oxygen, dispersion of feces and feed, and maximum current speeds. The aforementioned possibility of modeling full scale systems in detail makes it possible to understand the impact that different conditions have on these parameters. If needed additional models could be included in CFD codes to model oxygen disso- lution and consumption, particle dispersion etc. 1.1.2.1 Modeling Challenges At present date, the industrial standard for modeling of aquacultural systems is split into two parts. One option is to model the structural and rigid body dynamics of systems using software such as Aquasim and OrcaFlex [10][11]. Using this approach, a systems dynamic response to external conditions is analyzed through simplified models, generally based in Morison’s equation [12]. This allows for analysis of the large scale response of the systems, however local effects of structural components such as pressure drops, vortex shedding etc. are not captured. When these effects are of interest, one can turn to the second approach, CFD. This allows for high fidelity analysis of both the large scale response of systems and local effects of the system on the flow field, generating a solution that couples system dynamics with fluid dynamics. With the modern development of CFD, it has become possible to make accurate predic- tions of the motion of and forces acting on aquacultural systems without omitting too many important geometrical details. In modeling of hydrodynamic loads and motion of fish nets, there is however greater difficulty in conducting high resolution CFD simula- 3 1. Introduction tions at a reasonable computational cost. This is a result of the large range of scales involved, as commercial aquaculture systems can contain thousands of square meters of net that would have to be resolved at a millimeter scale. Because of this problem, previ- ously mentioned software based in simplified descriptions of the fluid loads are currently the only viable commercial options for fish net analysis. The issue with this approach is that in the net analysis, it neglects the possible local effects of structures surrounding the net, potentially omitting the flow effects that cause the largest hydrodynamic loads. Moreover, the CFD analysis omits downstream effects of the netting in the prediction of the dynamic system response. A big step towards accurate prediction of net loads, as well as dynamic behavior of full scale systems including nets, would therefore be a framework in which the simplified net models can be included in high fidelity simulations of large scale aquacultural systems. Establishing such a framework is the aim of the present work. 1.2 Aim The aim of the Master’s thesis project is to develop and validate a two-way coupled fluid-structure interaction (FSI) model for dynamic load prediction on fish nets, based in commercial software StarCCM+ and Abaqus [13][14]. Herein, the expression “two-way coupled” warrants the inclusion of the following phenomena: • The forces acting on an arbitrary net with certain parameters should be calculated using analytical models where the velocity field given by the Star-CCM+ CFD-solver is used as input. • The displacement and motion of the net structure resulting from aforementioned forces should be calculated using finite element (FE) solver Abaqus and fed back to the CFD-solver. • The corresponding force acting on the fluid should be included in the CFD-solution Conceptually, this gives a high level solution procedure as presented in figure 1.3. Solve relevant CFD-equations in Star-CCM+ Calculate forces acting on net Solve for FE-displacements in Abaqus Feed back deformed net geometry to Star-CCM+ Calculate effect of deformed net on flow field Figure 1.3: Conceptual solution procedure of the FSI-model, as defined by project aim Post development, the established framework should be validated using pre-existing ex- perimental data. Test results presented by Moe Føre et al. (2016) will be used to compare the measured forces and deformation obtained from the computational model to experi- 4 1. Introduction mental data [15]. In the tests, a circular net cage is exposed to steady currents ranging from 0.12 m/s to 0.93 m/s. The deformation of the net cage is then presented visually, hence validation of net deflection has been carried out through visual comparison. 1.2.1 Limitations The present work has been bounded by a set of limitations and delimitations, set either by external factors or internal decisions. • The project was carried out between January 16th and June 4th 2023. As a result of this time frame, some non-essential but interesting modeling options have not been considered in depth. Anything that has been omitted due to the time constraint is listed and shortly reflected upon in section 6.4. • For compatibility with existing frameworks, all CFD modeling has been carried out in Star-CCM+, no other options have been considered. Similarly, to maximize compatibility with Star-CCM+, all FE modeling has been executed in Abaqus and no other options have been considered. • In the modeling of the net structure and the forces exerted upon it, the only option considered has been to model the net as one continuous surface with a given solidity. In its turn, this delimitation has entailed the following; for the forces acting on nets only so called screen models have been considered, and for the FE modeling only different types of two-dimensional elements have been considered. • Because of the available experimental data (see chapter 4), the validation of the model will be restricted to steady state current velocities below 1 m/s. Inertial effects will be taken into account by the FSI model, however due to the steady state nature of the data these components of the model will not be validated. 1.2.2 Objectives From the above aims and limitations, one can identify a number of questions to which answers should be presented in the present work. First, relating to model development: 1. What is a suitable empirical model for determining the forces acting on a net subject to dynamic loads? What are necessary inputs to such a model? 2. How is the effect of said forces on the flow field best accounted for in the CFD solution in Star-CCM+? 3. When modeling net deformation, what is an appropriate level of complexity in FE and material modeling to capture relevant phenomena at a reasonable cost? Second, relating to model validation: 4. How accurately does the model predict the net deflection caused by the predicted forces? 5. How accurately does the developed model predict the forces acting on a fish net exposed to loads from steady currents? 5 1. Introduction 1.3 Societal, ethical and ecological aspects As touched upon previously in this chapter, there are societal, ethical and ecological effects to consider in work relating to aquaculture. On a societal level, aquaculture is as seen in figure 1.2 on of the possible methods available to increase the protein production for a continuously growing world population. From an environmental perspective, figure 1.1 shows that increased consumption of fish can play a part in stopping the global GHG emissions from increasing as the population increases. At its core, the new offshore part of the aquaculture industry is an attempt to solve many of the environmental and ethical issues present today. The ecological effects of sea lice resulting from large concentration of salmon in stagnant water could be resolved, the salmon welfare concerns relating to oxygenation could be resolved and general effects on natural ecosystems could be reduced. If aquaculture is to be moved offshore, developers do however have to ensure that this does not induce new issues such as poor fish welfare from large scale farm motion or release of large amounts of fish into natural ecosystems from cage breakage. As explained, the intended grand scheme purpose of this work is to contribute to the development of offshore aquaculture. This contribution mainly relates to the ability of developers to analyze the aforementioned issues in detail. With computational models, the hope is that issues such as sea lice infestation, low oxygenation, and net cage breakage could be predicted in a digital format and resolved before manifesting themselves in the real world. 6 2 Theoretical Basis Within the coming chapter, any pre-existing theoretical basis that is of relevance either to the project methodology or the succeeding discussion is presented. Theory is presented at the level of depth needed for understanding of the remainder of the thesis. 2.1 Computational Fluid Dynamics In CFD, the Navier-Stokes equations are solved using the finite volume method (FVM). For most cases of industrial scale, including those relevant to offshore aquaculture, solving the original equations would however lead to an unreasonably high computational cost. To solve this issue, common practice is to solve modified versions of the Navier-Stokes equations where some aspects of the fluid flow are modeled rather than solved for in full. 2.1.1 Reynolds Averaged Navier-Stokes Most commonly, the Navier-Stokes equations are modified through the assumption that a flow quantity ϕ can be split into a time averaged part ϕ̄ and a fluctuating part ϕ′. This assumption leads to the Reynolds Averaged Navier-Stokes (RANS) equations, consisting of the continuity equation (2.1) and the momentum equations (2.2). As a result of the modifications the RANS equations are less costly to solve than the original equations, however this comes at the price of having to model the turbulent fluctuations, resulting in extra work and loss of information. In the present work, a slight modification is made to the momentum equations where additions have to be made to account for the change in momentum caused by the forces exerted between the net and its surrounding fluid. ∂ūi ∂xi = 0 (2.1) ρ ∂ūi ∂t + ρ ∂ūiūj ∂xj = − ∂p̄ ∂xi + ∂ ∂xj ( µ ∂ūi ∂xj − ρu′ iu ′ j ) (2.2) 2.1.1.1 Closure Models To close the RANS equations, the unknown final term in equation 2.2, known as the Reynolds stress term, has to be modeled. For a majority of problems, the closure is achieved using so called two equation models where Reynolds stresses are modeled by addition of a turbulent viscosity νt to the kinematic viscosity. In these models, transport equations are solved for the isotropic turbulent kinetic energy k and dissipation of tur- bulent kinetic energy ε [16]. In the equation for turbulent kinetic energy, the production 7 2. Theoretical Basis term P k is written as [17] P k = νt ( ∂ūi ∂xj + ∂ūj ∂xi ) ∂ūi ∂xj (2.3) where the production of turbulent kinetic energy is defined using the gradients of the mean flow and the turbulent viscosity. 2.2 Finite Element Method When solving dynamic structural problems using FEM, one considers the dynamic equi- librium equation Mẍ(t) + Kx(t) = f(t) (2.4) where x contains the Lagrangian positions of the finite element nodes, M is the mass matrix, K is the stiffness matrix and f are any external forces affecting the system. 2.2.1 Two-dimensional finite-element formulation For structures that have a thickness that is significantly smaller than its other dimensions, which is the case in the present work, it is common for the full three-dimensional finite element formulation to be reduced to a simplified formulation using two-dimensional el- ements with degrees of freedom in three-dimensional space. When the description fits well, this approximation provides realistic solutions. In the FE-formulation for simplified two-dimensional elements, it can be shown that if the problem includes both in plane strains and bending of the structure’s mid plane, these phenomena are uncoupled [18]. Because of the lack of coupling, any material response to bending can simply be omitted from the problem without resulting in any further assumptions about the in-plane strains. This is useful in modeling of structures that have no considerable bending stiffness, which is the case for the nylon net structures considered in the present work. In Abaqus, this is implemented for two-dimensional structures that deform in three-dimensional space through the distinction between shell elements and membrane elements, where bending response has been omitted in the membrane formulation while being present in the shell formulation [19]. 2.2.2 Material Modeling For linearly elastic materials, the relation between material strains and material response is governed by the generalized Hooke’s law σ = Dϵ (2.5) where the stress vector σ, strain vector ϵ and constitutive matrix D are defined as [18] σ =  σ11 σ22 σ33 σ12 σ13 σ23  , ϵ =  ϵ11 ϵ22 ϵ33 2ϵ12 2ϵ13 2ϵ23  , D =  D11 D12 . . . D16 D21 D22 . . . D26 ... ... . . . ... D61 D62 . . . D66  8 2. Theoretical Basis For a general anisotropic material, with no material symmetry, it can be shown that the constitutive matrix D is symmetric (Dij = Dji). This leaves us with a constitutive relation between σ and ϵ that requires specification of 21 independent coefficients [20]. 2.2.2.1 Material Symmetry If two coordinate systems exist such that the constitutive matrix remains unchanged independently of which of the two coordinate systems is chosen, a symmetry plane is said to exist. Depending on the number of existing symmetry planes, the constitutive matrix can be simplified further from the 21 independent coefficients described above. In the simplest case, where every plane is a symmetry plane, material response can be described as isotropic [18]. This yields the isotropic constitutive matrix Disotropic =  2G + λ λ λ 0 0 0 λ 2G + λ λ 0 0 0 λ λ 2G + λ 0 0 0 0 0 0 G 0 0 0 0 0 0 G 0 0 0 0 0 0 G  (2.6) where Shear modulus G and Lamé parameter λ are defined as G = E 2(1 + ν) , λ = 2νG 1 − 2ν (2.7) where E is the material Young’s modulus and ν is the Poisson’s ratio. If there exists three symmetry planes and the coordinate planes of the chosen coordi- nate system are parallel to these, the material is considered orthotropic. The resulting orthotropic constitutive matrix consists of 9 independent coefficients [18] Dorthotropic =  D11 D12 D13 0 0 0 D12 D22 D23 0 0 0 D13 D23 D33 0 0 0 0 0 0 D44 0 0 0 0 0 D55 0 0 0 0 0 0 D66  (2.8) As with Disotropic, the coefficients of Dorthotropic can be expressed using Young’s moduli and Poisson’s ratios. However in the orthotropic material description, one can specify a separate Young’s modulus for each of the three orthotropic symmetry planes. Similarly the strain in each of the symmetry planes produced by stress in another symmetry plane is described by 6 separate Poisson’s ratios. [20]. 2.2.2.2 Rayleigh Damping When modeling dynamic phenomena using FEM, energy loss within the material has to be taken into account through damping. For this purpose, a damping term is added to equation 2.4. Mẍ(t) + Cẋ(t) + Kx(t) = f(t) (2.9) 9 2. Theoretical Basis With the damping term of equation 2.9, damping is accounted for through viscous damp- ing matrix C. In Abaqus, C is specified using Rayleigh damping, where it is approximated as a linear combination of the modal mass matrix M and the stiffness matrix K. C = αM + βK (2.10) Equation 2.10 results in a relative modal damping that approaches infinity for rigid body modes unless α = 0. Similarly, unless β = 0, the relative modal damping is increasing with increasing eigenfrequencies. [21] 2.3 Fish Nets Throughout the history of capture fisheries, aquaculture and offshore engineering, plenty of analytical models models describing the drag and lift forces acting on fish nets and similar structures have been developed. Models that are considered relevant to the present work, and the related terminology is described within this section. 2.3.1 Terminology In the present work, the technical terminology used when describing fish nets is the one described by Norsk Standard NS 9415:2021 [9]. All nets considered are knot less nylon nets with square meshes, which can be described by a solidity Sn, a twine diameter dt and a mesh width L. 𝑑t𝐴net 𝐴panel 𝐿 Figure 2.1: Parameters Sn, dt and L. Net of similar type as that used in Moe-Føre et al. (2016). The most important parameter in the modeling of net forces is the solidity Sn, defined as the ratio of the projected net area to the total area of a net panel Sn = Anet Apanel (2.11) The twine diameter dt is the thickness of an individual net twine, assuming that the twine is perfectly cylindrical. The mesh length L width of any given opening in the net, including the twine that separates one opening from the next one. 2.3.1.1 Inflow angle and force directions When tracking the motion of a net panel, its position will herein be described using (in addition to coordinates in the global coordinate system) its angle relative to the local flow 10 2. Theoretical Basis velocity vector u, denoted α. As seen in figure 2.2, this is defined as the angle between the fluid velocity relative to the net ur and net panel normal n in the plane defined by the two vectors. In most models used to calculate forces on net panels, the angle dependence is defined using θ, the adjacent angle of α. 𝒖r 𝑭D 𝜃 𝑭N 𝛼 𝑭L 𝑭T Net panel 𝒏 𝑭R Figure 2.2: Nomenclature describing the two-dimensional rotation of a differential section of a net panel and the forces acting on the section in relation to the local flow velocity vector u It is also in this plane that all relevant force vectors are defined. The resultant force acting on a net panel is split into two separate sets of components. Firstly, normal and tangential components where the normal component acts in the direction of the surface normal and the tangential component is perpendicular to to the normal component while being directed at an angle α from the relative velocity vector in the aforementioned plane. Secondly, drag and lift components where the drag component acts in the direction of the relative velocity vector and the lift component is perpendicular to the drag component and at an angle α + π/2 from the tangential force direction in the aforementioned plane . 2.3.2 Force Models In marine offshore hydrodynamics, the hydrodynamic forces acting on cylindrical objects are often obtained from Morison’s equation, describing the force exerted on a differential section dz of a cylinder by an accelerating fluid [12]. The equation reads dF = [ CM ( ρ πD2 4 ) ∂u ∂t + CD ρD 2 u2 ] dz (2.12) where CM and CD are added mass and drag coefficients, u is flow velocity and D is cylinder diameter. In addition to structural components such as pillars and pipes, equation 2.12 is widely used to describe the forces acting on fish nets, where individual net twines can be approximated as cylinders. To allow for simplifications such as lumping together of multiple twines and use of other measurements than cylinder diameters the force acting on an arbitrary body can be expressed as F = ∫ V CMρ ∂u ∂t dV + ∫ A 1 2CDρF̂D|ur|2dA + ∫ A 1 2CLρF̂L|ur|2dA (2.13) where V is the volume of the body and A is the characteristic area used when determining the force coefficients in use. In equation 2.13, drag force of equation 2.12 is complemented 11 2. Theoretical Basis with a lift force (3rd term) requiring its own force coefficient CL. Directions F̂D and F̂L are defined from convention as described by section 2.3.1.1. The force coefficients of equation 2.13 can be obtained from model tests or analytical models. In the present work, coefficients will be obtained from analytical models described in coming sections. Generally, analytical models for drag and lift coefficients on net panels describe the coefficients as functions of some combination of Reynold’s number Re, adjacent inflow angle θ and solidity Sn, all of which are explored below. 2.3.2.1 Reynold’s number dependence The effect of the Reynold’s number on the forces acting on a net panel is normally found in one of two ways, one of which is testing of specific nets. As described above, coeffi- cients will in the present work be obtained from analytical models, hence, testing is not considered. The other option is to analytically describe the forces acting on a net panel as an extension of theory regarding cylinder drag coefficients. Most commonly, the cylinder drag coefficients are obtained from the seventh order polynomial in log10Re in equation 2.14 [22][23]. The polynomial is fitted to the experimental data presented by Goldstein (1965) [24]. Ccirc.cyl D = −78.5 + 255(log10Re) − 328(log10Re)2 + 224(log10Re)3 − 87.9(log10Re)4 + 20.0(log10Re)5 − 2.45(log10Re)6 + 0.125(log10Re)7 (2.14) The Reynold’s number used is obtained by using the twine diameter as the characteris- tic length and the relative velocity magnitude between net and water. In the work by Kristiansen and Faltinsen (2012), there is also a scaling of the velocity to reflect velocity increase required to fulfill conservation of mass as water passes through a net. However, it is not certain whether or not it is reasonable to say that this scaling should be applied when the data used to determine equation 2.14 includes only free stream velocities. This leaves two options for calculation of the Reynolds number: Re = dtUr ν , ReKF = dtUr ν(1 − Sn) (2.15) After obtaining a drag coefficient from equation 2.14, it is scaled to be used in equation 2.13 where the characteristic length is a net panel area rather than a twine diameter. As this scaling is solidity based, it is described below. 2.3.2.2 Solidity dependence As with Reynold’s number, the effect of solidity can be included in multiple ways. One, is to include the solidity when converting cylinder force coefficients to screen force coeffi- cients. Kristiansen and Faltinsen (2012) use two scaling methods, given by equations 2.16 and 2.17. C KF [1] D (Re, Sn, θ = 0) = Sn (1 − Sn)2 Ccirc.cyl D (2.16) C KF [2] D (Re, Sn, θ = 0) = Sn(2 − Sn) 2(1 − Sn)2 Ccirc.cyl D (2.17) It is concluded by the authors (and confirmed by figure 2.3) that 2.16 is the more conser- vative option, resulting in larger net panel forces [22]. 12 2. Theoretical Basis Figure 2.3: Solidity dependencies for different models in literature. Assuming Ccirc.cyl D = 1. An alternative to these formulations is equation 2.18, proposed by Berstad (2012), which is used in Aquasim [23][10]. As seen in figure 2.3, this is the least conservative alternative presented here. CBerstad D (Re, Sn, θ = 0) = Ccirc.cyl D Sn( 1 − Sn 2 )3 (2.18) The other option is to use experimental data for nets of different solidities and fit an analytical model to the experimental data. This was done by Løland (1991), yielding the relation between drag coefficient and solidity given by equation 2.19 [25]. CLøland D (Sn, θ = 0) = 0.33Sn + 6.54Sn2 − 4.88Sn3 CLøland L (Sn, θ = π/4) = −0.05Sn + 2.3Sn2 − 1.76Sn3 (2.19) As one can note, the experimental data used to generate the analytical model includes no Reynold’s number dependence. 2.3.2.3 Angle dependence The dependency on the angle between the net panel and the flow velocity vector (inflow angle) could be considered the most important to capture accurately when modeling net screens subject to large inflow angles. This is commonly achieved by allowing the force coefficients to be represented by their Fourier series, using different values for coefficients ai and bi in equations 2.20-2.21 [26][22][25]. CD(Re, Sn, θ) = CD(Re, Sn, θ = 0) ∞∑ n=1 a2n−1cos ((2n − 1)θ) (2.20) CL(Re, Sn, θ) = CL(Re, Sn, θ = π/4) ∞∑ n=1 b2nsin (2nθ) (2.21) The choice of coefficients (and how many coefficients to include) however, is not straight- forward. In existing literature it is most common to use only a1 = b2 = 1 and set the remaining coefficients to zero [25]. As is made clear by figure 2.4 there is however a risk of over predicting the drag forces acting on circular net cages when using these values. This over prediction risk was shown by Kristiansen and Faltinsen (2012), by comparison with 13 2. Theoretical Basis experiments on a rigid circular net cage carried out by Zhan et al. (2006) [22][27]. They conclude that it is likely that, for their force model, a3 should be in the range [0, 0.1] but do not investigate the issue further. The only attempt to tune the coefficients, known by the author, was carried out by Martin et al. (2020). In this work, coefficients are fitted to the experimental data of Patursson et al. (2010) and Zhan et al. (2006) using the force model of Kristiansen & Faltinsen (2012), this results in the Fourier coefficients presented in table 2.1 [26][27][28]. Table 2.1: Fourier coefficents for equations 2.20-2.21 obtained by Martin et al. (2020) [26] a1 a3 a5 b2 b4 0.9725 0.0139 0.0136 1.2291 0.1116 Figure 2.4: Comparison of the angle dependence obtained using different Fourier coefficients in equation 2.20. Normalized by CD(Re, Sn, θ = 0). CD(Re, Sn, θ = π) = 0. Histogram shows the distribution of θ across a circular net cage at u∞ = 0.5 m/s as obtained from the model developed in the present work. One additional aspect to consider regarding angle dependence is whether or not to assume that the drag force acting on a net panel is zero when it is perfectly parallel to the flow velocity vector. Although it is trivial that there will be some drag force acting on a parallel net panel, there is no clear consensus on whether to include it in the existing literature [22][25]. To allow for the addition of a parallel drag coefficient, equation 2.20 can be modified: CD(Re, Sn, θ) = CD(θ = π/2)+ (−CD(θ = π/2) + CD(Re, Sn, θ = 0)) ∞∑ n=1 a2n−1cos ((2n − 1)θ) (2.22) In most existing literature, CD(θ = π) is set to zero. However Løland (1991) sets CD(θ = π) = 0.04 based on experimental results. 2.3.2.4 Added mass force In dynamic analysis of nets, containing unsteady motion, the effect of acceleration becomes important. This effect is included through the added mass coefficient CM in the 1st term 14 2. Theoretical Basis of equation 2.12. In the present work, the inclusion of added mass effects has however been simplified through the assumption of potential flow. It has been shown that for the motion of an ideal fluid, the added mass of a circular cylinder can be approximated as precisely equal to the displaced fluid mass [29]. This yields mtotal = (ρsolid + ρfluid)Vsolid (2.23) where mtotal is the solid mass used in the equation of motion of the solid. When the corresponding approach is used for an arbitrary body in an FEM-approximation, the addition of the fluid density will cause an alteration of the mass matrix in equation 2.4. 2.3.3 Reduced Velocity In most existing software used to analyze nets, the effect of the net on the fluid, causing a current reduction, is introduced through a reduced velocity model. In these models, a current reduction factor is applied to the flow affecting a net downstream of any other net as shown in figure 2.5. Figure 2.5: Current reduction factor as applied in AquaSim [30] In AquaSim, this factor is applied as suggested by Løland (1991) CRF = 1 − 0.46CD (2.24) where CD is the drag coefficient of the upstream net as defined by equation 2.19 [30]. The drawback of this modeling approach is that it does not account for conservation of mass, which leads to a risk of under estimating the current velocity of flow passing by the sides of a net panel. As net configurations for offshore aquaculture systems become increasingly complex, this could in its turn lead to under estimation of forces acting on net panels downstream of the most upstream net panels in the configuration. 2.4 Optimization In many specific engineering applications of theoretical models such as the net models presented in section 2.3.2, it can be necessary to tune model coefficients to fit the specific application. For this purpose, it can be useful to implement optimization algorithms, using some simplified model to find appropriate coefficients. In this process, the problem 15 2. Theoretical Basis is defined on so called negative null form as generalized by equation 2.25 [31]. min f(x) subject to h(x) = 0 g(x) ≤ 0 (2.25) On negative null form, the problem is defined as a minimization of some objective function f(x), subject to a set of equality constraints h(x) and a set of inequality constraints g(x). Mathematically, the problem can then be solved using a number of different methods. 2.4.1 Sequential Quadratic Programming One of the most efficient general purpose algorithms for solution of the minimization prob- lem is sequential quadratic programming (SQP). The algorithm utilizes the Lagrangian L(x, λ, σ) = f(x) + λTh(x) + σTg(x) (2.26) where λ and σ contain the Lagrangian multipliers for which the constraints in h(x) and g(x) are satisfied. This allows the optimum of f(x) to be found as the point that satisfies ∇L(x∗, λ∗, σ∗) = 0T (2.27) where subscript * denotes optimizers of the function. Equation 2.27 is then solved using Newton’s method, giving the optimum and optimizers of f(x). 16 3 Computational Model In the present work, an FSI model based in Star-CCM+ and Abaqus is developed and subsequently validated. In model development, modeling decisions are therefore not only based upon what is theoretically suitable, but also on availability of models and practical feasibility of implementation in aforementioned software. This decision making process, as well as the reasoning behind and consequences of certain decisions are described within the coming sections. As prescribed in section 1.2.1, only one option has been considered for the geometrical representation of the fish net within the computational model. The decision of represent- ing the net as a continuous surface has some fundamental consequences for the available modeling options. • There is no possibility of resolving the flow effects of individual net twines using the spatial discretization in Star-CCM+. Instead, the effect on the fluid caused by the net will have to be expressed per unit area of the continuous surface, effec- tively resulting in a spatial averaging of wake and pressure effects as is conceptually visualized in figure 3.1. High Pressure Low Pressure Individual net twines Continuous net surface Figure 3.1: Conceptual two-dimensional wake effect • Just as the flow around each individual twine will not be modeled in Star-CCM+, the structural response of each individual twine will not be modeled in Abaqus. Because of this, it is not solely the material response of the polymer making up the twines that has to be modeled, but also the response of the structure in which the twines have been sown together. • In the modeling of net forces, only so-called screen models are practical to implement using geometrical properties of the continuous surface. As a result of this distinction, the drag and lift coefficients of equation 2.13 have to be expressed in a way so that the area across which the integration is carried out is the area of the continuous surface, not the projected area of individual net twines or similar. 17 3. Computational Model 3.1 Modeling of Net Force For the net force modeling, a step-by-step approach has been taken to understand the impact of including different physical phenomena. As a staring point, the net model developed by Løland (1991) was chosen. As this model is based on experimental data of forces acting on net panels, it was chosen over the other models discussed in section 2.3.2, that use experimental data regarding cylinders. However as shown in figure 3.2, the experimental data used for validation includes a fairly wide range of Reynolds numbers within which the drag coefficient can be expected to take a large range of values. Because of this, a Reynolds number dependence based on equations 2.14 and 2.18 was added to the drag and lift coefficients given by Løland (1991). Figure 3.2: Drag coefficient of a cylinder with diameter dt = 2.42 mm (equation 2.14) in the Reynolds number range of the data of Moe-Føre et al. This Reynolds number dependence was added to the force coefficient formulation by addition of a scaling factor KRe = C Berstad (2012) D (Re, Sn, θ = 0) C Løland (1991) D (Sn, θ = 0) to the full expressions for drag and lift coefficients. This yields final force coefficients expressed as CD(Re, Sn, θ) = KRe ( CD(θ = 0)+ ( −CD(θ = 0) + C Løland (1991) D (Sn, θ = 0) ) ∞∑ n=1 a2n−1cos ((2n − 1)θ) ) (3.1) CL(Re, Sn, θ) = KReC Løland (1991) L (Re, Sn, θ = π/4) ∞∑ n=1 b2nsin (2nθ) (3.2) In the main model used in the present work, only one coefficient is used for each of the Fourier series representing the angle dependence (a1 = b2 = 1). For the parallel drag coefficient CD(θ = π/2), Løland’s value of 0.04 is used. This leaves the angle dependence of the original model of Løland (1991) unaltered. However as outlined in section 4.2, the impact of this modeling decision is later investigated by use of different coefficients. Variable Solidity To account for large deformations of the net, the solidity used in equations 3.1 and 3.2 is allowed to change as a function of the deformed net area. This is achieved by simply 18 3. Computational Model assuming that the change in solidity is the inverse of the change in area, yielding Sn = Sn[0] A[0] A (3.3) where subscript [0] denotes the parameter value in the undeformed state of the net. In practice this model assumes that as the area of the continuous surface changes, the total area of net material in each finite element Anet remains constant. 3.1.1 Added Mass Force To properly account for inertial effects, the added mass of displaced fluid is accounted for through modification of solid density as suggested by equation 2.23 yielding ρnet = ρnylon + ρfluid As the co-simulation framework of Star-CCM+ and Abaqus does not allow for export of solid acceleration, this force is however only included in the FE solver, and not in the momentum term added in the CFD-solver. As the present work is validated through steady state experimental results, this is not considered an issue. However, if inclusion of added mass forces in the momentum source term is considered necessary, this could be achieved through manual calculation of the acceleration from the velocity data exported from Abaqus. The added mass force would then be calculated using this acceleration and the estimated net volume in each differential section of the continuous net surface, which together with equation 2.23 would give the extra force that should be added to the momentum source term. 3.1.2 Implementation In Star-CCM+ the previously described model has to be mathematically implemented using the definitions of figure 2.2. At an arbitrary differential section of the net surface, this can be achieved using the local flow velocity u, the local net velocity ẋ and the local surface normal n. The initial step is to determine force directions F̂D, F̂L, F̂N and F̂T using the relative velocity between fluid and net ur = u − ẋ. The directions are obtained from equations 3.4-3.7. F̂D = unit (ur) (3.4) F̂T = − ( unit(ur × n) × n ) (3.5) F̂L = − ( unit(F̂D × F̂T) × F̂D ) (3.6) F̂N = − ( unit(F̂D × F̂T) × F̂T ) (3.7) where unit(a) = a |a| This ensures that the convention defined in section 2.3.1.1 is followed, meaning a) Drag and lift forces are perpendicular b) Normal and tangential forces are perpendicular c) All forces act in the plane defined by ur and n. Additionally it ensures that, if ur acts in the opposite direction of that displayed in figure 2.2, all resulting forces are accordingly 19 3. Computational Model reversed. Together with equations 3.4 and 3.6, equations 3.1 and 3.2 define the total force per unit area acting on the differential section as FR = FD + FL (3.8) with FD = F̂D 1 2CD(Re, Sn, θ)ρ|ur|2dA FL = F̂L 1 2CL(Re, Sn, θ)ρ|ur|2dA (3.9) When needed, e.g. when exporting pressure and wall shear stress to Abaqus as outlined in section 3.4.2, the total force FR can be split into its normal and tangential components by projection onto F̂N and F̂T: FN = FR·F̂N F̂N FT = FR·F̂T F̂T (3.10) 3.2 Computational Domains As described by figure 3.1, the small scale flow effects of the net twines will not be captured by the CFD model as long as the net is represented by a continuous surface. Therefore, there is no need to include the net geometry in the spatial CFD discretization. This allows for a discretization such as that displayed in figure 3.3, where the FE mesh and the CFD mesh exist independently of one another relieving any need to morph the CFD-mesh. CFD mesh Net surface Figure 3.3: Net surface together with two-dimensional slice of three-dimensional CFD-mesh This framework is advantageous for two main reasons: • It avoids the need to recalculate the positions of mesh vertices that are affected by net deformation each time the CFD solver receives a new surface geometry from the FE solver. • There is no risk of affecting the numerical solution obtained from the CFD solver by reduction of the mesh quality. Deforming the CFD-mesh would likely lead to skewed mesh cells, causing risk of numerical diffusion which could affect the final solution. Additionally, at presence of large deflections, there would be a risk of moving vertices far enough to create negative or zero volume cells, leading to numerical singularities. With the benefits of not including the net surface in the spatial discretization, there is however an additional demand for transfer of information between different discretizations. The procedure for this data transfer is described in section 3.2.1.1. 20 3. Computational Model 3.2.1 FSI Coupling When solving the FSI problem, the governing equations of the individual domains are solved sequentially. In practice this means that the relevant equations on one of the do- mains are solved first, and the resulting data is used to subsequently solve the relevant equations on the other domain. As recommended for mechanical FSI coupling, the se- quential solution of the governing equations is lead by solution of the equations on the solid domain [32]. The full coupling algorithm is presented in algorithm 1. Algorithm 1 FSI coupling algorithm 1: while Solution time < End time do 2: for n in range(1, max(1,N )) do ▷ If needed, this allows for N iterations of the FSI coupling within each time step 3: Solve time step in Abaqus 4: Map x, ẋ to Star-CCM+ using co-simulation interface 5: Calculate FR ▷ See equations 3.8 and 3.9 6: for Cell in CFD Mesh do 7: if abs(|xnn| − |c|) < ∆xsource then ▷ Find source term cells 8: Map FR using nearest-neighbor interpolation 9: end if 10: end for 11: Calculate si ▷ See equation 3.12 12: Solve time step in Star-CCM+ 13: if n = N then 14: Map u to net surface using nearest-neighbor interpolation 15: end if 16: Calculate FN, FT ▷ See equations 3.8-3.10 17: Map FN, FT to Abaqus using co-simulation interface 18: end for 19: end while For the purpose of the present work, the FSI coupling runs stably when exchanging infor- mation once per time step. If needed, the co-simulation capabilities of Star-CCM+ and Abaqus nevertheless allow for data exchange on a per-iteration basis. Within the present framework, there are however limited reasons to exchange data once per iteration. As the solid boundary of the net is not present in the CFD discretization, no true equilibrium between solid and fluid forces is solved for. There is therefore no reason to update the fluid velocity used to calculate forces more than once per time step. When the acceler- ation of the solid is large, it can however be useful to update ẋ and re-calculate the net forces more than once per time step. If the magnitude of ẋ is increasing quickly, this will help provide extra fluid damping to the FE solver which will complement the material damping provided by the Rayleigh damping added through viscous damping matrix C in equation 2.9. “Time step” in algorithm 1 denotes the FSI timestep. On lines 3 and 12, Abaqus and Star-CCM+ can solve the FSI time step using multiple internal iterations and/or sub time steps. 21 3. Computational Model 3.2.1.1 Mesh Interpolation The information transfer between the CFD mesh and FE mesh is carried out using Star- CCM+’s built in nearest-neighbor interpolation. In the interpolation scheme, the data received at the target coordinates is the data available at the nearest viable source coor- dinates. Both target and source coordinates can belong to either a mesh vertex, a face centroid or a volume centroid, depending on where the data in question is stored [32]. 𝑉nn 𝐶nn 𝐶t Figure 3.4: Nearest neighbor interpolation. When Ct is the target centroid, either Cnn or Vnn will be used as the data source, depending on whether the data is stored on vertices or centroids in the CFD code. Using this framework, necessary data is transferred between the FE domain and fluid domain inside Star-CCM+. Any data that has to be transferred between Abaqus and Star-CCM+ can then be transferred between the FE domain in Star-CCM+ and the FE domain in Abaqus using the built in co-simulation interface. All data transfer paths used by the FSI coupling algorithm are shown in figure 3.5. Co-simulation interface Nearest-neighbor interpolation Fluid DomainAbaqus Net Surface 𝑭R 𝑭N, 𝑭T 𝒙, ሶ𝒙 𝒖 17 4 8 14 Figure 3.5: Paths for data transfer between solid and fluid mesh. Numbers on paths denote line in algorithm 1 where path is used. 22 3. Computational Model 3.3 CFD-model For the modeling of the fluid, the RANS-equations closed by the k-ω SST model are solved. As described in section 3.3.2, the effect of the net on flow turbulence is disregarded and not properly included in the model. 3.3.1 Spatial Discretization The only geometry present in the domain is the continuous net surface, which as previously explained is not resolved by the CFD mesh. For this reason, the fluid domain can be coarsely discretized using a hexahedral mesh with a uniform cell size of 0.25 m. Figure 3.6: 2D section of 3D hexahedral CFD mesh used in fluid domain However, as one can observe in figure 3.6, a refinement is made in the area surrounding the net (visualized by the transparent grey surface). This refinement is not necessary in the one-way coupled models, however when using two way coupling it is considered useful to acquire a decent resolution of the flow gradients resulting from the source term in equation 3.11. In the refined volume, cell size is uniformly set to 0.0625 m. 3.3.2 Momentum Equation Source Term As the CFD solver solves transport equations for momentum, there is no need to imple- ment a reduced velocity model such as that outlined in section 2.3.3. Instead, a source term can be added to equation 2.2, ensuring conservation of momentum across the inter- face bewteen the CFD and FE domains. ρ ∂ūi ∂t + ρ ∂ūiūj ∂xj = − ∂p̄ ∂xi + ∂ ∂xj ( µ ∂ūi ∂xj − ρu′ iu ′ j ) + si (3.11) The addition of momentum source term si, of dimensionality N/m3, yields equation 3.11, which is used in any fluid cells within a certain distance ∆xsource from the nearest finite element of the continuous surface. An example of this cell selection can be seen in figure 3.7. 23 3. Computational Model (a) Side View (b) Top View Figure 3.7: Source term cells in CFD model at initial state. Cells in which source term si is added marked in red. As the source term is applied, flow velocity will be reduced as water moves past the coordi- nates of the continuous net surface. As the RANS equations are solved, this will however not only effect the downstream flow. In order to fulfill the mass conservation imposed by equation 2.1, the pressure field around the net will adapt and the fluid mass that no longer passes through the net will have to take another path. This will cause an increased velocity in the fluid surrounding the net structure. For complex net configurations, this will be advantageous compared to the reduced velocity models described in section 2.3.3 where fluid mass is essentially removed from the system as velocity is reduced downstream of a net panel. Additionally, as the source term is applied in the opposite direction of the resultant force FR acting on the net, it will not only alter the magnitude of momentum, but also the direction. To ensure correct magnitude and dimensionality of the source term in the Star-CCM+ implementation, a series of steps have to be taken. First, the force calculated on the continuous surface of the FE domain is transferred to the fluid domain as part of the interpolation procedure described in section 3.2.1.1. In each target cell, the force is then converted to N/m3 using the local cell volume. Depending on the size of ∆xsource and the FE and CFD mesh sizes, it is however possible (and likely) that the total force exerted the fluid will not match the total force exerted on the net. To ensure conservation of momentum, the source term magnitude in each cell is therefore normalized using the total volume of the source term cells, the total force transferred through the interpolation procedure and the total force acting on the net, resulting in equality between the final total force acting on the fluid and the total force acting on the net. Using this procedure, the source term vector is defined as si = −f int i ∫∫ AFE ∣∣∣fR i ∣∣∣ dA ∫∫∫ Vsource ∣∣∣f int i ∣∣∣ dV (3.12) Where Vsource is the total volume of the interpolation target cells, f int i is the interpolated force per unit volume, AFE is the total area of the continuous surface and fR i is the resultant force per unit area as defined by equation 3.8. In the developed model, the effect of the net on the fluid is only considered through the momentum equations. No consideration is given to the turbulence generated by the presence of the net. If the net geometry was resolved in detail, near wall gradients would generate turbulent kinetic energy in the k-ω SST model through the production term of equation 2.3. In the present model, turbulent kinetic energy will only be generated from 24 3. Computational Model large scale gradients resulting from the transport of momentum downstream of the source term cells. The consequences of this modeling decision are not investigated. 3.4 Structural Modeling As briefly described at the beginning of this chapter, the representation of a net as a continuous surface is a simplification that creates a need to model not only the behavior of the net material, but also the behavior of the net structure. Additionally, any volume loads acting in the net have to be appropriately scaled to account for the volume difference between the true net geometry and its FE counterpart. This modeling is mainly carried out through alterations to material models and material properties, and in part through the FE formulation. 3.4.1 Finite Element Model As described in section 2.2.1, the small ratios between a nets thickness and its other dimensions make it appropriate to use two-dimensional finite elements. In addition, this formulation simplifies the information transfer between different computational domains as the net surface used by the screen force model (figure 3.3) will share the geometry of the mid-plane surface in the finite element approximation. In Abaqus, the net surface is defined and discretized using finite elements of a constant size, which are given a constant thickness equal to the net twine diameter. As we will see in Chapter 5, the minuscule thickness together with small material stiffness leads to the structure having a close to negligible bending stiffness despite using a shell formulation. For this reason, shell elements were chosen over membrane elements as the FE model was significantly harder to converge using the latter option. For all results presented in chapter 5, the two-dimensional surface has been discretized using 4-noded quadrilateral shell elements with reduced integration. The target element size is set to 0.05 m, resulting in appproximately 30 elements in the z-direction and 110 elements in the θ-direction. Weights attached to the net are approximated as point masses in 16 uniformly distributed FE nodes along the bottom boundary. This results in the FE-discretization displayed in figure 3.8. Figure 3.8: Main finite element mesh. Nodes with attached weights on bottom circumference and nodes with translational degrees of freedom fixed on top circumference marked in red. The dynamic analysis of the problem was carried out using Abaqus’ implicit solver pro- 25 3. Computational Model cedure. As large deformations such as those displayed in figures 4.1 and 4.2 should be expected, the solver for non linear geometry is used. 3.4.2 Loads and Boundary Conditions To model the net attachment to the steel ring seen in figure 4.1, all finite element nodes marked by the red circle at the top boundary in figure 3.8 were fixed in space from translating motion. However, to capture the small shear and bending stiffness of the net the rotational degrees of freedom of these nodes were not fixed. In total the structure was subject to 4 separate loads, two body loads and two external fluid loads. Gravitational load was applied by giving the continuous surface a density ρFE using geometric scaling from the net density ρnet as described in section 3.4.3. To properly capture the inertial effects of the net mass, this density was not adjusted for buoyancy. Instead, an additional body load was added to include buoyancy effects. This load was defined as the sum of the displaced weight of water and the weight added from the added mass formulation described in section 3.1, acting in the opposite direction of the gravity vector. With this combination of body loads, it is ensured that the inertia, buoyancy and gravitational loads acting on the continuous surface are correctly modeled to best represent the true net geometry. The mass of the weight cylinders was accounted for by adding their submerged mass in the aforementioned point masses. As the replicated experiments measure a steady state condition, the inertial effect of the difference between submerged and non-submerged cylinder mass was not considered of importance. Figure 3.9: Fluid loads per unit area exerted on continuous net surface The external fluid loads were applied as a pressure and a surface traction acting on one of the two sides of the continuous surface. As the pressure is applied normal to the surface and the surface traction is applied tangentially to the surface, the drag and lift forces calculated from equation 3.9 are projected onto the correct directional vectors using equations 3.8 and 3.10. The normal and tangential forces are then imported into Abaqus from Star-CCM+ using the co-simulation interface and applied as a pressure and surface traction respectively as visualized by figure 3.9. 3.4.3 Geometric Scaling To account for the geometric differences between the true net geometry and its surface rep- resentation, certain properties are directly scaled based on the volume difference between 26 3. Computational Model the two. 𝐴net 𝐴FE (a) front view 𝐴net 𝐴net 𝐴FE 𝐴FE 𝑑t 𝑑t 𝑑t 𝑑t (b) Side/Top view Figure 3.10: Conceptual images of the geometrical differences between the true net and the finite element shell geometry The appropriate scaling factor is determined from the projected areas of the true net and the FE-model. From the front, as displayed in figure 3.10a, it is natural to scale using the solidity Sn (see equation 2.11). From the side/top, as displayed in figure 3.10b, one has to make an assumption about the cross-sectional shape of the true net. In the present work, the net is assumed to have perfectly cylindrical twines. This leads to the following derivation of scaling factor Anet AFE = 1 4πd2 t d2 t = π 4 In total, this gives a geometric scaling factor Vnet VFE = Snπ 4 . In this work, the factor is used to scale the material density of the net, the added mass density and the buoyancy force fb acting on the net. 3.4.4 Material Modeling To capture the response of the net structure consisting of twines and meshes using the continuous surface, some modifications have to be made to the material used. The main reason for this is that, as seen in figure 3.11, the in-plane shear response of the net structure differs significantly from what would be expected of a surface where material presence is more uniform. Moreover, adjusting the out-of-plane shear response allows for control of bending stiffness if needed. 27 3. Computational Model (a) Undeformed (b) Shear Deformation Figure 3.11: In plane shear deformation of net structure This behavior is controlled by components D44, D55 and D66 in equation 2.8. For the circular net cage used by Moe-Føre et al. (2016), a cylindrical coordinate system as described by figure 3.12 is used to define the material response. Figure 3.12: Cylindrical coordinate system used to define material behavior of circular net cage As a basis, the net structure is represented as an isotropic material described by equation 2.6. However to appropriately capture the structural behavior displayed in figure 3.11, factors are added to each shear component of the constitutive matrix to allow for reduction of the stresses created by shear strains in the individual directions. This results in the following constitutive relation, in Voigt notation  σr σθ σz σrθ σrz σθz  =  2µ + λ λ λ 0 0 0 λ 2µ + λ λ 0 0 0 λ λ 2µ + λ 0 0 0 0 0 0 krθµ 0 0 0 0 0 0 krzµ 0 0 0 0 0 0 kθzµ  =  ϵr ϵθ ϵz 2ϵrθ 2ϵrz 2ϵθz  (3.13) In this material model, shear factors krθ, krz and kθz and Young’s modulus E in equation 2.7 are tuned to give appropriate results. Due to the lack of numerical information available in the experimental data, this process is carried out using the visual deformation presented in figure 4.1. This is considered to give a sufficiently accurate material model for the purpose of the present work, however for quality analyses of specific nets material data would have to be obtained from experiments. 28 3. Computational Model 3.4.4.1 Material Damping To improve the numerical stability of the model, Rayleigh damping is implemented in Abaqus. As the intention of the model is to model net systems of large scale, and the surface approximation neglects some of the physics that are important to capture small scale effects, it was considered preferable to preserve large scale rigid body modes. For this reason, α of equation 2.10 was set to zero and all damping was imposed through the stiffness proportional term of the damping matrix C. To suppress any high frequency eigenmodes that were not considered important to the validation results in the present work, β of equation 2.10 was set to 0.1. In any analysis containing oscillating loads from waves, vortex shedding or other phenomena, one would have to ensure that this damping does not interfere with any of the eigenmodes of interest triggered by the periodic loads. 3.4.4.2 Compressive stiffness When modeling the structural behavior of the net structure, it is important to properly model the lack of bending stiffness of the structure. On a large scale, this is achieved to a sufficient level through the use of shell elements with a very small thickness. However at the scale of single meshes, the bending of individual twines will not be captured unless the mesh size is small enough to resolve bending at the millimeter scale. In practice this mesh size would generate models that are too large for practical use within the FSI framework. The most natural way to circumvent this issue is to allow the material to compress freely in both in plane directions. This would model the length change of a net section with small scale bending without resolution of the bending itself. Preferably, this would be achieved by alteration of the material response to strains ϵθ and ϵz of equation 3.13, where negative strains should give little or no resulting stresses. In Abaqus, the closest available option for implementation of this behavior is to remove all stress response to negative strains, removing all compressive response. Despite the material damping described above, this causes model stability too deteriorate to a point where it is not possible to converge the model using time steps of reasonable magnitude. For this reason, small scale wrinkling of the net structure has been fully neglected in the present work. 29 4 Model Validation For validation of the simulation model developed in the present work, the tests of Moe- Føre et al. (2016) were used. In the reported tests, 6 different net configurations were tested in a circular net cage with a height of approximately 1.50 m and a diameter of 1.75 m. For stability, 16 steel cylinders with a submerged weight of 4.48 N are distributed along the lower circumference of the net cage. Fluid forces acting on these cylinders were calculated and deducted from the final force measurements. Forces acting on the load cell structure and steel ring holding the net cage were measured beforehand and deducted from the final force measurements. The model, undeformed apart from deformation from weights, is displayed in figure 4.1. Figure 4.1: Circular net cage in its non-deformed state. Image courtesy of Moe-Føre et al. (2016) [15]. The circular net cage is subjected to constant currents, ranging from 0.12 m/s to 0.93 m/s. For each velocity, loads in the vertical and horizontal direction were measured over 30 seconds and subsequently averaged over the entire time span. Simultaneously, the deflection of the net cage is visually captured through images. Figure 4.2: Deformation of circular net cage of solidity Sn = 0.19 at current velocities [0, 0.26, 0.5, 0.76, 0.93] m/s. Image courtesy of Moe-Føre et al. (2012) [15]. 30 4. Model Validation For validation of the computational model, deflections were validated visually using figure 4.2. Although this does not result in any exact comparison, as image perspectives from tests and simulations will not be identical, it provides a sense of how accurately the net is modeled from a structural standpoint. Alongside the visual validation of deflection, force measurements of the computational model have been validated using the aforementioned vertical and horizontal force data, presented in figure 4.3. Figure 4.3: Force data obtained from tests carried out by Moe-Føre et al. "N19" denotes a net with Sn = 0.19 and so on [15]. In the present work, all model validation was carried out towards the net with solidity Sn = 0.19, “N19” in figure 4.3. 4.1 Material Tuning For the purpose of validating the computational framework developed herein, it was con- sidered adequate to tune the material model presented in 3.4.4 using the limited available experimental data. The aim of this process was to obtain material parameters E and kθz such that the FE model would give deflection patterns that are representative, but not necessarily identical to those of the experiments by Moe-Føre et al. The remaining material parameters were set to the constant values reported in table 4.1. Table 4.1: Constant material parameters Parameter Constant Value ν 0.45 krθ 1 krz 1 31 4. Model Validation From the experiments, the available data used to tune the material was the visual data of the net deforming at u∞ = 0, as well as the fact that the points where the weights were attached were reported to displace approximately 0.05 meters at this current veloc- ity. Using this data the curvature of the bottom circumference in the θz-plane could be matched sufficiently without having the deflection of the load points deviate to far from the reported 0.05 m. Figure 4.4: Side by side deformation comparison of net used by Moe-Føre et al. and FE representation of net in computational model The tuning gave the final parameters E = 150 000 Pa and kθz = 0.15, resulting in the deformation shown in figure 4.1. As stated in section 3.4.4, this tuning is considered sufficient for the purpose of the present work, and as will be shown in chapter 5 it provides more realistic deflection patterns than a fully isotropic model. However it should once more be stressed that for quality analysis of specific nets, a thorougher validation of the chosen material model should be carried out. 4.2 Model Sensitivity Studies To understand the impact of different input parameters on model output, a number of minor sensitivity studies have been carried out. These include studies to understand the impact of certain parameter values, as well as studies to understand the added value of certain modeling decisions. As with the model validation process, comparison between different cases was made using measured forces and visual deformation data. 4.2.1 Value of Model Components To understand the value of implementing the separate components going in to the FSI model, implementation was carried out in a step-by-step manner. The order in which the model components were tested is presented in table 4.2. As the orthotropic material model was considered the most realistic from a physical standpoint, this was used as the standard choice. However to understand its impact, an isotropic model was subsequently tested. 32 4. Model Validation Table 4.2: Order of implementation of FSI model components Model nr. Force Model FSI coupling Material model 1 Løland (1991) One-way Orthotropic 2 Løland (1991) w/ Re-dependence One-way Orthotropic 3 Løland (1991) w/ Re-dependence Two-way Orthotropic 4 Løland (1991) w/ Re-dependence Two-way Isotropic After implementation of all models, the two-way coupled model using an orthotropic ma- terial and the Re-dependent force model (denoted model nr 3 in table 4.2) was considered the main model, used as a baseline for remaining sensitivity studies. 4.2.1.1 Angle Dependence As there are clear uncertainties regarding the appropriate values of the Fourier coefficients and CD(θ = π/2) of equations 3.1 and 3.2, an additional set of values apart form the standard a1 = b2 = 1, CD(θ = π/2) = 0.04 was tested. The coefficients were determined using the circular net cage experiments of Zhan et al. where a rigid net cage of varying solidity is subjected to constant current [27]. As the net cage does not deform, it allows for simple analytic determination of the forces using the force model described in section 3.1 without having to use the FSI model. In its turn, this allows for efficient optimization using SQP or similar optimization algorithms. Moreover, the circular cage provides a uniform distribution of inflow angles in the range θ ∈ [0, 90]°, making it a good geometry for optimization of angle dependence. In the optimization, the objective function used was the relative root mean square error (RRMSE) between experimental results and model results across all testing points. The number of coefficients available to the model was set in line with the results of Martin et al., however as the net cage is symmetric and rigid, no lift force data is available and the lift force coefficients are not changed. This results in the following negative null form min RRMSE(a1, a3, a5, CD(θ = π/2)) subject to a1 + a3 + a5 − 1 = 0, − CD(θ = π/2)) ≤ 0 (4.1) When calculating the analytic forces acting on the rigid cage, it was assumed that the flow velocity was reduced by 15% when passing through the front half of the cage. Hence, a current velocity u = 0.85u∞ was used on half of the net cage area. The op- timization of equation 4.1 was carried out using the SQP algorithm available in SciPy’s scipy.optimize.minimize() functionality [33]. This gave the optimizers presented in table 4.3. 33 4. Model Validation Table 4.3: Optimum values of Fourier Coefficients and parallel drag coefficient Coefficient Optimal Value a1 0.830 a3 0.148 a5 0.022 b2 1 b4 0 CD(θ = π/2) 0.0235 With the optimizers of table 4.3, an RRMSE of 10.3% was achieved between the exper- imental data and analytical model. As seen in figure 4.5, these coefficients provide good correlation at low solidities, however at larger solidities accuracy starts to drift off slightly. Figure 4.5: Comparison of experimental results of Zhan et al. and analytical results obtained from equations 3.1 and 3.2 using coefficients presented in table 4.3. As shown in figure 4.6, these coefficients should lead to a significant force reduction for nets where large sections have an inflow angle between 10° and 80°. As this is the case for the deforming net cage used for model validation, model nr. 3 in table 4.2 was run an additional time with the optimized coefficients to understand the impact of the angle dependence on the force measurements. Figure 4.6: Optimized angle dependence of force model 34 4. Model Validation 4.2.2 Parameter Studies In addition to the study of model component impact, a parameter study testing the sensitivity to uncertainties in input parameters was carried out. The general procedure for this study was to test one parameter at a time, testing one larger and one smaller value than that used in model nr. 3 in table 4.2. The tested parameters and parameter ranges are presented in table 4.4. It should be noted that the results of the parameter study are valid mainly for the circular net cage in question and similar geometries. As the force measurements used are shape dependent, mainly due to the angle dependence in equation 3.1 and 3.2, it is not certain that the same conclusions could be drawn for geometries that are deflected in different ways. Table 4.4: Tested value for each individual parameter in the performed parameter sensitivity study Parameter Default Value Tested Range FE size 0.05 m 0.025-0.1 m CFD cell size 0.0625 m 0.03125-0.125 m Solidity 0.19 ±15% Twine diameter 0.00242 m ±15% Young’s modulus 150 000 Pa 75 000-300 000 Pa To understand the level of mesh dependence present in the model, multiple values were tested for both FE target size and CFD target cell size. For FE target size, the element size seen in figure 3.8 was simply doubled/halved to understand if the different levels of resolution resolve different physical phenomena and eigenmodes. For CFD cell size, the target cell size in the refined region seen in figure 3.6 was doubled/halved to understand what impact this might have on the effect of the source term in equation 3.11. In reality, measuring the net solidity and twine diameter accurately is not entirely straight- forward. Generally, these parameters are measured on a non-submerged net using image analysis. However given that the net twines consist of nylon fibers creating a rather soft structure, it is unclear whether these parameters should be measured on a net in a submerged, non-submerged, wetted or dry state. Moreover, this matters not only when measuring the parameters for the net subject to analysis, but also when making the cor- responding measurements on the nets used when developing the net force model. Because of the parameter uncertainties, the model sensitivity to both measurements is examined by increasing/decreasing each parameter by 15%. In addition to mesh and net parameters, the model sensitivity to the Young’s modulus used in equation 3.13 is tested. This provides an understanding of to what degree the results of the present work are effected by the material tuning described in section 4.1. 35 5 Results In this chapter, the key results from the model validation and subsequent sensitivity anal- yses are presented. Along with the key results follow brief commentary highlighting the valuable takeaways from the presented data. For in depth discussion of the consequences of the results and conclusions about the developed model, be referred to chapter 6. Results containing all validation data is available in appendix A. 5.1 Model Validation From the performed model validation, a number of interesting results arise. This mainly includes some differences in physical behavior and importance of different components in the force modeling, CFD modeling and FE modeling. 5.1.1 Full Model In figure 5.1, the horizontal force measurements show the accuracy of the computational model at different velocities. Figure 5.1: Horizontal forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, nr. 3 in table 4.2. The first clear takeaway from the comparison is that at current velocities above 0.5 m/s, there is some physical behavior causing the force increase to decline that is not properly captured by the computational model. For velocities up to 0.5 m/s however, the model accurately predicts the correlation between current velocity and horizontal force, albeit while slightly over predicting the forces. In this range, the horizontal force is consistently about 15% larger than the corresponding experimental value. Above 0.5 m/s however, 36 5. Results the relative difference between experimental measurements and measurements from the computational model goes as high as 74%, which is the difference at 0.93 m/s. Figure 5.2: Vertical forces measured in experiments by Moe-Føre et al. compared to results from the full computational model, nr. 3 in table 4.2. When examining the vertical forces presented in figure 5.2, the same lack of decline in slope at high velocities can be identified. In comparison to the measured horizontal forces, it should however be noted that the vertical forces are not consistently over predicted by the computational model. Here, the relative error between experimental and computational measurements range from -28% to 22%, with an average error of -8.5%. 0.00 0.26 0.50 0.76 0.93 Current Velocity [m/s] Figure 5.3: Deformation of net cage in experiments by Moe-Føre et al. compared to the deformation of the net cage in the full computational model, nr. 3 in table 4.2. In terms of net cage deflection, figure 5.3 shows that the structural response of the compu- tational model can be considered identical for current speeds up to 0.26 m/s. At 0.5 m/s, deformation is still close to identical. However if one observes the part most upstream (left in figure), it can be noted that the observable curvature differs between the exper- imental image and the image of the computational model. As current velocity increases further, this discrepancy becomes clearer and clearer, and can also be seen on the most downstream part. 5.1.1.1 Verification of Source Term Effects As there are no proper data available regarding the current reduction behind the net and related effects, detailed validation of the source term effects is not possible. Nevertheless, it is possible to verify that the source term results in the expected effects on the velocity field. 37 5. Results (a) Side View (b) Top View Figure 5.4: Normalized velocity magnitude |u|/u∞ of flow field surrounding net cage at u∞ = 0.5 m/s. Instantaneous snapshot. In figure 5.4 it is shown that the source term added in equation 3.11 results in the desired velocity reduction behind the net cage. Additionally, it is shown that as a result of the current reduction downstream of the net, there is a small (∼5%) increase in velocity magnitude to the sides of the net cage ensuring conservation of mass. Figure 5.5: Distribution of angle θ across all finite elements at u∞ = 0.50 m/s After the first passage through the net cage, the current velocity is decreased to approx- imately 0.90u∞. After the second passage, the current velocity reaches approximately 0.83u∞. The velocity reduction obtained using the source term model can be compared to the reduced velocity model of equation 2.24. Assuming an angle θ = 45◦, which is not unreasonably large considering the distribution shown in figure 5.5, equation 2.24 gives CRF = 0.91. As there is no Reynolds number dependence in equation 2.24, we get ve- locities 0.91u∞ and 0.83u∞ after the first and second passage through the net cage. This shows that the source term model gives velocity reduction values that are in line with existing models, while including additional effects of mass conservation. 5.1.2 Impact of individual model components In f