A COUPLED ACTUATOR LINE AND FINITE ELEMENT ANALYSIS TOOL

. Fluid-dynamic loading of many solid bodies can be simulated using geometry resolving computational ﬂuid dynamics methods, where the body’s shape is resolved in the mesh. In some cases though slender bodies, like ropes or cables, spars, turbine blades, foils and lattice structures would require prohibitively high cell counts, since the geometrical features to be resolve are much smaller than the overall domain. Such bodies are usually made up of generic cross sections like round, square or standardised technical proﬁles like the famous NACA digit series for which good parametrisations of reaction forces to incoming ﬂow exist. Actuator line methods apply inﬂow dependent reaction forces to the ﬂuid domain, thus allowing the computationally eﬃcient simulation of slender bodies and have been used extensively, for example in wind and tidal turbine simulations. Beam elements representing slender bodies are standard building blocks in structural ﬁnite element models. Combining actuator line theory with a ﬁnite element beam model allows the eﬃcient simulation of ﬂexible structures under ﬂuid loads, like turbine blades or nettings used in ﬁsh farms. This paper presents an implementation of such a coupled model in OpenFOAM based on the existing turbineFoam actuator line model. The underlying numerical method is detailed and ﬁrst test cases are provided.


Introduction
With sufficient spatial and temporal resolution even the most complex flow effects like flow separation and highly turbulent flows can be resolved using computational fluid dynamics methods (CFD) [1].
However, for practical design purposes, computational resources are limited and simplifications are still needed. Computational domains in marine engineering typically extend over several wave lengths or the water depth, and thus often range over 100s of metres. Challenges arise especially if details of the structure of interest are orders of magnitude smaller than the domain, as is frequently the case for lattice structures and foils, for example for tidal turbines, mooring systems, fish farms and wave energy converters.
Correctly resolving a NACA0012 profile in all operating conditions at a Reynolds Number of 6 million requires roughly 0.5 million cells in only 2 dimensions [2]. For most industrially relevant three dimensional problems the same level of resolution results in a very high cell count.
Slender structures have been studied extensively and parametric descriptions or lookup tables exist for most features like lift or drag coefficients. Actuator Line Models (AL) use such parametric descriptions of profiles to apply the resulting force to the flow field, without representing the structure in the mesh [3,4,5]. AL simplify the meshing, since only the fluid domain needs to be created. The resulting lower mesh size can reduce the computational burden and thus allows to focus on resolving far field effects like wakes in wind and tidal turbines [6,7,8,9,10]. Applications to other marine structures are still limited, one interesting exception is the work on the underwater kite system by Fredriksson et al. [11].
An interesting feature of slender bodies is their structural response. Hydrofoils are sensitive to the angle of attack (AoA) and thus torsional deformation. Mooring lines, nets and lattice frameworks can deform significantly, especially in dense fluids. Wind turbines are sometimes designed to deform into their ideal shape under wind loading. Tidal as well as wind turbines have been built to shed excessive loads by bend-twist coupling, see Nicholls-Lee et al. [12]. This paper describes a first version of ALFEA, implementing the coupling of an AL toolbox with Finite Element Analysis (FEA) to enable the simulation of fluid structure interaction of arbitrary slender structures. The focus in this paper lies on marine engineering applications. However, the methods should be applicable to a much wider range of physical problems. Similar implementations like Sale et al. [13], Neumann [14] are limited to turbine blades, do not seem to be maintained any more or are not available as open source code.
The toolbox is based on OpenFOAM R and extends the AL method originally implemented by Bachant et al. [5] with a custom written FEA model to take structural deformations into account.
The paper is structured as follows. The following section 2 provides an overview of the FEA module. Section 3 details the coupling with the AL method. A further section 4 provides details on test cases. Conclusions are drawn and an outlook on future applications in section 5.

