A new validated open-source numerical tool for the evaluation of the performance of wind-assisted ship propulsion systems

. Wind propulsion is envisioned as one of the solutions for the decarbonisation of maritime transport, as it offers high efficiency in terms of primary energy comsumption. Many wind propulsion systems already exist at various development stages, but the uncertainties over their performance is a strong obstacle to their adoption by ship owners. This paper presents a general method for the assessment of steady and unsteady performances of wind-propelled ships with 6 degrees of freedom, as implemented in the open-source program xWASP CN. Inspired by system-based modelling, the method consists in the independent modelling of the forces acting on the ship, as functions of the ship’s 6 degrees of freedom and environmental conditions. An original root-finding algorithm that leverages the specifics of the physical problem to find the steady equilibrium is presented. The method works either like a Power Prediction Program (PPP) or like a Velocity Prediction Program (VPP). As a PPP, the forward speed and course are fixed while the required propulsive power, leeway angle, heel, trim and sinkage are solved. As a VPP, only the course is fixed and the attained forward speed, leeway angle, heel, trim and sinkage are solved. This makes the method suitable for both hybrid propulsion and pure wind propulsion. The force models can be semi-empirical (usually requiring very little input data), based on preliminary experimental or numerical results (such as forward resistance curves, lift/drag coefficients and frequency-domain sea-keeping coefficients), or full-fledged flow solvers (e.g. potential theory, CFD). Thus the method is suitable for all design stages as each force can be modelled with several levels of accuracy depending on the input data available. Comparisons of an intermediate-level model with experiments on a 18-ft catamaran fitted with a Flettner rotor and a water turbine show good agreement for steady-state results.


