674 Journal of Marine Science and Technology, Vol. 21, No. 6, pp.

674-685 (2013 )
DOI: 10.6119/JMST-012-1121-1

PHYSICAL MODELS FOR SIMULATING SHIP

STABILITY AND HYDROSTATIC MOTIONS

Shyh-Kuang Ueng

Key words: ship motion simulator, physical-based animation, ship paper, we propose efficient physical models for simulating
stability simulation, physics engine. hydrostatic ship motions caused by sea waves, cargoes and
flooding water. Our physical models are based on the stability
theory of ship. We widen the scope of conventional ship mo-
ABSTRACT
tion simulation techniques by using the gravity, buoyancy and
In this paper, efficient physics models are presented for simu- meta centers of the ship to estimate the stability of the ship and
lating ship hydrostatic motions caused by sea waves, cargoes compute the forces and torques acting on the ship so that the
and flooding water. Based on the stability theory of ship, a ship motions can be simulated in real time.
ship is regarded as a floating body whose stability is indicated Since the ship shape is irregular, it is impossible to compute
by its gravity, buoyancy and meta centers. These stability the stability centers analytically. To overcome this problem,
centers are influenced by cargoes, sea waves, flooding water we split the ship body into cells and distribute the ship mass,
and the ship mass. Variations of these stability centers create cargoes and flooding water into the cells. Thus, the stability
torques and forces which cause the ship to heave, pitch and centers, forces and torques can be computed by using finite
roll. Because of the irregularity of the ship shape, analytical volume methods. Furthermore three distinct coordinate sys-
solutions of these stability centers do not exist. In this work, a tems are employed in the proposed models to specify the po-
finite volume approach is employed to overcome this problem. sition and orientation of ship, the geometry of ship body, the
At first, we split the ship body into cells by using a regular grid external forces and ship motions so that the formulations of
and distribute the ship mass, cargoes and flooding water into forces and ship motions are simplified and all essential entities
these cells. Then numerical methods and physics laws of can be efficiently calculated by using fundamental physical
floating body are utilized to compute the coordinates of these laws. The proposed models are efficient, realistic and flexible.
centers. Subsequently, torques and forces are calculated and Users are allowed to compose various simulation scenarios
utilized to create ship motions. by loading cargoes into the ship, generating severe sea waves
and creating holes on the ship hull. The resulted ship motions
will be automatically animated in real time with decent visual
I. INTRODUCTION
effects.
Ships have been utilized for fishing, transportation, recrea- The rests of this paper are organized as follows: Related
tion and adventures for thousands of years. Learning the researches on ship motion simulation are presented in Section
knowledge of ship-handling is essential to our daily lives. In II. In Sections III-VI, the proposed physical models are sys-
modern maritime training courses, ship motion simulators are tematically formulated. To establish the theoretical back-
widely used to teach people how to operate ships under dif- ground, the stability theory of ship is introduced in Section III.
ferent sea conditions. To become an effective learning tool, a Then the coordinate systems and the wave model employed in
ship motion simulator must be equipped with a physics engine the proposed ship motion models are described in Section IV.
which can simulate the maneuverability of ship as real as The computational grid and the algorithms for computing the
possible under a severe time constraint [10]. Scientists have stability centers are presented in Section V. The procedures
developed some mathematical models for computing ship for calculating ship motions are formulated in Section VI.
motions [1, 9, 16]. However, realizing these models requires Experimental results are presented and discussed in Section
tremendous computational costs, and thus these model are too VII. The conclusion of this paper is drawn in the last section.
slow for building real-time ship motion simulators. In this
II. RELATED WORK
Paper submitted 01/04/12; revised 09/07/12; accepted 11/21/12. Author for
correspondence: Shyh-Kuang Ueng (e-mail: skueng@mail.sju.edu.tw).
When floating on sea surface, a ship is surrounded by two
Department of Computer Science and Engineering, National Taiwan Ocean types of fluids, the air and the sea water. The material prop-
University, Keelung, Taiwan, R.O.C. erties of air and water are different, and the boundary between
 S.-K. Ueng: Physical Models of Ship Stability 675