Frame Analysis Implementation
A wide range of open source FEA methods are available and could be used for coupling. However, most FEA methods offer many more features than could possibly be used in this application and are not written in C++. To ease maintenance of the code base and integration with OpenFOAM R , it was decided to write a stand alone C++ library with only the required functions.
This section provides a basic overview of the code design and implementation. The method used closely follows the formulation in Przemieniecki [15]. Beam elements are based on 2 nodes describing the end positions of each element and can evaluate bending along and torsion around their length axis and extension/compression due to normal forces according to linear theory. The current implementation employs a steady state formulation.
A first implementation in GNU octave [16] was used for initial development and validation. In a second step the code was reimplemented in C++ making use of the Armadillo library for matrix manipulation [17]. OpenFOAM R [18] offers its own syntax for matrix operations but since the Armadillo library was specifically written to provide a similar interface to GNU octave, its use simplified the reimplementation of the code signficantly. The entire frame analysis class definition consists of less than 800 lines of code.
All frame analysis code is located in a new folder called FEA, located inside the original turbinesFOAM library inside the src folder.
The FEA mesh is defined by lists of node locations and connections. The constructor of the novel FrameAnalysis class calls the beam1() function, which iterates over all elements and by calling the elemstiff() function creates the respective element matrices K e in local coordinates.
Here E is the Young modulus, A the cross section area, L the length, I z and I y the cross section inertias, G = E 2(1+P ) the shear modulus with P the Poisson coefficient. elemstiff() also creates the transformation matrices required to transform from local (element) to global coordinates. The function call assem() then assembles the overall system matrix K in global coordinates. K is then reordered by the function neworder() according to the boundary conditions provided, yielding separate submatrices for free ( f ) and restrained ( r ) nodes, with u being the deformations and F the forces: Algorithm 1 provides a pseudo code to help understand the code structure.
The armadillo library's solve() function is then used to obtain the deformations at the free nodes by solving and the reaction forces are then obtained by evaluating

Algorithm 1 FrameAnalysis Workflow
Import nodes, elements.. for all elements do elemstiff() create element stiffness and transformation matrix K e and RotM at assem() Assemble system stiffness matrix K end for neworder() find free and restrained nodes reorder() Partition K according to free and restrained nodes, such that The armadillo library employs its own matrix type definition which differs from the lists used in the AL model. Helper functions to convert lists of floats or integers into armadillo matrices were implemented and are called List2Mat() and List2intMat(). The function Mat2List() allows to reconvert data from armadillo matrices to lists.

Coupling with AL model
Details and first validation cases of the AL method are provided by Bachant et al. [5]. This section focusses on the coupling methodology.
As shown in the schematic in Fig 1, forces are evaluated at the midpoint of each AL element, while in the FEA model forces and restrains must be defined at the nodes or endpoints of each element. The FEA mesh is thus generated with twice as many elements as the AL definition. Forces evaluated by the AL for a given node position and orientation are passed to the FEA model. These fluid forces are then applied to every second node in the FEA mesh. The use of equal numbers of FEA and AL elements would also be possible, but require the evaluation of equivalent loads and thus a significant increase in the codebase and runtime. While it is unclear whether the additional operations would outweigh the benefit of a smaller FEA system matrix, it should be noted that in typical application scenarios the FEA matrix is expected to be orders of magnitude smaller than the matrices involved in the fluid solver. The FEA simulation are thus not believed to be a bottleneck for the overall simulation and efficiency considerations less relevant.
Resulting deformations from the FEA run are then applied to the AL method by updating the node positions and re-evaluating the AL elements' centre pitch and fluid force properties. Restraints can only be defined for the first or last element of a section and are applied at the first or last node respectively.
All required FEA inputs are defined through the AL model definition in system/fvOptions as shown in listing 1 which was extended to allow for the definition of material and geometrical properties like the Young modulus, section area/inertias and restraints. § The original AL model interpolates section values linearly along a geometric section for any wanted number of elements, this same interpolation is also applied to FEA properties except the restraints, which are only applied to the first or last node of each section.
The AL evaluates source terms for the momentum equation. The force acting on the flow is thus applied as a discrete source term every time the velocity equation is solved. For each AL element, the last available flow field is used to evaluate parameters like the AoA and velocity magnitude. Each actuator line element also saves the applied forces and moments. The differential between new and current time step is applied in the FEA model, such that a constant load in time will not lead to further deformation after finding the equilibrium state. The FEA model is solved at every new time step and the positions and rotations returned to the AL model. The AL then updates element positions and profile pitch is adjusted by the negative cross product of the elements' span direction and FEA rotation vector. A high level pseudo code of the algorithm is provided in Algorithm 2.

