Multidisciplinary Design Analysis and Optimisation of Floating Offshore Wind Turbine Support Structures - Coupled Model Solution Strategies

Floating Offshore Wind Turbines (FOWT) can be installed at the sites of the most abundant wind resource. However, the design uncertainties and risks must be reduced to make them economically competitive. The design and optimisation methodologies for FOWT support structures adopted up to date tend to follow a sequential analysis strategy. Since the FOWT system involves multiple distinct, highly coupled disciplines, its analysis and design are challenging. This paper presents an efficient implementation of a coupled model of dynamics in an optimisation process by applying a Multidisciplinary Design Analysis and Optimisation (MDAO) methodology. The coupling effects studied include the interdependence of the mean offset of the platform and the aerodynamic and mooring loads, as well as the velocity of the platform and the viscous damping. The trade-off between the solution accuracy and efficiency for the coupled and uncoupled models was quantified, and a range of iterative solvers were compared. The study showed that the coupling between the platform offset and the mooring and thrust loads has a significant influence on the values of the responses, converging at higher surge and pitch offsets, higher mooring loads, and at lower thrust. These non-conservative results demonstrated the criticality of the two-way coupling between the platform excursion and the mooring loads. Notably, the coupled solution was achieved at a relatively low increase in the total solution time (+16%), due to the high efficiency of Broyden’s method.


Introduction
Floating Offshore Wind Turbines (FOWT) have the potential to capture untapped wind resources in deep-water conditions.However, reduction of design risk through an improvement of the analysis and design processes is still necessary for FOWT to become economically competitive [1,2].The design and optimisation methodologies for FOWT support structures still largely rely on the oil and gas industry legacy, with a sequential analysis strategy.However, the dynamics of the FOWT system involve a larger variety of distinct and highly coupled disciplines, including the interaction between the vast aerodynamic and hydrodynamic loads.Although many of the optimisation frameworks successfully reflected the most critical coupling effects [3], some of the important interactions are often either not considered or modelled as one-way relationships due to prohibitively long simulation times.
A fully coupled model can be implemented in an optimisation process by adopting a Multidisciplinary Design Analysis and Optimisation (MDAO) approach, whereby the individual disciplines are solved by highly specialised routines, with the global system solved by variables connection and appropriate direct or iterative solution strategy [4].Coupled multidisciplinary analysis (MDA) can be seen as a single analysis that evaluates the feasibility and fitness of candidate designs [5].Partitioning a large coupled model into a series of specialised components has multiple advantages, including the use of purpose-selected solvers, exploitation of the system structure, and modular architecture that facilitates component replacements and system extension.
We aim to demonstrate that: i) MDA is crucial to quantify the FOWT design objective and constraints; ii) the choice of iterative solvers influences the efficiency of the coupled analysis.The research is a part of an R&D project towards the development of the MDAO framework FEDORA, currently conducted at the University of Strathclyde and Delft University of Technology.It aims to develop an unbiased, efficient and flexible FOWT support structure optimisation framework.

State of the art
The physical environment and dynamics of FOWT are complex, with strong interdependencies between the wind-and wave-driven responses [6].The aerodynamic loads significantly increase the mean and maximum surge and pitch responses [7,8].Likewise, the amplitude of the heave response increases due to the vertical aerodynamic load component [8].This leads to a significant change in the operational position of the mooring system, which in turn leads to a shift in the 6 Degrees of Freedom (DOF) natural frequencies [9,10,11].Conversely, the wave-induced rigid motion of the platform substantially affects the tower-top motion, hence influencing the magnitude and variance of the rotor thrust and torque [7,12,8].Aerodynamic loads also amplify the coupling among the 6 rigid DOF, with greater pitch-heave coupling observed due to a positive mean pitch angle through altered hydrodynamics and changes in mooring lines tension [11,13].Indeed, large displacements of the platform can induce a non-linear behaviour (force-displacement relationship) of the mooring lines [14,7] and alter the potential and viscous hydrodynamic coefficients and forces on the floating substructure [15,11].
The coupled multidisciplinary frameworks developed up to date were extensively reviewed in [3].Due to the high computational cost, fully-coupled time-domain models of dynamics are not usually applied in highly iterative studies such as the parameter sweep or design optimisation.A frequencydomain coupled model of dynamics is therefore desirable for these applications.However, as concluded in [11] and demonstrated in the above review, its development is challenging if the coupling effects mentioned above are to be included.
Very recently, a few comprehensive multidisciplinary optimisation frameworks have been developed.For instance, the frameworks developed by Pollini et al. [16] and Ferri et al. [17] employed coupled aero-hydro-elastic models of dynamics by assembling the relative effects into the linear frequency domain equation of motion.In [16], all nonlinear terms were linearised with respect to the undisturbed floater position, hence failing to account for the coupling between the loads and the mean position of the platform.Ferri et al. [17] improved upon that model by linearising the loads around the equilibrium position of the platform using the linearisation functionality of the time-domain code FAST [18] at each step of the optimisation process.The computational effort was reduced by using lookup tables for the hydrodynamic coefficients based on the dimensions of the platform.However, the ability to recalculate the coefficients at a new equilibrium position has been lost.
Ideally, the coupled model should represent all important coupling effects with two-way (or feedback) connections, where appropriate, forming the so-called cyclic system.If such a model can be conceived, an efficient iterative solution strategy will also be necessary.Both aspects are the objectives of this study.