the air and the sea water may fluctuate with time. Therefore use the Navier-Stokes equations to compute liquid motions.
ship motions are difficult to model. Furthermore, the ship They compute the time varying fluid surface based on the
shape, draft and displacement also affect the motions of the marked particles approach. Once the velocity field and the
ship. These factors make the modeling of ship motion even pressure field have been calculated, the forces acting on the
harder. Currently, most numerical algorithms for computing ship floating objects are estimated. Subsequently, the motions of
motions are based on the strip theory, developed in [9, 16]. the floating objects are calculated by using the Lagrange
These algorithms calculate ship motions by solving highly equations of motions. In [19], a numerical model is pro-
nonlinear partial differential equations. Users may spend posed to generate waves. Then the drag, lift and buoyancy
hours of CPU time just to obtain a set of solutions for one or forces acting on the floating bodies are computed and used
two ship motions [1]. These traditional numerical procedures to estimate the interactions between waves and the floating
are too slow for constructing ship motion simulators. bodies. In [3], Carlson et al. present the Rigid Fluid method
Some researchers proposed simplified governing equa- to compute the coupling between rigid bodies and fluids.
tions and boundary conditions to model ship motions. In [20], They regard rigid bodies as fluids. Thus the motions of
Zhang et al. developed mathematical models for simulating fluids and objects can be modeled by using the Navier-Stokes
the motions of a ship sailing inside a harbor. In their method, equations. However, the rigid body motion constraints are
the forces acting on the ship are estimated first. Then, sim- enforced in the governing equations such that rigid objects
plified differential equations are derived to model the relations are not allowed to deform. In [8], Kim et al. apply the diver-
between the forces and the accelerations of ship motions. gence theorem to convert volume integration into surface
Subsequently, the differential equations are solved by using integration and use graphical hardware to speed up the com-
a Runge-Kutta method to compute the ship motions. Since putations. Then they utilize their hardware integration pro-
their models focus on the ship’s motions inside harbors, only cedure to compute the gravity forces and buoyancy forces
the physical models for surge, sway and yaw are developed. acting on floating bodies. In turn, the motions of floating
Another simple ship motion model is reported in [4]. The objects are calculated by using these forces.
authors use the sea wave model proposed in [7] to generate In [18], six mathematical models are proposed by Ueng
sea waves. At each time step, the height field of the sea sur- et al. to simulate ship motions. In their method, the height
face under the ship body is computed. Then, the tangent field of the sea surface under the ship body is computed to
plane of the sea surface is calculated and the ship is rotated estimate the external force acting on the ship. Then the ex-
such that its orientation is aligned with the tangent plane. ternal force is reduced by the damping effect of sea water
Their method can simulate heave, pitch and roll of a ship. to produce the net external force. The magnitudes of heave,
However, sea waves are the sole force of the ship motions. pitch and roll are computed by using the net external force.
The ship will stop moving immediately as soon as sea waves They also proposed mathematical models to compute the
come to rest. This phenomenon is in conflict with the be- gross internal force produced by the ship’s propellers and
haviors of a real ship. rudders. Then the resistance of sea water is subtracted from
Other researchers propose statistical methods to predict ship the gross internal force to generate the net internal force. In
motions. In [17], Kalman filters are utilized to estimate ship the following step, the magnitudes of surge, sway and yaw
motions. Their models rely on many ship parameters to tune are computed by using the net internal force. Finally, the
the governing equations. If these parameters are absent, the six motions are super-imposed to generate the resulted ship
simulation cannot be performed. In [11], an improved method motions. Their method allows users to tune ship motions
for computing ship motions is proposed. Their method re- according to the ship shape, draft and mass. Therefore ships
quires less ship data. However, the accuracy of their model of different shapes can behave differently under the same
is still dependent on the availability of ship parameters. In sea condition. The simulation results produced by using their
ship motion simulations, the ship models are usually fabric- models are physically sound and visually realistic.
cated. Real ship parameters are hard to estimate. In [21], However, all the aforementioned ship motion simulation
another statistical method is proposed to predict ship motions. methods are not based on the stability theory of ship. Ship
In their method, the status and motions of the ship at some motions caused by unbalanced cargo loading, body damages
previous steps are recorded and converted into a tensor field. and severe sea conditions are not or only partially modeled.
Then the eigenvalues and eigenvectors of the tensor field are In their models, a ship can always retain its initial stability
computed. Finally, a minor component analysis is conducted under any circumstances. It will not list, overturn and sink
to predict the ship motion at the next time step. Their method even it is loaded with heavy cargoes, encounters severe sea
is useful for short-term ship motion prediction. For a long conditions or bears serious ship hull damages. Therefore,
term ship motion simulation, the accuracy of their method may these models are in conflict with the real hydrostatic behav-
be decreased by accumulated numerical errors. iors of ships. In this article, a new set of physical models are
Recently, in the computer graphics society, several CFD proposed to supplement these conventional ship motion
algorithms have been proposed to model the interactions be- models such that the hydrostatic behaviors of ships can be
tween fluids and floating bodies. In [6], Foster and Metaxas better simulated.
676 Journal of Marine Science and Technology, Vol. 21, No. 6 (2013 )

gravity
M force
gravity G gravity
waterline gravity Z
force buoyancy force
G force M
G Z force
G M
B B
buoyancy B1 M buoyancy
force G Z force
B B B1
B B1 B1
vertical line vertical line buoyancy
(a) a ship at rest (b) an linclined ship force
(a) righting lever (b) capsizing lever (c) neutral lever
Fig. 1. The gravity center G, buoyancy center B and meta center M of a
ship. (a) When the ship rests on a flat water surface, G and B are Fig. 2. Three different levers created by the gravity and buoyancy forces,
collinear. (b) As the ship inclines, B1 becomes the buoyancy cen- (a) G is under M, a righting lever, (b) G is above M, a capsizing
ter and M is located at the intersection of the line passing B and G lever, (c) G = M, a neutral lever.
and the vertical line passing B1.

G, the ship is in an unstable equilibrium. It will capsize to

gain a new stability. On occasions, M coincides with G and
III. THE STABILITY THEORY OF SHIP
the ship is in a neutral equilibrium. Any external force will
Our physical models are based on the stability theory of change its posture forever.
ship. The fundamental knowledge of ship stability is briefly
presented in this section. Detail description about this topic 2. The Levers of Stability
can be found in the books written by Lewis [12-14] and These three equilibriums are illustrated in Fig. 2. Assume
Derrett [5]. that the perpendicular projection of G in the vertical line
passing through B1 is Z. Then the edge GZ serves as a lever
1. The Stability Centers which creates a torque to change the posture of the ship. In
The stability of a ship is decided by its gravity, buoyancy the part (a), M is above G and the ship is in a stable equilib-
and meta centers [5]. The gravity center G is the point at rium. The edge GZ is a righting lever which tends to bring
which the mass of the ship is concentrated. The gravity force the ship back to its original posture. In the part (b), M is below
is considered to act downward at G. The magnitude of the G and the ship is in an unstable equilibrium. The edge GZ is
gravity force is equal to the weight of the ship. The buoyancy a capsizing lever which makes the ship rotate further until the
center B is the geometrical center of the underwater portion of ship obtains a new balance. A neutral equilibrium is shown in
the ship. It is also the geometrical center of the displaced the part (c). In this case, G is coincident with M. Any force
water body. The buoyancy force acts vertically upward at B. will change the posture of the ship.
According to Archimedes principle, the magnitude of the
buoyancy force is equal to the weight of the displaced water
body.
IV. THE COORDINATE SYSTEMS AND
When the ship floats at rest on flat sea surface, the gravity
WAVE MODEL
center and the buoyancy center are located at the same vertical In a ship motion simulation, many variables have to be
line [5], as shown in the part (a) of Fig. 1. The magnitudes updated at each time step, for examples, the orientation and
of the gravity force and the buoyancy force are equal, but position of the ship, the three stability centers, the forces and
these two forces act in opposite directions. If the magni- torques and the ship motions. Since the ship may vary its ges-
tudes of these two forces are different, an external force is ture and position constantly, specifying and computing all
generated and the ship will heave upward or downward. If these variables in the world coordinate system will complicate
the ship inclines or the water surface varies, the buoyancy the simulation. In the proposed physical models, three distinct
center will move to a new position B1, as depicted in the part coordinate systems are employed to define these entities such
(b) of Fig. 1. Since the gravity force and the buoyancy force that essential computations are localized and simplified. Sea
no longer act in the same vertical line, a rotational torque is waves produce forces and vary the position, posture and ori-
created and causes the ship to roll or pitch. entation of the ship. Sea waves are an important cause of ship
The vertical line passing the new buoyancy center B1 will motions. The wave model used in our work is also presented
intersect the line passing B and G. The intersection point is in this section.
called the meta center M of the ship. The meta center is the
indication of stable equilibrium. If M is above G, the ship is 1. The Coordinate Systems
in a stable equilibrium. The ship will retain its posture after The three coordinates systems used in the proposed physi-
making successive motions. On the other hand, if M is below cal models are shown in Fig. 3. The first coordinate system is
 S.-K. Ueng: Physical Models of Ship Stability 677