Test Cases
AL methods have been tested extensively for a wide range of applications [5,3,19,20,11], the test cases presented here thus focus on the FEA part.
The fluid domain is identical for all cases and based on the original static actuatorLine tutorial case from the turbinesFoam library. In all cases the pimpleFoam solver was used with OpenFOAM R -8 in transient mode with 2 outer correctors. The kEpsilon turbulence model was employed and cases were run for at least 2 seconds. Resulting deformations are always presented for the last timestep unless explicitly defined differently. A steady-state simulation might save some computational effort and would also have sufficed for demonstrating the steady state FEA implementation. However, to ease testing for users already familiar with the AL methods and because we plan to extend the method to transient cases, we chose to keep the fluid solver settings of the test cases and setup unchanged.
The domain is a cube with a bounding box with extreme positions (-1.52, -1.83, -1.22) and (2.16, 1.83, 1.22). The basic mesh is refined with 32, 32 and 24 cells in x, y and z direction respectively. snappyHexMesh is used to create a zone with refinementLevel 2 within a cylinder with 0.5 m radius between points (0 0 -1.5) to (0 0 1.5). An example of the resulting mesh is shown in Fig 2. Boundary conditions used in the fluid simulation are presented in Table 1, numerical values used are given in listing 2.
4.1. Test Case 1: Bending. A first test was devised to evaluate deformation due to bending. Case files can be found under turbinesFoam/tutorials/actuatorLine/bendingtestFEA. A blade profile of L = 2 m length is fixed in rotation and deflection in all directions at the position (0,0,-1) and points in z direction as shown in and the analytical solution for deflection δ under constant load is [21] δ(x) = qL 4 24EI z End effect models, used to correct for the loss in lift close to wing or blade tips were not used. Although unphysical this setting allows better comparison with a simple analytical load model of constant load along the beam.
and a deformation angle of Variable names and values used are given in Table 3  Table 4 lists resulting deformation angles from Eq 4.4 and the ALFEA simulation (AL). Eq 4.4 yields to 5.5% less deformation than the simulation. The cause for this deviation is again the difference between analytically derived force and ALFEA results. Results obtained using Eq 4.4 with the loads from the ALFEA model (Eq Sim ) yield excellent agreement with the ALFEA simulation. 4.3. Test Case 3: Changing angle of attack by a deforming substructure. A last test scenario demonstrates the change of AoA by the deformation of the supporting structure. Two cases are simulated. Case 1 simulates a deforming support structure with the inertia of the supporting beam set to 7.E−11 m 4 . Case 2 simulates a virtually rigid support structure with the inertia of the supporting beams set to 7.E−6 m 4 . The two simulations can be found under tutorials/actuatorLine/changeAoA and fixedAoA respectively. Key input parameters of the cases are listed in Table 5.    The upward force causes the deformation δ(x) along the length of the supporting beams and a deflection angle at the profile position [21] tan 2θ = F L 2 S 2EI z (4.7) The factor 2 in front of the angle θ and δ(x) accounts for the number of supporting beams. For a rigid structure the AoA remains constant and the C L only varies during the transient startup while the flow field develops. The soft support structure causes a reduction of AoA due to the deflection angle and thus limits the loading. Fig 7 shows the deflection of the beam under loading equivalent to the undeformed structure, with an AoA of 10 • (Eq) according to Eq 4.6. Deformation by a load corresponding to a reduction in C L by this maximum deformation is shown as Eq Iter and yields much lower deflections. ALFEA simulations (AL) yield to the mean deformation of those extreme assumptions. Eq Sim , obtained using the loading from the ALFEA model show excellent agreement with the simulation results.

Conclusions
This paper presents details on a first version of the ALFEA toolbox. Coupling of AL and FEA methods now enables the efficient simulation of flexible slender structures in OpenFOAM R . Details on the implementation are provided. Two test cases demonstrate correct application of the forces evaluated by the AL model to the newly implemented FEA toolbox for bending and torsional deformation. A third test case demonstrates interactions affecting the AoA and the change of fluid loading caused by the structure's deformation. Figure 7. Resulting deflection along the supporting beam according to simulation (AL) and three analytical solutions. Eq is the result assuming no deformation and using a C L corresponding to the AoA of the undeformed structure at 10 • . Eq Sim presents results with the load taken from the last simulation step and shows best agreement. Eq Iter assumes a load corresponding to the AoA after deformation by the load of the undeformed beam.