=Paper=
{{Paper
|id=Vol-2331/paper3
|storemode=property
|title=Dynamic Simulation Model for an Autonomous Sailboat
|pdfUrl=https://ceur-ws.org/Vol-2331/paper3.pdf
|volume=Vol-2331
|authors=Moritz Bühler,Carsten Heinz,Simon Kohaut,Helmi Abrougui,Samir Nejim
}}
==Dynamic Simulation Model for an Autonomous Sailboat==
https://ceur-ws.org/Vol-2331/paper3.pdf
Dynamic Simulation Model for an Autonomous Sailboat
Moritz C. Buehler Carsten Heinz Simon Kohaut
Sailing Team Darmstadt e.V.
{moritz.buehler, carsten.heinz, simon.kohaut}@sailingteam.tu-darmstadt.de
Abstract
Sailing without any human intervention generates a great fascination,
since it is challenging while enabling many opportunities. Moving with-
out external energy supply, possibly transporting goods, collecting plas-
tic waste or recording scienti c measurement data, new attractive sce-
narios become possible. As rst milestone, we aim to send a robotic
sailboat across the Atlantic ocean, coping with bad weather including
storms, other vessels or oating waste.
For testing, designing control, sophisticated route planning, and man-
aging algorithms, we are interested in dierential equations for a sim-
ulation model. Therefore, we theoretically derive relations, describing
the sailboat motion as those of a rigid body having six degrees of free-
dom (6 DOF). This allows us to reproduce signi cant nonlinear eects
like the ones created by ow separation, speed dependence of the dy-
namics and oscillations by waves, while maintaining a comprehensible
structure of the external forces of the sailboat. In contrast to standard
velocity prediction programs (VPP), we are interested in the actual
dynamical behavior, the eect of sail angles, rudder and waves while
the precise prediction of the actual reachable velocities are of minor
interest.
The dynamic model can serve for controller design and for the genera-
tion of test data. Additionally, it is used to feed a hardware model of
our boat, providing an intuitive way of demonstrating the boat move-
ments.
1 Introduction
For large parts of the oceans, accurate measurements of the environment are missing. More detailed data would
bene t climate research and environment monitoring. A way to perform the ocean sensing is the usage of
autonomous sailing boats. This reduces the cost of the required hardware from a big research vessel to a small
boat and eliminates personnel costs as there is no crew on the ship. Further use cases are the collection of trash
and cheap, energy e cient transportation of goods.
These long term missions through rough weather conditions require a high con dence in the hardware of the
boat and in the algorithms. Therefore, a way to test the implemented controller is required. In this paper,
we present a simulator to provide a basis for testing our developed algorithms and systems of an autonomous
Copyright c by the paper's authors. Copying permitted for private and academic purposes.
In: S. M. Schillai, N. Townsend (eds.): Proceedings of the International Robotic Sailing Conference 2018, Southampton, United
Kingdom, 31-08-2018
31
ROBOTIC SAILING 2018
sailboat. This allows a short development cycle with fast feedback without the need to use a real boat. The
focus is on a six degrees of freedom (6 DOF) model for a sailing boat. It is used for simulating the dynamic
boat movement. Further it is important for controller design and a test pattern generator for the whole software
architecture as well as for visualization purposes.
1.1 Related work
For modeling the dynamics of motorized vessels there is rich literature, e.g. (Fossen, 2005), however for wind
propulsion, special eects have to be considered. The speci cs of sailing are treated in models for yachts (Philpott
et al., 1993) in the development of Velocity Prediction Programs (VPP). However, they aim to optimize the race
performance and therefore focus on the achievable speed, while we are interested in the dynamic motion eects.
Dynamic models for sailing vessels are derived in (Saoud et al., 2013), (Alves, 2010), (Masuyama and Fukasawa,
2011) and (Roncin and Kobus, 2004). However, (Masuyama and Fukasawa, 2011) concentrates on tacking,
(Saoud et al., 2013) and (Alves, 2010) reduce the model to 3 DOF, neglecting the in uence of roll and waves.
(Roncin and Kobus, 2004) relies on the identi cation of a larger set of parameters based on experiments for which
we would need a setup respectively CFD model. In contrast, we outline a 6 DOF model including wave eects
while transparently formulating the acting forces using fewer parameters to obtain a rich but comprehensible
simulation setting.
1.2 Assumptions
For our model, we take several assumptions. First, we neglect the change of the wind over the height and assume
a constant wind speed as well as angle over the whole sail. A similar assumption is made for the lateral parts
and the water speed. Also the interference between both sails is neglected. Each sail will change the wind locally
which will have eects on the other sails. Instead, lift and drag coe cients of the sails are approximated by
the results of thin airfoil theory, which additionally neglects the shape change of the sail. The ow separation
is estimated in a simple model interpolating between unseparated and fully separated ow. For the waves, it is
assumed, that the wave length is large compared to the ship length.
2 Mathematical model of the boat movement
To simulate the movement of the sailboat, we need a model in the form of dierential equations, describing the
temporal relations of external forces acting on the sailboat through dynamic and kinematic relations. For our use
case, we are interested in a rich but comprehensible model, covering important eects of the sailboat's motion.
We model the dynamics and kinematics for a 6 DOF boat model, that depend on the geometries and inertia. A
sailboat operates at the boundary between two media, water and air and uses speed dierences to move forward.
As important forces, we derive equations for the following components:
• sail, keel and rudder force
• buoyancy forces including eects from ocean waves
• wave resistance
• damping forces
2.1 Dynamics and kinematics
The dynamics describe how the forces acts on the boat speeds and turn rates due to the boat inertia. It is wise
to derive the relations of forces and dynamics in dierent reference frames, that are connected through kinematic
relations, see also Figure 1.
We describe the position xg , yg , zg and heading ψ of the boat in our globally xed or navigational frame, with
z = 0 at the undisturbed water surface.
At the center of gravity of the boat, we locate the second, the heading frame, also parallel to the water surface
with x axis in heading direction. The transformation from the global frame consists of the translation of position
(xg , yg , zg ) and rotation according to heading ψ . In this heading frame, we formulate the buoyancy Fhs , the
translation dynamics (change of speeds vx , vy , vz ) and the wave resistance Fwr . The basis vectors in the heading
frame are denoted by {ex , ey , ez }.
32
Vessel Development & Modelling
−ϕ
Fae,1
αae
γ y αhy
Fae,1
x vhy Fhs
ψ
x3 z
vae y
x2
yg Fhy
Fhy
xg mg
(a) (b)
Figure 1: Top view (a) and rear view (b): Forces of sail(s) Fae , keel Fhy and buoyancy Fhs . Position (xg , yg , z)
and heading ψ is measured in a global frame, buoyancy and lateral dynamics in a frame centered at the ship but
oriented to the water surface, and the angular dynamics in a ship oriented coordinate frame. The aerodynamic
angle of attack αae results from sail angle γ and apparent wind direction βW A = γ + αae , hydrodynamic αhy
from the leeway drift and water stream.
The third frame is the body frame, xed to the boat main axes. Here, it is intuitive to describe the rotational
dynamics, sail, rudder and keel forces. The transformation from the heading frame is achieved by rotating due
to roll angle ϕ, while the pitch angle θ is often negligible small. The basis vectors of this frame, we denote as
{e1 , e2 , e3 }.
Altogether, we have 12 states, the translations xg , yg , zg , heading ψ , roll ϕ and pitch θ and the change rates
in positions vx , vy , vz and angles p, q, r, where v = vx ex + vy ey + vz ez is the speed at center of gravity and
Ω = pe1 + qe2 + re3 the angular velocity.
From the kinematics, we obtain the rst part of dierential equations, the position (ground based) changes
according to heading and speeds,
xg vx cos(ψ) − vy sin(ψ)
d
yg = vx sin(ψ) + vy cos(ψ) . (1)
dt
z vz
Similar holds for the angles, where it is used, that the pitch angle is small.
ϕ p
d
θ ≈ cos(ϕ)q − sin(ϕ)r (2)
dt
ψ cos(ϕ)r + sin(ϕ)q
The speed and angular rates change due to forces and moments according to the dynamic equations. The boat
33
ROBOTIC SAILING 2018
speed change due to overall applied forces F and boat mass m,
d
m v = F. (3)
dt
For the left side we need to respect the angular velocity of the heading frame, the change of heading, m dtd
v=
m((v̇x − ψ̇vy )ex + (v̇y + ψ̇vx )ey + v̇z ez ).
When the boat accelerates, surrounding water is accelerated, resulting in a transient hydrodynamic force that
is usually modeled as an added mass. Additionally, this eect introduces cross terms, coupling roll and sway. We
plan to include added masses soon in our description, since they have an important in uence, especially for roll
motion (Korotkin, 2009). Therefore, (3) and (4) have to be slightly adapted to include the coupling of inertia.
For the angular velocities, we have the relation of the angular momentum ΘΩ and overall torque T,
d X
(ΘΩ) = T = ri × Fi , (4)
dt i
where Θ is the angular inertia, Ω the angular rate vector and ri the vector from center of gravity to the ith force
point of application. Due to approximate symmetries of the boat, the principal inertia axes do coincide with the
boat main axes, leading to a diagonal inertia Θ (apart from added masses). For the derivative, it is important
to respect, that the ship itself is moving and the ship based coordinate frame itself has to be dierentiated,
dt (ΘΩ) = ṗΘ1 e1 + q̇Θ2 e2 + ṙΘ3 e3 + Ω × ΘΩ. The formulation in ship based coordinates is however necessary
d
to avoid a full and time dependent inertia matrix.
In the following we derive equations for the force components describing the motion of the sailboat, needed
for equations (3) and (4).
For the aero and hydrodynamics, we have to dier apparent wind (stream) and the true wind (stream). The
true wind is the wind referred to the ground, that is equal for all vessels at a given position. For the hydro-
and aerodynamic forces, the relative ow is relevant, meaning that the ship speed has to be subtracted from the
absolute. Hence the apparent or relative wind (and water ow) is calculated as the vectorial sum of the true
wind (water stream) and the negative boat velocity. With raising boat velocity, the apparent wind changes its
directions towards coming from ahead with varying strength. Additionally, the rotation of the boat changes the
relative ow on the foils. Therefore, we chose the point of load rf oil of the foil and assume constant ow speed
over the whole foil (which seems to be ne for not too large rotation rate). The relative velocity is calculated as
vrel = vf low − v − Ω × rf oil . For the rudder this might be problematic due to the nearby keel but this remains
for further re nement.
2.2 Sails and lateral foils
A main impact results from the sails and lateral plans, keel and rudder. Both, we model analog as ideal foils in
a uid stream. We calculate the resulting uid forces based on results from theoretical aerodynamics, providing
equations for lift and drag forces and the points of attack.
In dierence to a hard foil where the actual shape is xed, the shape of a sail varies according to load and
design. However, we want to stay at a simple description combined with fast computation and neglect these
eects.
2.2.1 Lift and drag
Lateral foils and sails are placed in a uid stream that reaches them under a speci c angle of attack α. For the
keel, this is the relative angle of the water ow. For the sails, the apparent wind angle is referred to the sail
angle γ (see Figure 1).
For the forces, we use the results from theoretical aerodynamics of thin foils leading to the relation for the lift
coe cient as
cL = cL,0 + cL,α α, (5)
where cL,0 is the lift coe cient at zero angle of attack that is zero for a symmetric foil, and for thin foils a
theoretic value of the slope coe cient of cL,α = 2π is known (while neglecting the foil shape) (Spurk and Aksel,
2010).
34
Vessel Development & Modelling
The drag force of a foil consists of dierent components: Friction at the surface, induced drag due to nite
foil height and pressure loss when it comes to ow separation that we treat later. At small angles of attack, the
drag coe cient is
cD = cf + cD,i . (6)
The induced drag depends quadratically on the lift coe cient cD,i = c2L /(πΛ) with the geometric foil
√ stretching
Λ, height h and area A. The theoretical friction coe cient of a plate in laminar ow is cf = 2.66/ Rel (Spurk
and Aksel, 2010), where Rel = vl ν is the Reynolds number with length l of the plate in direction of the ow and
ν is the kinematic uid viscosity.
With drag and lift coe cients we can calculate the forces by multiplying the dynamic pressure q = 0.5ρvrel 2
and the area A,
Ff low = qA [(sin(βW A )cL − cos(βW A )cD )e1 − (sin(βW A )cD + cos(βW A )cL )e2 ] , (7)
that is transformed to the ship based frame according to the relative ow angle βW A .
The point where the aerodynamic force acts, is for the at symmetric wing at one forth length, but varies until
one third for asymmetric shapes where the location depends on the angle of attack (Spurk and Aksel, 2010).
2.2.2 Flow separation
At larger angles of attack (e.g. sailing in front of the wind for the sail, departing from low speeds for the keel),
the ow separates from the foil, the lift reduces and the drag by the pressure loss gains signi cant in uence.
In the extreme case of an angle of attack of 90 degrees, there is no lift (due to symmetries), while the drag
coe cient is slightly above cD ≈ 1 and the point of attack is the center of the foil. At smaller angles, we model
the drag coe cient as cD = sin(α)2 representing the dynamic pressure of the ow perpendicular to the foil. For
most foil shapes, ow separation leads to stall angles, the angle where cL is maximal, between 15 and 25 degree.
We interpolate
between not separated case and ow separation heuristically according to a separation factor
s = 1 − exp −(α/αsep )2 and use a characteristic angle αsep = 25◦ , resulting in a stall angle of 15◦ .
2.3 Wave resistance
The velocity of a displacement vessel is limited by the wave resistance at maximum hull speed.√This resistance
results from the waves generated by the ship movement. At about a Froude number of F r = v/ glwl = 0.4, the
wave resistance reaches a local maximum that cannot be exceeded by our keel boat. √
The maximum hull speed depends on the square root of the waterline length lwl of the ship vhull = 0.4 glwl ,
resulting to a maximum speed for a 4 m water line length of about 2.5 m/s.
We approximate the wave drag by a 6th order polynomial introducing it as speed limiting factor, neglecting
the hull shape dependent interference eects and model the wave resistance according to
4
vhy
Fwr = − sign(vx )cwr qhy ALK e1 , (8)
vhull
with lateral area ALK and a wave resistance weight parameter cwr .
2.4 Hydrostatic buoyancy and wave in uence
The hydrostatic force is central for the boat to oat on the water and the moment stabilizes roll and pitch angles.
The hydrostatic bouyancy respectively Froude-Krylov force (Fossen, 2005) is the integrated water pressure
over the hull, also described as the weight of displaced water. For the calculation of the pressure, the virtual water
surface without the vessel has to be used. In calm water and zero speed, the hydrostatic force acts vertically
(along zg ), compensating the vessels weight at the stationary point of oating. With some elevation of the ship
zg or similarly an actual wave elevation η , the hydrostatic force Fhs is calculated as
Fhs ≈ mg + ρhy gAW (η − zg ), (9)
with gravity acceleration g , water density ρhy and the water plane area at equilibrium AW . The integration of
the waterline surface changes is neglected leading to the above linear relation.
35
ROBOTIC SAILING 2018
The point where the force acts depends on the height dierence of the buoyancy point above the center of
gravity hhs,0 and the second moments of water plane area IL/T .
For the side shift, we get
ρ ρ
hy hy
rhs ≈ − IL + hhs,0 sin(ϕef f )ey + IT + hhs,0 sin(θef f )ex , (10)
m m
again neglecting the change of the water plane due to roll and heave. The rst part of each component respects
the stabilizing in uence of the hull shape of the ship, the second the in uence of the static deviation from the
center of gravity. The second moment of water plane area IT for the pitch movement is at a larger order than
IL for the roll movement due to the lengthiness of common boat shapes. Therefore, pitch angles stay small
especially in the absence of waves.
The eective hydrostatic roll angle has to be referred
to the water surface including
waves. With the wave
elevation eld η(x, t) it becomes ϕef f = ϕ − arctan ∂η
∂y and θ ef f = θ + arctan ∂η
∂x , where it is assumed, that
the wave length is large compared to the ship length (4 m), which should be the case for oshore ocean waves
with signi cant amplitudes.
The buoyancy force changes the direction according to the water surface, and partially acts in horizontal
∂η
directions, again linearized to Fhs = Fhs ez − Fhs ∂x ex − Fhs ∂η
∂y ey .
2.5 Damping
Lastly, we model linear constant damping of the boat, though the addition of a quadratic term will be more
appropriate. The used damping force and moment are FD TD = −d v Ω ,where we use decoupled, rough
estimates for the damping parameters d.
2.6 Inputs
For a simulation, we have to specify all external variables like wind and waves as well as control variables, rudder
and sail angles. The rst are part of the simulation environment, while the latter have to be chosen by a boat
controller.
For the rudder, a limitation (±35◦ ) is used. For the choice of the rudder angle, the ow separation at high
angles of attack of the rudder has to be avoided or respected during controller design. For the sails, the possible
angle is geometrically restricted by shrouds (±90◦ ). A stronger constraint is imposed by the limited rate of
change for which we introduce a rst order model in the next part. Additionally, for practical reasons of energy
management, the sail angle should not be changed too often. An adaption seems to be necessary mainly at
maneuvers and changes in heading and wind. In contrast, the rudder is important to control the heading and
reducing the impact of disturbances like waves or wind gusts. We choose a sample time of 100 ms, but its value
may be adapted according to energy resources and control requirements.
2.6.1 Limited actuator speeds
It is not possible to set sail and rudder angles immediately, instead changing them requires time and energy.
The actuators have to change the rope length (sail) or the rudder position. To respect this in a simple way, we
chose a rst order linear behavior determined by some time constant τ describing the typical time for the rudder
respectively sail to be set. This temporal behavior results from the actuator position controllers. Additionally,
the absolute change rate may be limited due to a maximum actuator speed.
The time constant of the rudder is assumed to be short (0.1 s), while for the sail we expect a larger time
constant due to hardware design (e.g. 3 s). For the sails, the actuation by ropes implies another eect: Its angle
can change very fast during jibing or tacking. It is rather the maximum absolute angle that is chosen by the
controller setting the length of the rope. The side respectively sign of the sail angle is chosen by the relative
wind.
2.7 Implementation and numeric integration of the dierential equations
For the simulation, the dierential equations of the dynamics model (1)-(4) have to be integrated departing from
the initial state values x0 . Therefore we use the Dormand Prince method, a standard Runge-Kutta method with
adaptive step size. For a simulation time of 120 seconds, 2 seconds computation time is needed on 2 GHz single
CPU computation.
36
Vessel Development & Modelling
1.0
Forward speed in m/s
Leeway in / ms
0.8
Flow separation 10
0.6 0
−10
yg / m
0.4
−20
0.2 −30
0.0 −40
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Time /s 0 10 20 30 40 50 60 70
xg / m
Figure 2: Boat accelerating from zero speed. For-
ward speed vx and leeway vy are both rising, until Figure 3: Simulated ship trajectory: Gaining speed
the overall speed augments and the keel compen- tacking and jibing; boat starts at blue cross, wind
sates the side force. The high initial leeway drift direction is marked by the green arrow.
angle leads to ow separation.
Forward Forward
2.0
Leeway Leeway
1.5 Vertical Vertical
1.5
Speed / m/s
Speed / m/s
1.0 1.0
0.5
0.5
0.0
0.0 −0.5
0 25 50 75 100 125 150 0 25 50 75 100 125 150
Time / s Time / s
Figure 4: Speed of the boat in the dierent direc- Figure 5: Same simulation scenario but with waves:
tions: Horizontal speed is almost decoupled, hence Ship speeds for a wave eld with 1 m wave height
close to zero in the absence of waves. at a wave length of 100 m, arriving from right.
The implementation of our simulation setup is available at https://github.com/simko96/
stda-sailboat-simulator.
3 Simulation
The simulator should be able to reproduce main eects of the sailboat dynamics to deliver useful data for testing
software components like controllers or path planning. One interesting situation is the departure from zero speed,
where the eects of rudder and keel dier from normal operation (here the inertia dominates). Other interesting
situations include state changes during maneuver execution and the in uence of bigger waves. The simulation
serves as data generation to test our modules like heading controller, path planning or lter algorithms. However,
so far we do not have empirical data to validate its performance. Instead we qualitatively examine for plausibility.
3.1 Scenario: gaining speed, tacking and jibing
Our scenario consists of the boat accelerating at departure from zero speed, while the wind comes from south
with 5m/ s (right angle to initial heading), followed by tack and jibe maneuvers (Figure 3). Therefore, we use a
feedback control for the heading through the rudder, that is derived from a simpli ed nonlinear dynamics model
for the yaw motion, and a feed forward control of the sails.
The rst part of the scenario is to gain speed (Figure 2). At departure, there is no water ow which leads to
a high initial leeway (acceleration according to the sail force), and the starting ow separates. While the boat
37
ROBOTIC SAILING 2018
360 Heading 7.5 Roll
315 Reference Pitch
Drift 5.0
270
Angles / degree
Angles / degree
2.5
225
180 0.0
135 −2.5
90
−5.0
45
0 −7.5
0 25 50 75 100 125 150 0 25 50 75 100 125 150
Time / s Time / s
Figure 6: Plot of the controlled heading angle, its Figure 7: Corresponding angles in roll and pitch,
reference and the leeway drift angle compensated due to the high longitudinal stability (geometric
by the controller. inertia), the pitch angle remains small.
speeds up, the keel force raises until it compensates for the side force of the sails.
At tacking, the speed reduces due to the missing propulsion when crossing the wind. Afterwards, the boat
speeds up at the new heading. For jibing, the speed decrease is less pronounced (Figure 4). During both
maneuvers, the roll angle changes the sign. The pitch angle remains negligible during the scenario (Figure 7)
and becomes larger but still small when introducing waves (not shown). The heading (controlled by use of the
rudder) is shown in Figure 6.
With waves of 1 m height coming from east (−xg -direction), the base trajectory looks similar. Oscillations
induced by the waves are superposed on all states, exemplary shown for the speeds in Figure 5.
3.2 Visualization by hardware model
Another application is the data generation for our
movable hardware model of the boat shape that
serves for demonstration and visualization issues
(Figure 8). It is eligible when presenting our project,
e.g. at fairs, and useful to test interfaces between
the dierent components, software, electronics and
actuators.
4 Conclusion
In this paper, we derived the dierential equations
describing the dynamics of a sailboat, by considering
the forces generated by both uid ows, wind and
the relative water stream.
This model is used, to simulate the sailboats mo-
tion, particularly for maneuvers like tacking or jib- Figure 8: Demonstration model of a sailing boat.
ing. The simulation results seem plausible and well
suitable for our demonstration model. Further, it is useful for testing our control and path planning algorithm.
However, we still need to record real data, to identify parameters like damping and compare and evaluate the
simulation quality.
References
Alves, C. (2010). Sailbot: Autonomous marine robot of eolic propulsion. Master's thesis, Instituto Superior
Tecnico, Universidade Tecnica de Lisboa.
Fossen, T. (2005). A nonlinear uni ed state-space model for ship maneuvering and control in a seaway. Journal
of Bifurcation and Chaos, 15:27172746.
Korotkin, A. (2009). Added Masses of Ship Structures. Springer.
38
Vessel Development & Modelling
Masuyama, Y. and Fukasawa, T. (2011). Tacking simulation of sailing yachts with new model of aerodynamic
force variation during tacking maneuver. Journal of Sailboat Technology, 119.
Philpott, A. B., Sullivan, R. M., and Jackson, P. S. (1993). Yacht velocity prediction using mathematical
programming. European Journal of Operational Research, 67:13 24.
Roncin, K. and Kobus, J. (2004). Dynamic simulation of two sailing boats in match racing. Sports Engineering,
7:139152.
Saoud, H., Hua, M., Plumet, F., and Amar, F. (2013). Modeling and control design of a robotic sailboat. Robotic
Sailing.
Spurk, J. and Aksel, N. (2010). Stroemungslehre. Springer.
39