Y Xb
Ys Yb
Heave
X
O Yaw Roll Xb
Z Xs

Yb
Surge P

Xb G Zb P Zb
Pitch (a) bounding box (b) top view of the grid
P Yb Yb
Zb Zs
Sway
Fig. 3. Three coordinate systems and the six ship motions.

P Zb P Xb
(c) rear view of the grid (d) side view of the grid
the world coordinate system which is used to designate the
Fig. 4. The computational grid. (a) The bounding box, (b) top-view of
ship’s position and orientation and model sea waves. In the
the grid, (c) rear view of the grid, and (d) side view of the grid.
world coordinate system, the sea surface is spanned by the
X and Z axes. The Y axis is pointed vertically to the sky. The
second coordinate system is the body coordinate system. This where Ai, λi and ωi are the amplitude, wave length and speed
coordinate system is attached to the bounding box of the
of the i-th sinusoidal wave. The variable θi is the angle be-
ship body and moved with the ship. The three axes, Xb, Yb
tween the X axis of the world coordinate system and the in-
and Zb, are aligned with the bounding box boundaries. This coming direction of the i-th sinusoidal wave, and t is the
coordinate system is used to define the geometrical mesh of
time variable. We employ this wave model to generate sea
the ship body and the positions of the three stability centers.
waves, because it offers straight forward mechanisms for
The third coordinate system is the sea-keeping coordinate generating various sea waves. Our ship physical models are
system. It is used to specify and compute the forces and ship
independent of this wave model. If it is necessary, this wave
motions. Its origin is the gravity center G. The three axes of
model can be replaced by other wave models.
this coordinate systems are Xs, Ys and Zs. The Xs axis is aligned
with the forward direction of the ship. The Zs axis is hori-
zontally pointed to the starboard (right side) of the ship, and V. COMPUTING THE STABILITY CENTERS
the Ys axis is directed vertically upward. In this coordinate The three stability centers affect the ship’s stability. Varia-
system, heave is the translate motion along the Ys axis, pitch
tions of these stability centers generate forces and torques and
is the rotation about the Zs axis, and roll is the rotation about
trigger ship motions. In this section, the numerical methods
the Xs axis. When the ship rests on a flat water surface, the for calculating the gravity centers are presented. Since the
axes of the body and sea-keeping coordinate systems are par-
ship changes its posture and position frequently, specifying
allel, but the two origins are separated by a fixed distance.
these centers in the world coordinate system is a difficult job.
However, as the ship moves, the axes of these two coordinate Thus they are computed in the body coordinate system.
systems will not coincide. For example, when the ship in-
clines, so does the bounding box. The Xb is rotated about the 1. The Computational Grid
Zs axis, but the axis Xs is not changed. The axis Xs always Since the ship shape is irregular, the stability centers cannot
points to the forward direction of the ship and is parallel to
be computed analytically. In this work, they are calculated by
the flat water surface.
using a finite volume method. At first, the ship body is split
2. The Wave Model into cells by using a regular grid, as shown in Fig. 4. The size
of a cell is 1meter × 1meter × 1meter (= 1 cubic meter). The
Sea waves are usually modeled as a time-dependent height
equivalent volume of water weighs about 1 ton. Then those
field. A widely used sea wave model is proposed in the papers cells which are outside the ship body are excluded. To do so,
of [7, 15]. To increase the fidelity of sea waves, we assume
we associate each cell with a flag. If a cell is fully or partially
that sea waves are composed of several sinusoidal waves of
contained in the ship body, its flag is set to 1. Otherwise its
different frequencies, amplitudes, speeds and directions [18]. flag is set to 0 and this cell is excluded from the following
In our work, sea waves are defined by:
computations. A cell with a flag value of 1 is called a sig-
nificant cell.
n
2π
h ( x, z , t ) = ∑ Ai sin( (( x cos θi − z sin θi ) − ωi t )), (1) Once the computational grid has been built, the ship mass,
i =1 λi cargoes and flooding water are distributed to all the significant
678 Journal of Marine Science and Technology, Vol. 21, No. 6 (2013 )

Xb cell centers

Yb
mk
Yb
mi

Xi sea
surface
Zb Zk Zb
Xb
(a) beam structure in Xb direction (b) beam structure in Zb direction
Fig. 6. The immersed portion of a cell is decided by the height of sea
Fig. 5. The ship body is divided into cross-sections and treated as beam surface and the position of the cell center. The immersed volume
structures when computing G. of the ship body is equal to the sum of the immersed portions of
the cells.