Introduction
Maritime transport currently accounts for 2.89% of human carbon dioxyde emissions [1], for around 70% of the global freight [2]. As global maritime freight is expected to triple by 2050 [3] and as avoiding global warming reaching catastrophic levels requires achieving net zero emissions by 2050-2070, the International Maritime Organization voted a resolution in 2018 to cut in half the global carbon emissions of the sector [4] by 2050. This will lead the industry to change the means of propulsion of ships.
Several alternatives are being investigated to replace diesel engines on ships: combustion engines with alternatives fuels (e.g. methane or hydrogen), electric motors (with either batteries or hydrogen fuel cells) and wind propulsion. While it comes with a set of challenges (space * e-mail: moran.charlou@ec-nantes.fr on deck, port operations, crew training...), wind propulsion has by far the best efficiency because it does not suffer from conversion losses along a long production chain.
However, and as opposed to when wind propulsion was still used -before is was gradually abandoned for the benefit of steam propulsion -, today's transport industry is ruled by strict timing. Alternative fuels require little change to the current operation of ships, but wind propulsion necessitates precise knowledge of the ship's behavior and good estimation of the environmental conditions along the route in order to ensure on-time delivery.
While this knowledge can be acquired through experimental means, that is prohibitively expensive during design stages. On the other hand, physical modelling and numerical simulation have been used successfully in the past decades to deal with physical problems closely related to wind-assisted ship propulsion.
For example, the estimation of a ship's propulsion requirements is a standard procedure in modern ship design. It mainly relates to two aspects: ship forward resistance and propeller (and eventually engine) performance. Both have been extensively studied, and many models allow for a reasonably accurate prediction of a ship's propulsive performance at various design stages. Such models include empirical models suitable for early design stages, while potential flow solvers and eventually CFD (based on viscous flow solvers) are becoming an industry standard for later design stages. These methods may be called Power Prediction Programs (PPP), which aim to estimated the propulsion power needs and/or the propeller dimensions and/or the engine speed at a set forward speed.
However, ship power prediction is usually focused on a single degree of freedom. While this is fine for conventional propulsion where side forces are negligible compared to propulsion and forward resistance, it is not sufficient for wind propulsion where the propulsion systems usually produce side forces of the same order of magnitude than the propulsion forces. In case of wind propulsion, the side forces are balanced by leeway resistance, which is out of the scope of conventional PPPs.
It can be noted that ship design is also subject to strong safety regulations, which has led to the development of numerical and experimental methods for the assessment of stability, sea-keeping and maneuvering. Most notably, the latter includes models for side forces which may be considered for wind propulsion modelling.
In yachting, Velocity Prediction Programs (VPP) have been in use for several decades for racing (e.g. ORC's VPP [5]). Because very different boats may compete against each other, their performance are rated and a handicap attributed in order to limit the bias and make the races more dependent on the crew's performance than the boat's performance. Similar methods are also used to optimize racing yachts, albeit with more accurate models.
VPPs include more degrees of freedom (3 to 6) than PPPs, and their goal is very different: they aim to compute the velocity a given ship may reach in a set of environmental conditions. In this, they may seem closer to the needs of wind propulsion. They were developed with rather small sailing boats in mind, and thus the physical models they use are usually not suitable for very large cargo ships.
More advanced versions of VPPs include the estimation of dynamical performance in waves or gusty wind: they are referred to as Dynamical Velocity Prediction Programs [6].
Thus, one can see that conventional methods used in the naval industry and for yachting are ill-suited to deal with the specifics of wind propulsion for maritime transport. Therefore, dedicated methods have been developed, e.g. [7]. However -to the knowledge of the authors -, none of these methods is readily available to the industry (commercially or in open-source), open about the physical modelling and adaptable enough to be suitable for all design stages. Therefore, the present study proposes a modular method that can accommodate several control strategies and give steady and dynamic results in 6 degrees of freedom, with multiple levels of modelling depending on the design stage. The method is made available as the open-source program xWASP CN, developed at Ecole Centrale de Nantes. Results are obtained for a set range of wind speeds and wind angles, that can later be used in a weather-routing program. Steady-state results were validated against a prototype catamaran.

Model
The ship's position in the fixed frame R 0 is denoted X = [x, y, z] R0 , its velocityẊ = [ẋ,ẏ,ż] R0 projected in the ship's frame R b is denoted U = [u, v, w] R b , and its angular velocity (in the ship's frame) is denoted Ω = [p, q, r] R b . The angular position of the ship is defined using ordered intrinsic yaw, pitch and roll rotations [φ, θ, ψ] (respectively). In the following, the generalized position X u = [x, y, z, φ, θ, ψ] and velocity V u = [u, v, w, p, q, r] are used for convenience. Note that U =Ẋ and V u = dXu dt because of the composition of rotations. There is still a derivation relation between X u and V u such that there are only 6 independent states.

System-based modelling
The method used in this study is inspired by system-based modelling [8]. This approach is based on Newton-Euler's laws of motion for a rigid body (2 × 3 equations): where R 0 is an inertial (fixed) frame, a(G) is the acceleration of the center of mass G, m is the ship's mass, O is an arbitrary point and M O R0 is the dynamic moment (rate of variation of the angular momentum), which can be expressed as a function of the inertia matrix, the angular acceleration and the distance between O and G (omitted here for brevity). In the ship's reference frame, fictitious forces must be included, depending on the angular velocity and acceleration.
In system-based modelling, forces acting on a body or a set of bodies are modelled independently. Each force is computed as a function of the ship's position X u and velocity V u and environmental conditions. This approach has several major advantages: -Simplicity: the core is only the equation of motion -Modularity: each force can be modelled differently, without affecting the rest of the force models -Adaptability: suitable for any design stage depending on the level of accuracy of the force models -Extensibility: force models can be added and developed on top of existing ones It also has some drawbacks: -Decoupling of the forces: because the forces are modelled independently, coupling effects are ignored -Uncertainties: since the force models can have various and different levels of accuracy, uncertainties are harder to track The limitation of decoupled forces can be mitigated by merging two models into one when a strong physical coupling is expected (e.g. interaction between sails, or interaction between rudder and propeller).

Steady state
The steady state is the static equilibrium of forces when the ship is subject to steady environmental conditions. The equations (1) are simplified to: Additionally, assuming a rectilinear motion of the ship in the horizontal plane, some of the states in X u and V u are null or fixed: Since Ω = 0, all fictitious forces are null and equations (2) may be solved in the fixed frame or in the ship's frame equivalently. The forces equation and the moments equation may even be solved in different frames.
If there are less unknowns than equations, the corresponding number of equations (forces or moments static equilibrium) must be ignored.

Forces modelling
It should be noted that only a limited subset of the force models available in xWASP CN is presented and used in this study.

Forward speed resistance
The forward speed resistance is classically expressed as a drag coefficient (3) and broken down into a wave-making part and a friction part (4). The wave-making part scales with the Froude number F r, while the friction part scales with the Reynolds number Re. The wave-making part can be computed using a potential-flow theory (strip method or 3D Rankine panel method). The friction part can be estimated using the ITTC-1957 formula (5) and a form factor k [9].
Alternatively, the forward speed resistance can be estimated using semi-empirical models. Such models are established using regression over a limited set of hull shapes, and are thus only suitable for hull shapes close enough to the regression set. The most used model of this type is the Holtrop & Mennen model [10].

Wind modelling and wind forces
A logarithmic wind profile is classically used to describe the mean wind field in the surface atmospheric boundary layer. For neutral stability conditions, the mean horizontal velocity writes: Where u * is the friction velocity, κ is the Von Karman constant (κ ≈ 0.41), d is the zero-speed plane and z 0 is the roughness length. At sea, d = 0 and z 0 is between 0.0002m and 0.05m depending on the sea state. Since u * is not easy to estimate, it is more convenient to use this equation relative to the velocity u r at a reference height z r : Alternatively, a power law wind profile may be used instead of the log wind profile [7].
The main contribution of the relative air flow generated by the wind and the ship's velocity is quadratic with regards to the flow velocity V , and can be expressed using lift and drag coefficients: The lift and drag coefficients are usually functions of the incident angle of the relative air flow in passive systems, but control parameters can also influence them in active systems (e.g. Flettner rotors, suction wings). Such a model neglects the added mass caused by the disturbance in the air flow, which is not an issue for steady state assessment, but may be a mistake for dynamic cases.

Hull lift
One of the major challenges of wind-propelled ship simulation is the modelling of side forces acting on the hull. In naval design, hydrodynamic side forces are traditionally studied in maneuvering problems, in which they play an important role. Widely used models such as the Abkowitz [11] or MMG [12] models use Taylor series to express the forces acting on the hull in terms of products between hydrodynamic derivatives and ship horizontal velocity components (u, v and r). The hydrodynamic derivatives can be identified using a Planar Motion Mechanism or advanced computational methods based on viscous flow theory (CFD), or estimated using empirical formulae [13]. According to Tillig and Ringsberg [14], the empirical formulae are not well suited to the estimation of side forces on the hull in the context of ship propulsion. They propose instead to use low-aspect-ratio (AR = T /L P P ≤ 4) wing theory to express the side force with lift and drag coefficients (10,11), and Inoue's empirical formulation for the center of efforts [13].

Rudder and daggerboard
Daggerboards and rudders can be modelled using classical high-aspect-ratio wing theory (9). The 3D lift and drag coefficients (C L and C D ) can be estimated from 2D coefficients (c L and c D ) with the inclusion of the liftinduced drag (12). The 2D coefficients are easily found in tables for well-known profiles such as NACA, or can be approximated using empirical formulae.
The effective aspect ratio AR e is taken as twice the aspect ratio (AR e = 2AR) to account for double-body flow condition [14]. The local flow velocity V (and ensuing incident flow angle) is projected in a plane normal to the profile's span-wise direction.

Hydrostatics
Hydrostatic forces can be computed using a fully nonlinear model, integrating the hydrostatic pressure over a 3D surface mesh of the hull. Alternatively, and for small motion around the equilibrium position, hydrostatic forces can be computed using a 6 × 6 linear hydrostatic stiffness matrix K HS (13).
Where non-zero coefficients can be easily computed with (14) and all other coefficients are equal to zero for a symmetric hull.

Numerical method for steady ship behavior
The present method was implemented in the open-source program xWASP CN, licensed under the weak-copyleft Eclipse Public License v2.0 (EPLv2). It relies heavily on the ship simulator xdyn -developed by Sirehna (Naval Group) -for the computation of forces and the time integration (for dynamic computations). xWASP CN was developed in Python while xdyn, in which computationally expensive force models may be evaluated, was developed in C++.

Control strategy
In xWASP CN, two approaches are available:

A. Fixed speed (PPP) B. Fixed propulsion (VPP)
Strategy A corresponds to the Power Prediction Program approach: the ship's forward speed is fixed, the wind propulsion system produces some propulsive and side forces and the goal is to find the remaining propulsion required to maintain the set speed. As for a conventional PPP, this strategy allows to estimate either the required propulsion power or the engine speed if a propeller is modelled. This strategy allows the time of arrival to be ensured.
Strategy B corresponds to the Velocity Prediction Program approach: the propulsion is fixed, and the goal is to find the forward speed reached by the ship. In conventional VPPs, the fixed propulsion if null (only the sails propel the ship). But for hybrid wind-assisted ship propulsion, it may be set so that a minimal speed is guaranteed.

Root-finding algorithm
The solution of the steady-state problem (2) requires finding the root of an equation with several unknowns. Unfortunately, global multivariate optimization algorithms are very computationally expensive, multivariate root-finding algorithms are very sensitive to the starting point and local optimization algorithms are prone to falling into local minima. Overall, existing root-finding methods are not robust enough to solve the present problem with arbitrary force models. A simple but robust root-finding algorithm that leverages the physical couplings was developed to circumvent these shortcomings.
It works by approximating the coupled solution with decoupled root-finding procedures on each problem, and then using the approximate decoupled solution as the starting point for a coupled root-finding procedure. If the coupled root-finding procedure fails to find the coupled solution, the process is repeated until a coupled solution is found (or a maximum number of iterations is reached).
The method leverages the fact that fortunately, the system (2) is generally diagonal-dominant in the case of a ship, that is to say each force or moment equation is mainly influenced by one of the unknowns. In other words the equations are weakly coupled. With this assumption, force and moments equations can be associated with unknowns (Table 1).
These relationships serve a dual purpose. Firstly, it allows to choose which force equilibrium to ignore when there are less than 6 unknowns: without a rudder, the yaw moment equilibrium M z can be ignored; without a propeller model in PPP (fixed speed), the propulsion equilibrium can be ignored. Secondly, it is useful to build a robust root-finding algorithm. Couplings between the equations still exist at different degrees, which allows to separate the global problem into two sub-problems: the horizontal problem (equations F x , F y and M z ) and the hydrostatic problem (equations F z , M x and M y ). These two problems are still coupled, but the coupling between them is weaker than inside each of them.
In practice, an inner loop (Fig. 1) is used to solve the horizontal and hydrostatic problems (separately). In each of these two problems, the decoupled problems are scalar (1 equation and 1 unknown, e.g. F y (ẏ) for the leeway problem). It uses well-known scalar rootfinding algorithms (e.g. Newton, Secant or TOMS748) for the decoupled problems and well-known multivariate root-finding algorithms (e.g. Anderson, Broyden or Powell).
An outer loop (Fig. 2) solves the horizontal problem and the hydrostatic problem as decoupled problems, and the global problem as the coupled problem.

Description of the experimental setup
The Farwind energy ship prototype [15] was built from a Hobie Tiger catamaran and equipped with a Flettner rotor and a water turbine (Fig. 3). Among the two daggerboards and the two rudders originally installed, one of each was kept on the port hull for the tests. The boat was instrumented with 2 anemometers, an accelerometer, a GPS receiver, and a magnetic compass. The anemometers were placed on poles at the front of each hull at a height of z r = 1.26 m, thus avoiding perturbations from the rotor from upwind to beam reach. The main particulars of the prototype ship are summarised in Table 2. The rotor rotational speed ω r , the rudder angle δ r and the water turbine voltage U wt were also recorded during the tests. The ship speed (forward speed and leeway) were extracted from the GPS trace, the ship attitude (heel and trim) was recorded from the accelerometer readings and the True Wind Speed (T W S) and True Wind Angle (T W A) were deducted from the ship speed and the anemometer readings.

Flettner rotor model
The Flettner rotor installed on the Farwind prototype is made of two stacked segments for a total length of 2.67 m and a main diameter of 0.45 m, with bottom and top end plates. The two segments have the same length and are joined by a plate, effectively forming a 'middle plate'. The rotor is placed on the port hull, aligned with the front beam [15].
The rotor can be modelled using the wing theory formulae (8) with S = D r L r , but the lift and drag coefficients (C L and C D ) do not depend on the incident angle of the flow. Instead, they depend on the spin ratio SR = ωr * Dr 2V for a given aspect ratio λ = Lr Dr [16] (see Fig. 4). For the Farwind prototype, the rotor can be modelled either as a single rotor with λ = 6 -thus neglecting the middle plate -or as two separate rotors. In the latter case, the top rotor has λ = 3, but the bottom rotor has λ = ∞ because it is vertically bound and has no tip [16]. In practice, the two options give similar results in the range of SR that was tested.
An additionnal drag force with a coefficient of C D = 2.5 (regardless of wind angle) was applied with a reference area of S = 2.1 m 2 (same as the rotor) to account for the windage on the rest of the structure [15].

Water turbine model
The water turbine installed on the Farwind prototype is a Watt&Sea Cruising 600 with a 28 cm turbine diameter. It is mounted on a vertical swivel in order to always align itself with the relative water flow velocity. The force along the turbine axis can be modelled using a thrust coefficient C T (15) -which is actually rather a drag coefficient is the case of a turbine -, while the torque can  reasonably be neglected.
In general, C T depends on the tip speed ratio T SR = ωwtDwt

2V
. Here, the turbine rotational speed ω wt is controlled by the water turbine input voltage U wt . C T (U wt , V ) was established experimentally from towing tank testings [17], and interpolated for the simulations.

Other input data
The wave-making part of the resistance curve (5a) was obtained through the potential flow 3D rankine panel method REVA [18]. The friction coefficient C f was obtained using the ITTC-57 formula (5). A value of k = 0.3 was used for the hull form factor, which was validated from CFD computations by Farwind.
The daggerboard and rudder are assumed to have a NACA-0009 profile, with the corresponding lift and drag coefficients (5b).

Experimental data points
A set of static data points was obtained from tests with the Farwind prototype. Each static data point was obtained by averaging times series over a dynamic steady state period. Unfortunately, the tests did not produce results for a range of TWA at fixed TWS or for a range of TWS at fixed TWA, so the 15 data points have to be studied separately. These data points are presented in Table 3.
HDG and TWA denote the two control strategies used on the prototype. The HDG strategy is for heading control while TWA is for True Wind Angle control. Note that the standard variation of the True Wind Angle is not significantly lower in TWA control mode.

Computation results
The results of the simulations are shown next to experimental results on Figure 6. The ship speed is well predicted for all experimental data points. The leeway angle is also well predicted for port winds (data points 1-4 on Fig. 6b), but shows discrepancies for starboard winds (data points 12-15 on Fig. 6b). As for rudder angles, they are again better predicted for port winds.
The difference between data points sets 1-4 and 5-9 is of particular interest: they are both for port winds, but the former was obtained with heading control and the latter with TWA control. The simulation gives much better results for heading control. This could be interpreted as the influence of the PID controller's calibration.
The fact that the sideway motion is not as well predicted as the forward motion is expected: the hull lift forces were modelled using a semi-empirical method, which is in general less reliable than the other models used in the simulations.

Conclusion
The numerical method presented in this study was able to predict the steady behavior of a 18-ft catamaran equipped with a Flettner rotor and a water turbine, using separate models for the various forces acting on the ship, for beam reach conditions. As long as the domain of validity of the force models is respected, the global simulation model is expected to perform similarly for other wind angles.
This study also outlined the difficulty to predict side forces compared to longitudinal forces. This can be explained by the relative abundance of research dealing with conventional ship propulsion in constrast with the literature on wind propulsion, but also by the more complex physical phenomena at play with larger leeway angles (most notably vortex-shedding along the keel).
Future work shall include the validation of xWASP CN for unsteady cases, such as maneuverability tests and sea-keeping over regular and irregular waves.