Aero-hydrodynamic model
This section describes the aero-hydrodynamic model developed within the current study, starting with the global coupled model and the flow of information between the subdisciplines (or components; Section 3.1), followed by the particular disciplinary models (Section 3.2), and the cyclic subsystems (Section 3.3).

Coupled model
The global coupled model was constructed with disciplinary submodels connected through input-output relationships, as presented in the N-squared diagram [19] in Figure 1, leveraging the functionalities of the open-source OpenMDAO package [20].In the N-squared diagram presented here, the green boxes represent the disciplinary components (analysis tools), while the parallelograms represent the design variables (white) and the coupling variables (grey) passed among the components.The inputs are passed vertically and the outputs horizontally.The variables above the diagonal feed forward, while those below feed backward.For instance, the damping module receives information from the potential problem module (potential damping Bh), from the aerodynamics module (aerodynamic damping Ba), and from the viscous flow module (viscous damping Bv), and combines them into the total damping B, after which it passes the output forward to the RAO module.The information "flows" through the statistics and viscous flow components, to be finally passed back to the damping module through the viscous damping coefficient Bv.Feedback connections were introduced to model the two-way coupling between the chosen components, forming two cycles, as also illustrated in Figure 2.

Fundamental disciplinary model
The floater and tower are assumed rigid structures with 6 DOF.A linear frequency domain approach is adopted, where the response of the platform to wave excitation is computed based on the solution of the 6 DOF forced oscillation problem.The complex amplitude of motion in the j − th direction due to excitation in the i − th direction is obtained by the following equation [21]: where a is the wave amplitude, M ij and C ij are the mass/inertia and stiffness matrices, respectively, A ij , B ij frequency-dependent added mass and radiation damping matrices, respectively, and X i is the frequency-and direction-dependent complex wave excitation matrix.
Inertial and restoring models The inertia matrix (M ij ) is computed through standard analytical formulations (see [22]).Gravitational stiffness is calculated using standard analytical formulations (see [21]), hydrostatic stiffness is calculated by an external module Capytaine [23], while the mooring stiffness is computed by MoorPy [24].This quasi-static mooring model is used twice: i) to obtain the total load on the platform from the three mooring lines -a part of cycle 1 in Figure 1, ii) to obtain the tension in the mooring lines and the mooring stiffness matrix in the equilibrium position.

Potential coefficients
The potential module includes an interface to the Boundary Element Method 3D panel code NEMOH [25], developed at the École Centrale de Nantes (more specifically, a Python version of the code, Capytaine [23]), which computes the first-order potential wave loads (added mass, radiation damping, and diffraction loads) for a floater in a given position and orientation.

Aerodynamic damping
The aerodynamic thrust load was precomputed for a range of wind speeds and tower inclination angles, varying the wind speed in small increments to numerically compute the derivatives of thrust with respect to wind speed, i.e., damping [26].The coefficients obtained this way were stored in a lookup table for the use of the framework.The state-of-theart aerodynamic time domain code employed to calculate the thrust derivatives, AeroDyn [27], uses Blade-Element Momentum Theory with the Pitt and Peters skewed wake correction model, as well as the Prandtl tip-loss and hub-loss corrections.The effect of varying blade pitch angle and generator speed as control actions has been taken into account, by setting the two values at each operational point according to the control strategy.

Short-term statistics
The evaluation of design constraints requires knowledge of the system motion response statistics.By definition, the standard deviation of a response α is obtained based on the zeroth moment of the respective response spectrum, which is obtained based on the product of the sea spectrum and the square of the module of the respective RAO, as in Equation 2: (2)

Extended dynamic model
This section presents the models implemented to introduce the most critical coupling effects into the global model.