cells via a graphical user interface. In the ship motion simu-

lation, users can add cargoes and flooding water into signifi- vertical position of the cell center. If the cell center is lower
cant cells to alter the stability of the ship. Users can also re-
than the sea surface by more than one half of the cell width, the
move the flooding water or cargoes from significant cells to
cell is totally immersed. If the cell center is higher than the
reduce the mass of the ship. sea surface by more than one half of the cell width, the cell is
2. Computing the Gravity Center G entirely above the sea surface and its immersed volume is zero.
Otherwise, the volume of the immersed portion of the cell is
To compute the Xb coordinate of G, we group the signifi- computed by:
cant cells with the same Xb coordinate into a slice and divide
the ship body into L cross-sections in the Xb direction, as
vim = (( H sea − C y ) / Cw ) + 0.5)vc , (3)
shown in the part (a) of Fig. 5. Then the ship body is treated
as an 1-D beam structure. By applying the physical law for
computing the gravity center of a beam, the Xb coordinate of where vim is the immersed volume of the cell, Hsea is the height
G are calculated by: of the sea surface computed by using Eq. (1), Cy is the ver-
tical position of the cell center, Cw is the width of the cell, and
L vc is the volume of the cell.
X bG = ( ∑ mi xi ) / M s , (2) The method of computing the immersed volume of the
i =1 ship body is illustrated in Fig. 6. The ship are divided into
cells. The cells shaded with the dark gray color are totally
where X bG is the Xb coordinate of G, mi is the mass of the i-th immersed while the cells shaded with light gray colors are
cross-section, xi is the Xb coordinate of the i-th cross-section partially immersed.
center, and Ms is the ship mass. The value of mi is calculated Once the immersed volumes of all significant cells have
by summing up the masses of the significant cells of the i-th been calculated, the ship body is divided into L cross-sections
cross-section. The same approach is adopted to compute the in the Xb direction. The Xb coordinate of B is computed by:
Zb coordinate of G. In this circumstance, the ship is divided
X bB = ( ∑ i =1 ρVi xi ) / M s ,
L
into cross-sections along the Zb axis, as shown in the part (b) (4)
of Fig. 5. The Yb coordinate of G is calculated by following
this procedure too.
where X bB is the Xb coordinate of B, ρ is the density of sea
3. Computing the Buoyancy Center B water, and Vi is the volume of the immersed portion of the i-th
The buoyancy center B is the geometrical center of the cross-section. Vi is computed by summing up the immersed
immersed portion of the ship body. If the immersed portion portions of the significant cells of the i-th cross-section. The
has been figured out, the algorithm for computing G can be other two body coordinates of B are computed by using this
used to calculate the body coordinates of B. method too.
To estimate the immersed portion, the grid of the ship body
is vertically projected onto the X-Z plane of the world coor- 4. Computing the Meta Centers
dinate system to create a 2D grid. Assume that the ship is The meta centers of pitch and roll are at different posi-
not present, the wave model defined in Eq. (1) is employed tions. Let’s denote the two meta centers as Mp and Mr. Ac-
to generate a height field in the 2D grid. Then the ship is cording to the stability theory of ship [5], the vertical distances
brought back. The immersed portion of each significant cell between B and the two meta centers can be analytically cal-
is computed by comparing the sea surface height with the culated by:
 S.-K. Ueng: Physical Models of Ship Stability 679

BM p = (breadth * length 3 ) /(12 M s ), (5) ah = Fh / M s , (9)

v h = v h + a h ∆t , (10)
BM r = (length * breadth 3 ) /(12 M s ), (6)

where ah is the acceleration of heave and ∆t is the time step

where breadth and length are the width and length of the ship size. Finally, the magnitude of heave is calculated by:
body, and BMp and BMr are the vertical distances between
B and the meta centers of pitch and roll. Since the meta cen-

∆heave = vh ∆t , (11)
ters are always vertically above B, their positions can be
computed once B, BMp and BMr are available. heave = heave + ∆heave. (12)

VI. PHYSICAL MODELS OF SHIP MOTIONS 2. The Physical Model of Pitch

As floating on the sea surface, a ship possesses six degrees If the gravity center and the buoyancy center are not col-
of freedom, which include heave, pitch, roll, surge, sway and linear, the buoyancy force exerts a rotational torque and forces
yaw, as shown in Fig. 3. These motions can be divided into the ship to pitch and roll, as shown in the parts (a) and (b) of
two groups. The first group of motions are triggered by Fig. 2. Let Z be the perpendicular projection of the gravity
the ship’s propellers, rudders, sea currents and winds. Surge, center G in the vertical line passing the buoyancy center B1,
sway and yaw belong to this group. Through these motions, as shown in the parts (a) and (b) of Fig. 2. The edge GZ acts
the ship navigates from one place to another. The other group as the lever of rotation. By calculating the moments at G,
of motions are induced by cargoes, flooding water inside the the magnitude of the rotational torque is defined as:
ship and sea waves. This motion group includes heave, pitch
and roll. These motions vary the ship’s posture and cause the τ = Ww * GZ , (13)
ship to oscillate, capsize and sink. Thus the hydrostatic prop-
erties of the ship under different sea conditions are revealed. where Ww denotes the weight of the displaced water body and
This paper mainly focuses on the development of physical GZ is the length of the lever.
models for animating the hydrostatic behaviors of ships.
In the sea-keeping coordinate system, pitch is the rotation
Therefore, only the physical models of heave, pitch and roll
about the Zs axis and triggered by the Xs component of the
are considered. The physical models of other ship motions can
rotational torque. Thus, the force of pitch is equal to the Xs
be found in the paper of [18].
component of the rotational torque. As the ship pitches, the
1. The Physical Model of Heave resistance of sea water will prevent the ship from inclining or
tilting further. Assume that the resistance force is propor-
The difference of the gravity force magnitude and the
tional to the speed of pitch but exerts in the opposite direction.
buoyancy force magnitude causes the ship to heave. As the
Then the magnitude of the net force of pitch is determined as
ship heaves, the viscosity of sea water produces a resistance
follows:
force to keep the ship from moving upward and downward.
Assume that the resistance force is related to the moment of

