© M. Charlou et al., Published by EDP Sciences, 2023

1 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 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.

2 Model

The ship’s position in the fixed frame R₀ is denoted X = [x, y, z]R₀, its velocity 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 and 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.

2.1 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):

(1)

where R₀ 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 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).

2.2 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:

(2)

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:

• No rotational velocity:

• No vertical velocity:

• Arbitrary horizontal position: x = y = 0

• Fixed course (to enforce TWA): ψ =0

• With PPP strategy: (forward speed) is fixed

Since , 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.

This leaves 4 or 5 unknowns: (leeway speed), z (sink-age), φ (heel) and θ (trim), and with the VPP strategy. Some other unknowns may be added from control parameters (e.g. engine speed w_m or rudder angle δ_r) up to 6 unknowns.

If there are less unknowns than equations, the corresponding number of equations (forces or moments static equilibrium) must be ignored.

2.3 Forces modelling

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

2.3.1 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 Fr, 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].

(3)

(4)

(5)

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].

2.3.2 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:

(6)

Where u_* is the friction velocity, κ is the Von Karman constant (κ ≈ 0.41), d is the zero-speed plane and z₀ is the roughness length. At sea, d =0 and z₀ 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

(7)

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:

(8)

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.

2.3.3 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_PP ≤ 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].

(9)

(10)

(11)

2.3.4 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 lift-induced drag (12). The 2D coefficients are easily found in tables for well-known profiles such as NACA, or can be approximated using empirical formulae.

(12)

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.

2.3.5 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).

(13)

Where non-zero coefficients can be easily computed with (14) and all other coefficients are equal to zero for a symmetric hull.

(14)

3 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++.

3.1 Control strategy

In xWASPJCN, 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.

3.2 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. for the leeway problem). It uses well-known scalar root-finding 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.

Fig. 1 Root-finding algorithm inner loop flowchart: η is the vector of unknowns, ηi and Fi are associated according to Table 1.

4 Validation

4.1 The Farwind energy ship prototype

4.1.1 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 dag-gerboards 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 w_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 (TWS) and True Wind Angle (TWA) were deducted from the ship speed and the anemometer readings.

Fig. 2 Root-finding algorithm outer loop flowchart: and ηhor make the horizontal problem, and ηHS make the hydrostatic problem.

Fig. 3 The Farwind prototype with its instrumentation. The rudder, daggerboard and water turbine are not visible through the water surface.

4.1.2 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_rL_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 for a given aspect ratio [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² (same as the rotor) to account for the windage on the rest of the structure [15].

4.1.3 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.

(15)

In general, C_T depends on the tip speed ratio . Here, the turbine rotational speed w_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.

Fig. 4 Lift and drag coefficients of a Flettner rotor with several aspect ratios, from [16].

Fig. 5 Input data for the force models, courtesy of Farwind Energy.

4.1.4 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).

4.2 Results

4.2.1 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.

4.2.2 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.

Fig. 6 Results of the simulation compared to experimental data.

5 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.

Nomenclature

R₀: Inertial (fixed) frame

R_b: Local ship frame

X: Position of the ship in the fixed frame,

: Velocity of the ship relative to the fixed frame, expressed in the fixed frame,

U: Velocity of the ship relative to the fixed frame, expressed in the ship’s frame,

Ω: Rotational velocity of the ship relative to the fixed frame, expressed in the ship’s frame,

φ: Heel angle (3^rd intrinsic rotation)

θ: Trim angle (2^nd intrinsic rotation)

ψ: Yaw angle (1^st intrinsic rotation)

X_u: Generalize position vector, X_u = [x, y, z, φ, θ, ψ]

V_u: Generalize velocity vector, V_u = [u, v, w, p, q, r]

a(G): Acceleration of point G

G: Ship’s center of mass

ρ_α: Air density

ρ_w: Water ensity

∇: Ship’s volume displacement

A_wp: Ship’s waterplane area

L_PP: Length between perpendiculars

S_w: Hull wetted surface

T: Ship’s draft

δ_r: Rudder angle

w_m: Engine rotation speed

AR: Aspect ratio (span/cord)

AR_e: Effective aspect ratio

GM_L: Longitudinal metacentric height

GM_T: Transverse metacentric height

k: Hull form factor

TWA: True Wind Angle

TWS: True Wind Speed

V: Local flow velocity

xB: Longitudinal position of the center of buoyancy

w_r: Flettner rotor rotation speed

A_T: Turbine rotor swept area

C_T: Thrust coefficient

U_wt: Water turbine command tension

Fr: Froude number (length-based)

Re: Reynolds number (length-based)

This work was supported by the French Environment and Energy Management Agency (ADEME) and the Region Pays de la Loire. The authors thank all the partners who have contributed ideas for the development of xWASP-CN, and Farwind Energy for providing experimental results and input data for their prototype. Finally, the authors would like to acknowledge the support of the Association Francaise de Mecanique (AFM) for the publication of the present article.

References

All Tables

All Figures

	Fig. 1 Root-finding algorithm inner loop flowchart: η is the vector of unknowns, ηi and Fi are associated according to Table 1.
In the text

	Fig. 2 Root-finding algorithm outer loop flowchart: and ηhor make the horizontal problem, and ηHS make the hydrostatic problem.
In the text

	Fig. 3 The Farwind prototype with its instrumentation. The rudder, daggerboard and water turbine are not visible through the water surface.
In the text

	Fig. 4 Lift and drag coefficients of a Flettner rotor with several aspect ratios, from [16].
In the text

	Fig. 5 Input data for the force models, courtesy of Farwind Energy.
In the text

	Fig. 6 Results of the simulation compared to experimental data.
In the text