Viscous damping and response statistics
The potential flow-based hydrodynamic approach presented above cannot account for the viscous hydrodynamic loads.The current work reintroduces the viscous load by adding a viscous damping coefficient to the potential and aerodynamic damping coefficients.This is achieved through a stochastic linearisation of the drag term of the Morison's equation [28].For each horizontal slice of the floater (dL) the viscous drag (dF d ) is given by: dependent on the drag coefficient (C D ), characteristic dimension (D), relative normal velocity u and its standard deviation (σ u ).Viscous damping is computed according to its definition: Since the drag load given by this approximate expression depends on the platform response statistics (standard deviation of velocity due to platform motion -σ u , in Equation 3), and the calculation of the response requires viscous drag as input, the viscous damping and the motion response form a cycle (cycle 2 in Figure 1).

Aerodynamic and mooring loads, and mean offset
The static equilibrium of the system in 6 DOF is defined through Equation 5: Aerodynamic thrust load (F t) is obtained from a pre-computed lookup table (see Section 3.2).The total load on the platform from the mooring system (F m) is computed through a quasi-static mooring model, as described earlier in Section 3.2.Since both the aerodynamic and mooring loads are to be computed at the mean offset of the platform (η), and the calculation of the mean offset relies on the values of these loads, an iterative solution approach must be adopted, as indicated in Figure 1 through cycle 1.
The Equilibrium, Aero, and Mooring components form a group.The implicit Equilibrium component is solved by driving the residual of Equation 5to zero, by varying the state vector (η), while keeping the states of components Aero (F t) and Mooring (F m) fixed.At the group level, all three states (η, F t, F m) are found iteratively, following one of the strategies defined in Section 4.
As demonstrated, this approach allows updating the values of the loads based on the rotor plane inclination and platform offset in surge and heave.Hence, the two-way coupling of the aerodynamic and mooring loads with the attitude and position of the floating system is approximately captured.

Model verification
The dynamic model has been verified by comparison with the results obtained from the time domain (TD) solver OpenFAST [18], set up to closely resemble the capabilities of the current model.The simulated scenario included the case of the OC3 Hywind spar platform [29] with the NREL 5MW reference wind turbine [30], subjected to 11.4 m/s wind and a moderate sea state (significant wave height of 3.66m and peak period of 9.7s).
The results of OpenFAST and FEDORA (the current approach) are reported and compared in Table 1.A very close agreement has been observed between the mean and standard deviation of the platform motions, line tension, and nacelle acceleration obtained with the two models.As presented in Figure 3, the frequencies and magnitudes of the responses are all sufficiently well captured.Therefore, the model is considered suitable for design optimisation purposes, within the conditions assumed.A generalisation of the verification results to a wider range of floating platforms and environmental conditions requires further study.Note that the current frequencydomain model's execution time is in order of seconds, while the higher-fidelity time-domain model used as a benchmark executes in a matter of minutes or hours, depending on the number of random seeds for stochastic sea states generation, as well as the time series length and sampling rate.

System Solution Strategies
In this work, five different iterative solvers are considered, differing in how the coupling variables (i.e., the responses passed from one component to another) are handled.Nonlinear Block Jacobi solver (NLBJ) [31] (i.e., fixed-point iteration), starts with an initial guess of all coupling variables (û (0) ), and then runs all components simultaneously to update the complete set of coupling variables (û i ), all at once.This is repeated until convergence is achieved.Since the input to the components is provided from the previous iteration, all components can be run in parallel without communication.Nonlinear Block Gauss-Seidel solver (NLBGS) [32] solves each component in turn, using the latest values of the coupling variables available, hence usually converging in fewer iterations than the Jacobi solver.However, since the components must be executed in sequence, it cannot benefit from parallel computing, and hence each iteration may take longer.A comparison of the two solution strategies for a two-components coupled system is shown in Figure 4. Nonlinear Block Gauss-Seidel solver with Aitken relaxation (NLBGS+A) [33] is a solver where the coupling variables are updated not simply based on the latest values available, but on the latest values available modified by the relaxation factor θ, as in: where û(k) i is the next estimate of the variable, ûtemp is the value of the component's response based on the values of all other components (as would normally be computed in the standard Gauss-Seidel solver), and ∆û (k) i is the step (difference between the current and the previous update for the component).The relaxation factor is adapted on each iteration.Through dampening of the divergence, quicker convergence is expected [5].
Monolithic Newton's solver can only be applied to the residual form of a coupled system, i.e., the states and residuals of the component must be available.It takes the residual of the component's governing equation and iteratively searches for the states that satisfy the residual equation.At each iteration, the solution is estimated by approximating the residual as a linear function and finding the step that makes it zero (this requires a linear solver).
Broyden solver [34] modifies Newton's method by using a finite-difference approximation of the Jacobian (or derivatives, in the 1-dimensional case).It also uses an approximated inverse of the Jacobian instead of solving the linear system, therefore, in some cases, it can be more efficient.This solver often takes more iterations to converge, but each of these iterations is cheaper (i.e., it requires fewer operations).
The convergence of cycle 1 does not depend on the solution of cycle 2 (and vice versa).Hence, the best strategy for each cycle can be studied independently.For cycle 1, two tests were performed, each with one of the two solvers suitable for the implicit system: Newton and Broyden.For cycle 2, five nonlinear solvers: NLBJ, NLBGS, and NLBGS+A, Newton and Broyden were applied, tested, and compared.Note that the NLBJ and NLBGS solvers cannot be applied to cycle 1, as they require the equation of equilibrium (Equation 5) to be solved for η, which cannot be achieved since the hydrostatic stiffness matrix is singular and does not have an inverse (this is due to the lack of stiffness in the surge direction, which makes the first row of the matrix a zero vector).