heave and the ship shape. The net force for heave is estimated τ p = Ww * < GZ , X s >, (14)
by:
Fp = τ p − bp I pω p . (15)


Rh = bh M s vh , (7)
In Eq. (14), τp is the magnitude of the torque of pitch, <> is


Fh = (Ww − Ws ) − Rh . (8) the inner-product operator, X s is the directional vector of the
Xs axis, and GZ is the vector from G to Z. In Eq. (15), the
In Eq. (7), Rh is the resistance force, bh is the damping co- magnitude of the net force of pitch, Fp, is obtained by sub-

tracting the resistance from the magnitude of the torque of
efficient of heave, Ms is the ship mass, and vh is the current
heave velocity. In Eq. (8), Fh is the net force for heave, Ww pitch. The resistance is determined by the damping coefficient
is the weight of the displaced water body and is equal to the bp, the inertia of moment Ip and the current pitch speed ωp.
buoyancy force, and Ws is the ship weight which is equal to the Then the magnitudes of acceleration, speed and angle of
gravity force. The net force Fh is obtained by subtracting the pitch are computed by:
resistance force Rh from the surplus of the buoyancy force
Ww over the gravity force Ws. Once the net force has been α p = Fp / I p , (16)
obtained, the acceleration and velocity of heave are computed
and updated by: ω p = ω p + α p ∆t , (17)
680 Journal of Marine Science and Technology, Vol. 21, No. 6 (2013 )

1 α r = Fr / I r , (23)
ri 2
ωr = ωr + α r ∆t , (24)
Xb mi G
1 2 G L Xb
∆θ r = ωr ∆t , (25)
rj
Zb Zb
mj
K θ r = θ r + ∆θ r , (26)
(a) beam structure for pitch (b) beam structure for roll
Fig. 7. Beam structures for computing the inertia moments. (a) For
where the angle of roll is represented by θr.
computing the inertia moment of pitch, the ship body is divided To estimate Ir, the ship is treated as a beam with K
into cross-sections along the Xb axis. (b) The ship body is divided cross-sections as shown in the part (b) of Fig. 7. By using
into cross-sections along the Zb axis for computing the inertia the same method for computing the inertia moment of pitch,
moment of roll. Ir is approximated by:

∆θ p = ω p ∆ t , (18) I r = ∑ m j rj2 , (27)

j =1

θ p = θ p + ∆θ p , (19)
where mj is the mass of the j-th cross-section, and rj is the
difference of the Zb coordinates of G and the j-th cross-section
where αp is the angular acceleration magnitude of pitch and center.
θp is the angle of pitch.
The inertia moment Ip in computing pitch plays the same 4. The Damping Coefficients
role as the ship mass in calculating heave [2]. To compute Ip, The damping coefficients bh, bp and br are critical parame-
the ship is divided into L cross-sections and regarded as a ters in our physical models. They decide the converge rates
beam, as shown in the part (a) of Fig. 7. Based on physical and magnitudes of the ship motions. Smaller damping coeffi-
laws, the inertia moment Ip is calculated by: cients result in longer and larger oscillations, while larger
damping coefficients produce shorter and smaller oscillations.
L
In our models, these coefficients are confined to the range of
I p = ∑ mi ri2 , (20)
i =1
{0~2}, but users are allowed to modify these coefficients to
enhance or reduce the magnitudes and durations of ship mo-
tions. Some default values of these coefficients are predefined
where ri is the difference of the Xb coordinates of G and the
in our physical models for special types of ship. For example,
i-th cross-section center and mi is the mass of the i-th
the damping coefficients of a massive and streamlined war
cross-section.
ship are set to 1.3. Therefore, the war ship is relatively insen-
3. The Physical of Roll sitive to sea waves. On the other hand, these coefficients of a
raft are set to 0.5 so that the raft can fluctuate with sea waves.
The physical model of roll is similar to that of pitch. To
compute the angle of roll, the Zs component of the rotational
torque is computed first. Then the magnitude of the net force VII. EXPERIMENTAL RESULTS
of roll is obtained by subtracting the resistance of sea water AND ANALYSIS
from the magnitude of torque of roll:
Based on our physical models, a ship motion simulator is

created to animate hydrostatic behaviors of ships. The em-
τ r = Ww * < GZ , Z s >, (21) bedded computer of our simulator is a desk-top machine
equipped with an Intel Core2 CPU of 1.86 GHz clock rate and
Fr = τ r − br I rωr , (22) an Nvidia GeForce 8800 GTS GPU. The rendering subroutine
is implemented by using OpenGL libraries. The computa-
where τr and Fr are the magnitudes of the torque and the net tional modules are designed by using C-language. Some basic
force of roll, and br, Ir and ωr are the damping coefficient, and compound ship motions are simulated by using our system.
inertia moment and speed of roll. The results are presented and discussed in this section.
In turn, the net force magnitude is used to compute the The main Graphical User Interface (GUI) of the system is
acceleration magnitude of roll αr. The speed of roll ωr is shown in the part (a) of Fig. 8. The main GUI composes of
updated once the acceleration magnitude is available. The the main screen which shows the motions of the ship and
increment of roll angle is decided by multiplying the roll speed a sub-window which displays the path of the ship. The other
with the time step size: two key GUI are depicted in the parts (b) and (c). In the part
 S.-K. Ueng: Physical Models of Ship Stability 681

Force

(a) Time = 0 (b) Time = 6

(a)
(c) Time = 10 (d) Time = 90

Heave (meter)
2
Time (sec)

-2
Ti me
Time = 0