Results and Discussion
The case considered here is a spar-supported FOWT with a catenary mooring system [29] with the NREL 5-MW reference wind turbine [30].The environmental conditions are as described in Section 3.4.Subsection 5.1 compares the uncoupled and coupled models, while Subsection 5.2 reports the results of the solvers benchmarking exercise.

Influence of the introduction of the feedback loops
The coupled model, as introduced in Section 3.1 and sketched in Figures 1 and 2, includes two cyclic subsystems: cycle 1, coupling the mean position of the platform in 6 DOF (η avg ) with the mooring and aerodynamic loads (F m, F t), and cycle 2, coupling the viscous damping with the significant deviation of the platform motion velocity (B v , σ u ).The comparative study presented here is related to the model with both cycle 1 and cycle 2 present (coupled model), and the models with each cycle removed in turn, one at a time.The influence of eliminating coupling effects on the efficiency of the solution was studied by measuring the average simulation time over 100 runs.2, the introduction of the coupling between the platform offset, thrust and mooring loads increased the solution time by 16%.Although this is a relatively insignificant increase in the case of a single simulation, it becomes important when performing an optimisation involving a vast number of iterations or when exploring a large design space for parametric analysis.Moreover, the computations performed within the cycles were relatively inexpensive, the most time-expensive block of computations (BEM solver) performed before entering the first cycle.The assessment of the time cost of introducing the feedback loops highly depends on the fidelity (and cost) of the models involved within the cycles.The coupling effects had a noticeable impact on the mean surge offset and pitch inclination angle: the uncoupled model underestimated these responses by 31% and 11%, respectively.This was mainly due to the inability to account for the nonlinearity of the mooring loads with respect to the platform displacement (in the uncoupled model, the offset is computed by equating the thrust load with the product of the linear stiffness and displacement, evaluated at the undisplaced position of the platform).These non-conservative errors should generally be avoided, as excessive displacements of the platform may induce augmented mooring line and power export cable loads.In fact, the uncoupled model underestimated the mean mooring tension by approximately 30%.The uncoupled model also underestimated the heave offset by more than 50%, as it did not account for the vertical component of the thrust on an inclined rotor, and due to the lower estimated vertical force from the mooring system.The result can be substantially improved by performing just one iteration of the Newton solver, reducing the errors to below 10% for all offsets and forces considered (except for the heave, for which the improvement was small).
The significant amplitude of the nacelle acceleration was almost unaffected by the representation of the coupling effects.For a given FOWT design, it only depends on the integral over a function of the RAO in surge and pitch, and the sea spectrum S ζ (which is fixed).These, in turn, were found to vary only slightly, except in the narrow region near the peak frequencies (Figure 5).The differences were due to the differences in the aerodynamic damping (higher damping in the coupled system) and in the potential coefficients (slightly lower wave excitation and higher hydrodynamic added mass in the coupled system).
Overall, the judgment of the trade-off between the gains in accuracy and the increase in computational effort is highly dependent on the fidelity of the models used within and outwith the cycle.In the case currently considered, the errors due to the uncoupled representation of the system were unacceptably high, but a vast improvement in the accuracy of the solution has been achieved with a relatively small increase in computational effort.

The effect of cycle 2
The differences between the coupled and uncoupled solutions were marginal (below 0.5%) for the structure and the load case studied.This was due to the relatively low contribution of viscous damping to total damping (Figure 6).The introduction of this cycle increased the solution time by approximately 1%.Considering the low additional cost, the cycle should be kept in the system.Although its influence on the results presented in this paper was negligible, this may not be the case when different environmental conditions (e.g., higher sea states) and floating structures (particularly those with small-volume members) are considered.