Ti
Ti

m =

m
e=6

e=
10

90
(e)
Fig. 9. (a)-(d) Snapshots of the heave motion of a ship, (e) the magnitude
of the heave motion at each time step.
(b) (c)
Fig. 8. (a) The main GUI of the simulator, (b) the window to the show
wave height and motion magnitudes, (c) the GUI to load and
unload cargoes and create body damage. immersed portion of the ship is increased, the buoyancy force
magnitude is larger than the gravity force magnitude. When
the external force is released, the buoyancy force lifts the
(b), the wave height under the gravity center and the magni- ship upward immediately, as illustrated in the part (b). As the
tudes of motions of the elapsed time steps are shown so that ship moves upward, the buoyancy force magnitude is de-
the effects of waves on the stability of the ship can be re- creased. The resistance of sea water and the gravity force
vealed. In the part (c), the GUI in which users load or unload gradually reduce the heave velocity to zero. At this moment,
cargoes and water into the significant cells is illustrated. In the ship starts to fall because of the gravity force. The result
our implementation, only some significant cells are allowed is displayed in the part (c). As the ship falls, the buoyancy
to store cargoes and water. These cells are shown by gray force magnitude is increased again and becomes greater than
cubes in this image. The stability centers are also displayed the gravity force magnitude. Thus the ship is pushed upward
to demonstrate the influence of the loadings toward the sta- again. After a series of oscillations, the resistance of sea water
bility of the ship. wears down the force of heave. Finally, the ship retains its
balance, as shown in the part (d). The magnitude of the heave
1. Basic Motions: Heave and Pitch motion at each time step is shown in the part (e). As shown by
In the first experiment, a simple ship model is loaded into the graph, the heave magnitude oscillates and gradually con-
our simulator to perform heave and pitch motions. Four im- verges to zero as time elapses.
ages of the heave motion simulation are shown in the parts Another four images are displayed in Fig. 10 to illustrate
(a)-(d) of Fig. 9. In these images, the gravity center G and the pitch motion of the ship. Initially, a force is exerted on
the buoyancy center B are represented by a red ball and a the ship’s front end. The force causes the ship to incline for-
blue ball respectively. The ship body and sea water are ren- ward as shown in the part (a). Since the gravity center and
dered as transparent objects such that the variations of G and the buoyancy center are not collinear, a rotational torque is
B can be monitored. created. As the force is released, the rotational torque makes
Initially, we exert an external force on the ship and press the ship tilt upward, as shown in the part (b). Then the
the ship downward, as shown in the part (a). Since the damping force of sea water cancels the rotational torque and
682 Journal of Marine Science and Technology, Vol. 21, No. 6 (2013 )

Table 1. Geometrical and damping coefficients of the speed

Force
boat.
Length = 40 m Width = 10 m Height = 15 m
Geometrical coef.
Weight = 523 t Draft = 5.05 m
Damping coef. bh = 1.6 bp = 0.6 br = 0.4

(a) Time = 0 (b) Time = 6

wave direction

(a) Time = 0 (b) Time = 2

(c) Time = 10 (d) Time = 90

Pitch Angle (degree)

5
Time (sec)

-5
Ti me
Time = 0

Ti
Ti

m =

m
e=6

e=
10

90

(e)
(c) Time = 4 (d) Time = 6
Fig. 10. (a)-(d) Snapshots of the pitch motion of a ship, (e) the magnitude
of the pitch motion at each time step. Wave height at G (meters)
2

-2
prevents the ship from tilting upward further. However, G and Heave
(meters)
B are not collinear either at this moment. The gravity force 2
creates another rotational torque to incline the ship, as dis-
played in the part (c) of this figure. Finally, after pitching -2
Pitch angle
(degrees)
back and forth for several rounds, the damping effects of sea 3
water diminish the rotational torque and the ship recovers its
stability, as illustrated in the part (d). -3
Roll angle
The variations of the pitch angle are depicted in the part (e). (degrees)
3
Initially, the external force causes the ship to pitch by 5 de-
grees. Then the gravity force and the buoyancy force make the -3
pitch angle fluctuate periodically. As time elapses, the damp- 0 2 4 6 Time (sec.)
ing effects of sea water reduce the rotational torque to zero, (e) Motion magnitudes
and the pitch motion gradually converges to rest. This graph Fig. 11. (a)-(d) Snapshots of the wave-induced motions, (e) the wave
shows a typical behavior of a inclined floating body. It verifies height at G and the magnitudes motion at each time step.
the validity of our pitch model.

2. Compound Motion: Wave Induced Motions the visual effects of the simulation, the damping coefficients
When encountering sea waves, a ship will heave, pitch and of pitch and roll are set to 0.6 and 0.4. Therefore, the motions
roll simultaneously. The resulted motions are called wave- are exaggerated and last longer. The sea waves are composed
induced motions in sea-keeping literatures. In another ex- of three sinusoidal waves of different frequencies, directions
periment, the motions of a speed boat encountering sea waves and amplitudes. The major sea wave approaches the ship from
are simulated by using our simulator. Some results are pre- the north-east direction.
sented and analyzed in this section. The fundamental data of Four snapshots of the simulation are shown in Fig. 11.
the speed boat are listed in Table 1. The ship weighs about The initial scene is displayed in the part (a). At this moment,
523 tons. Its length is about 40 meters. In order to enhance the ship head is in a trough and a crest just passes the rear end.
 S.-K. Ueng: Physical Models of Ship Stability 683

Table 2. Computational cost profiles of two ship motions

(by seconds).
Wave-induced motion sinking Cargoes
Wave generation 0.00838 0.00838
Rendering 0.00258 0.00258
Motion computation 0.00161 0.00180
Frame rate 66 fps 64 fps
(a) Time = 0 (b) Time = 12

The ship heaves slightly downward and inclines forward be-