Influence of the choice of the iterative solver
The performance of the two iterative solvers for cycle 1 and the five iterative solvers for cycle 2, as defined in Section 4, was compared on the basis of the average execution time over 100 runs.When studying cycle 1, the solver for cycle 2 was set to NLBGS+A, and when studying cycle 2, the solver for cycle 1 was set to Broyden's method.
As demonstrated in Table 3, the choice of iterative solver influenced the solution time to a degree.Between the Newton and Broyden solvers applied to cycle 1, the latter was more timeefficient (−4%).For cycle 2, the NLBGS solver with Aitken relaxation correction was the most efficient, converging over twice as fast as the same solver without the correction, and 60% faster than the NLBJ solver.The superior performance of the Broyden and NLBGS+A solvers can also be seen in Figure 7, which presents the convergence history of the two cycles.All solvers but the NLBJ achieved convergence in 10 or fewer iterations.For cycle 1, only minor differences were observed between the Newton and Broyden solvers, the former converging three iterations earlier.This agrees with what is usually observed (the Broyden solver tends to take more iterations, but each of these iterations tends to be less time-expensive, as explained in Section 4).For cycle 2, the NLBGS+A and Newton solvers only required four iterations to converge.While the NLBGS solver without the Aitken relaxation correction required just one more iteration, the NLBJ solver took considerably more, 25 steps to converge (note that only the first 10 iterations are presented in the figure).Parallel computing was not sufficient to offset the drawback of a higher number of iterations required.Again, the Broyden solver took two iterations more but converged the system considerably quicker than the Newton solver.
A single iteration of all solvers but NLBJ vastly improved the accuracy of the solution at an almost negligible increase in the solution time.

Conclusions
The article presented a coupled, multidisciplinary model of dynamics developed to efficiently quantify the design constraints within FOWT optimisation studies.
The framework of interconnected components (disciplinary analyses) was demonstrated, following the principles of the Multidisciplinary Design Analysis and Optimisation approach.The most important coupling effects were modelled, including the interdependence of the mean offset of the platform and the mooring and thrust loads, as well as the velocity of the movement of the platform and the viscous damping, which formed cyclic subsystems.The trade-off between the solution accuracy and efficiency for the coupled and uncoupled models was analysed, and a range of iterative solvers were benchmarked.
The coupling between the mean offset of the platform and the aerodynamic and mooring loads had a non-negligible influence on the values of the responses, converging at 31% higher surge and 11% higher pitch offsets, and at lower thrust and mooring stiffness.The uncoupled model underestimated the mooring tension by approximately 30%.These non-conservative results demonstrate the criticality of the two-way coupling between the platform excursion and the mooring loads.The dynamic characteristics of the system's global motion response, such as the nacelle horizontal acceleration, were relatively unaffected by the introduction of the feedback loops.The numerical implementation of these coupling effects was achieved at a low increase in the total solution time (+16%), due to the application of the efficient Broyden's method.
Regarding the coupling of the viscous damping with the statistics of the motion velocity, the discrepancies were marginal, mainly due to the low contribution of the viscous damping to the total damping and its minor effect on the response.Although all solvers were very efficient, the bestperforming Gauss-Seidel solver with Aitken relaxation massively outperformed the other solvers tested.
Since the results of benchmarking of the different solvers are highly problem-specific, but relatively inexpensive, it is recommended to always perform such an investigation before attempting a more comprehensive study such as design space screening or design optimisation, especially if time-intensive models are included within the cyclic subsystems.

Figure 3 :
Figure 3: Response Amplitude Operators of surge, heave and pitch motions, comparison between the results from OpenFAST and FEDORA.

Figure 5 :
Figure 5: Comparison of the magnitudes of the surge and pitch RAO for the coupled and uncoupled models of dynamics.C1 -cycle 1, C2 -cycle 2.

Figure 6 :
Figure 6: Surge RAO amplitude for different levels (sources) of damping: Bh -radiation damping, Ba -aerodynamic damping, B add -additional linear damping to calibrate with tank test, Bv -viscous damping.

Figure 7 :
Figure 7: Coupling variables convergence history for different solvers.

Table 1 :
OC3-Hywind Spar case study -mean and standard deviation of responses.

Table 2 :
Performance of the uncoupled and coupled models -the effect of cycle 1.
5.1.1.The effect of cycle 1 As reported in Table