cause of the passing crest. Two time steps later, another crest
reaches the front end and the ship starts to tilt upward. Since
the water level under the ship body is still low and the buoy-
ancy force magnitude is too small, the ship does not heave
upward. The snapshot is presented in the part (b). As the (c) Time = 27 (d) Time = 80
new crest passing the ship body, the buoyancy force and the
gravity force disturb the stability of the ship. The ship heaves, Heave (meter)
pitches and rolls simultaneously, as shown in the part (c). 6
When this crest passes the rear end, the ship inclines forward.
The ship’s posture is displayed in the part (d). -6
In the part (e), the height of the sea level at the gravity Pitch Angle (degree)
30
center and the magnitudes of heave, pitch and roll at each
time step are displayed. By examining these graphs, we can
conclude that the variations of the sea waves indeed induce -30
Roll Angle (degree)
the motions. As the sea surface varies periodically, so do the 30
Time (sec.)
motion magnitudes. However, there exist latencies between
the variation of wave height and the ship motions, as shown -30
Ti Ti Ti Ti
in the 1st, 2nd and 4th graphs. These latencies are typical m m m
e= e= e=
m
e=
phenomena of wave induced motions and have been men- 0 12 27 80
tioned in the literature of [12]. Our physical models are able
(e) Motion magnitudes
to produce these latencies in the experiment.
The profile of the computational costs of a frame is de- Fig. 12. (a)-(d) Snapshots of the capsizing motions, (e) the wave height at
G and the magnitudes of motion at each time step.
picted in the 2nd column of Table 2. The costs of computing
the sea waves and rendering the scene are contained in the 2nd
and 3rd rows. The costs of computing the ship motions are
listed in the 4th row. The frame rate of this simulation is re- served as the ship to simulate the ship motions caused by
corded in the last row. Our system needs only about 0.00161 unbalanced cargoes. The damping coefficients of the ship are
second to compute the ship motions. The rendering process decreased to 0.1 such that the ship motions are exaggerated.
consumes about 0.00258 second. (The ship model is com- Four images of the simulation are presented in Fig. 12.
posed of 7,000 polygons, and the sea surface contains 32,768 Initially, some heavy cargoes are loaded into the right front
triangles.) The costs for computing the sea waves are about part of the ship. These cargoes are represented by red cubes,
0.00838 second. Since the sea waves are composed of three as shown in the part (a) of this figure. Since the right front
sinusoidal waves and the range of the sea surface is large, part of the ship gains extra masses, the gravity center G (the
computing the sea waves requires more efforts than comput- red ball) is shifted ahead. The ship inclines forward and
ing the ship motions. In total, the frame rate of the simulation rolls right. At time T = 12, the heavy cargoes cause the ship to
is about 66 fps. pitch downward. The ship head is totally immersed into sea
water, and the buoyancy force magnitude is greatly increased,
3. Compound Motion: Cargo Induced Motions as shown in the part (b). The buoyancy force creates a right-
If a ship is loaded with unbalanced cargoes, its gravity center ing lever to bring the ship back from the inclination. How-
may be higher than its meta center. The ship will capsize to ever, the heavy cargoes prevent the ship from recovering its
reach a new stable posture. (In sea-keeping literatures, the initial stability. Instead, the ship can only tilt back by a few
term, capsize, does not mean overturning. Instead, this term is degrees. The scenario is displayed in the part (c). At this
used to describe the process of obtaining a new stability.) moment, the ship head reaches its highest position. Then
In another simulation, a ship of 1,200 tons of displacement is the gravity force causes the ship to incline forward and heave
684 Journal of Marine Science and Technology, Vol. 21, No. 6 (2013 )

cannot support the ship weight and the ship will sink. In the
last simulation, we create holes in the right-front side of the
speed boat and let sea water leak into its body. As sea water
floods into the ship, the ship weight exceeds the maximum
buoyancy force magnitude that the ship body can produce.
Subsequently, the ship sinks into the sea.
Four snapshots of the sinking process are shown in Fig. 13.
At T = 18, a small volume of sea water leaks into the boat. The
(a) Time = 18 (b) Time = 24
boat pitches forward, rolls to the right and heaves downward,
as shown in the part (a). As more water floods into the ship
body, the ship stops rolling but its head is immersed into the
sea. The result is displayed in the part (b). At T = 40, more
cells of the ship body are occupied by sea water and the speed
boat rolls to the left and heaves downward, as illustrated in
the part (c). At T = 95, too much water has flooded into the
ship, and thus most part of the ship sinks into the sea. The
(c) Time = 40 (d) Time = 95 result is shown in the part (d).
In the part (e) of this figure, the motion magnitudes, the
Heave (meter) ship displacement and the buoyancy force magnitude at each
8
time step are depicted. The heave magnitude at each time step
-8 is displayed in the first graph. The heave magnitude is always
Pitch Angle (degree)
30
negative (downward). The ship displacement and the buoy-
ancy force magnitude are represented by the red curve and
-30 the blue curve respectively in the last graph. Before T = 40,
Roll Angle (degree)
30 the displacement and the buoyancy force magnitude increase
simultaneously. At T = 53, the displacement and the buoy-
-30 ancy force magnitude reach their maxima, but the displace-
red : Displacement
blue : Buoyancy (ton)
ment is larger. Therefore, the ship is doomed to sink.
3000
The costs of the sinking simulation are described in the 3rd
Time (sec)
0 Ti Ti Ti Ti column of Table 2. It takes about 0.00180 second to com-
m m m m pute the motions for each time step. As sea water flooding
e= e= e= e=
18 24 40 95 into the ship body, the ship gets extra weights. Our simulator
(e) Motion magnitudes has to re-calculate G in each frame. Therefore the costs for
simulating the sinking process are slightly higher than those
Fig. 13. (a)-(d) Snapshots of the sinking motions, (e) the wave height at G
and the magnitudes of motion at each time step. for simulating the wave-induced motions. The frame rate of
this simulation is about 64 fps.

downward again. In this manner, the ship oscillates for several

VIII. CONCLUSION
rounds. Then the motion magnitudes are gradually reduced
by the damping effects of sea water, and the ship obtains a In this article, efficient methods for computing ship mo-
new stability at T = 80, as shown in the part (d). tions are presented. In our methods, we use the gravity,
The variations of the motion magnitudes are shown in the buoyancy and meta centers to predict the stability of the ship
part (e). In this experiment, there is no sea wave. The only and compute the external forces causing ship motions. Based
external force is created by the heavy cargoes. At T = 12, on our models, we implement a ship motion simulator to
the ship heaves downward by about 6 meters and pitches simulate various ship motions. Experimental results verify
forward and rolls right by about 30 degrees. The ship motion that our methods are physically sound and capable of simu-
magnitudes reach their maxima at this time step. When lating compound ship motions in real time with decent visual
reaching the new stability, the ship’s position is lowered by effects. In traditional ship motion models, sea waves are the
2 meters and inclined forward and rolled right by about 12 sole factor which triggers ship motions. Our models cope with
degrees. not only sea waves but also flooding water and cargoes. Users
are allowed to compose complex ship motions by modifying
4. Sinking Caused by Flooding Water the environment and the ship parameters. For examples, they
If a ship bears severe ship hull damages and too much water can change the damping coefficients, add cargoes into the ship
floods into the ship, the ship weight will exceed the weight of and create holes on the ship hull. Our ship motion simulator
the maximum displaced water body. Thus the buoyancy force can generate the entire sequence of the ship motions auto-
 S.-K. Ueng: Physical Models of Ship Stability 685

matically in real time to reveal the hydrostatic behaviors of the 9. Korvin-Kroukovsky, B. V. and Jacobs, W. R., “Pitching and heaving mo-
ship. Beside computing the ship motions, our simulator also tions of a ship in regular waves,” Transactions of the Society of Naval
Architects and Marine Engineers, pp. 590-632 (1957).
exhibits the variations of key parameters in ship motion 10. Kuhl, J., Evans, D., and Papelis, Y., “The Iowa Driving Simulator: an im-
simulations. By examining the output data, users can learn mersive research environment,” Computer, Vol. 28, No. 7, pp. 35-41 (1995).
more knowledge about sea-keeping and hydrostatic properties 11. Lainiotis, D. G., Charalampous, C., Giannakopoulos, P., and Katsikas, S.,
of ships. “Real time ship motion estimation,” Proceedings of OCEAN’92 Con-
ference, pp. 283-287 (1992).
12. Lewis, E. V., Principles of Naval Architecture, Volume I Stability and
ACKNOWLEDGMENTS Strength, Society of Naval Architects (1988).
13. Lewis, E. V., Principles of Naval Architecture, Volume II Resistance, Pro-
This research is supported by NSC of Taiwan under the pulsion, and Vibration, Society of Naval Architects (1990).
grant numbered as NSC 98-2221-E-019-025. 14. Lewis, E. V., Principles of Naval Architecture, Volume III Motions in
Waves and Controllability, Society of Naval Architects (1990).
15. Peachy, D., “Modeling waves and surf,” Proceedings of SIGGRAPH’86,
REFERENCES pp. 65-74 (1986).
16. Salvesen, N., Tuck, E. O., and Faltinsen, O., “Ship motions and sea loads,”
1. Aryanpour, M. and Ghorashi, M., “Heave and pitch motions of a ship due Transactions of the Society of Naval Architects and Marine Engineers,
to moving masses and forces,” Journal of Sound and Vibration, Vol. 241, pp. 250-287 (1970).
pp. 185-195 (2001). 17. Triantafyllou, M. S., Bodson, M., and Athans, M., “Real time estimation
2. Benson, H., University Physics, John Wiley & Sons, Inc. (1996). of ship motions using Kalman filtering techniques,” IEEE Journal of
3. Carlson, M., Mucha, P. J., and Turk, G., “Rigid fluid: Animating the in- Oceanic Engineering, Vol. 8, No. 1, pp. 9-20 (1983).
terplay between rigid bodies and fluid,” ACM Transactions on Graphics, 18. Ueng, S.-K., Lin, D., and Liu, C.-H., “A ship motion simulation system,”
Vol. 22, No. 3, pp. 377-384 (2004). Virtual Reality, Vol. 12, No. 1, pp. 65-76 (2008).
4. Cieutat, J. M., Gonzato, J. C., and Guitton, P., “A new efficient wave 19. Yuksel, C., House, D. H., and Keyser, J., “Wave particles,” ACM Trans-
model for maritime training simulator,” IEEE Spring Conference on Com- actions on Graphics, Vol. 26, No. 3, Article No. 99 (2007).
puter Graphics, pp. 202-209 (2001). 20. Zhang, X., Jin, Y., Yin, Y., and Li, Z., “Ship simulation using virtual
5. Derrett, D. R., Ship Stability, Butterworth-Heinemann (2000). reality technique,” Proceedings of ACM Siggraph International Confer-
6. Foster, N. and Metaxas, D., “Realistic animation of liquids,” Graphical ence on Virtual Reality Continuum and its Applications in Industry, pp.
Models and Image Processing, Vol. 58, No. 5, pp. 471-483 (1996). 282-285 (2004).
7. Fournier, A. and Reeves, W. T., “A simple model of ocean waves,” Pro- 21. Zhao, X., Xu, R., and Kwan, C., “Ship-motion prediction: algorithms and
ceedings of SIGGRAPH’86, pp. 75-84 (1986). simulation results,” Proceedings of IEEE International Conference on
8. Kim, J., Kim, S., Ko, H., and Terzopoulos, D., “Fast gpu computation of Acoustics, Speech, and Signal Processing, Vol. 5, pp. 125-128 (2004).
the mass properties of a general shape and its application to buoyancy
simulation,” Visual Computer, Vol. 22, No. 9, pp. 856-864 (2006).