ARTICLES

Sci. Rep. Kanazawa Univ.
Vol. 62, pp. 51–68, 2018

An ODE control system of a rigid body on an ocean wave
for a surfer simulation in the SPH method
Reza Rendian S EPTIAWAN∗
Graduate School of Natural Science and Technology, Kanazawa University,
Kakuma, Kanazawa, 920-1192, Japan(Received June 26, 2018 and accepted in revised form August 7, 2018)

Abstract In this work we use a smoothed particle hydrodynamics (SPH) method coupledwith a rigid body simulation to simulate a surfing board on top of an ocean wave. External
forces are applied to the board to represent a surfer trying to control a surfing board. An
ordinary differential equation (ODE) control is used to manipulate the external forces basedon a position, velocity, and an inclination angle of the surfing board. The control system
successfully helps the surfing board to move to and maintain its desired position.
Keywords.
simulation.

1

ODE, rigid body simulation, smoothed particle hydrodynamics, surfing board

Introduction

Surfing on ocean waves poses some interesting mathematical and physics problems. One of themis a modeling of a surfer that controls the movement of the surfing board. In what we call a surfing
problem, the goal is to maintain the position of the surfing board to be on the upslope part ofthe ocean wave as long as possible. The surfer maneuvers the surfing board by controlling theforce given to the board by pushing or shifting their feet on the board, effectively changing the
distribution of forces exerted on the board and modifying its inclination angle.
Since the upslope part of the ocean wave is the domain of interest in the surfing problem, inthis work we choose a frame of reference to be located there. This frame moves together withthe wave. The goal of this work is to control the position of the surfing board to be at the desired
point with respect to the wave. Since the only thing that we can control in an attempt to move thesurfing board toward the desired position is its inclination angle, we propose an ODE control of
the inclination angle that represents a surfer in a surfing problem. The control takes into accountonly the position, velocity, and the inclination angle of the surfing board. The inclination angle
from the ODE control acts as a “target angle” for the surfing board. To reach the target angle wegive two forces at the tips of the board, mimicking the surfer giving forces to the board via theirfeet. To find suitable parameters for which the control system is stable, we perform a stability
analysis of a linearized simplified one-dimensional ODE model of the surfer on an ocean wave.
∗ Email:za22061991@gmail.com

52

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

To verify the capabilities of our ODE control in controlling the dynamics of the rigid bodyunder interactions with the fluid, we simulate the full system using the smoothed particle hydrodynamics (SPH) method. SPH method is one of the most popular Lagrangian solvers of the fluid
equations. It was introduced by Lucy [1], and Gingold and Monaghan [2] to solve astrophysicalproblems. Since then, the SPH method has been used to simulate computational fluid dynamics
problems, such as shockwave problems [3, 4], heat transfer problems [5], multi-phase fluids [6],and also water waves [7, 8].
Another interesting problem which can be solved by the SPH method is a fluid-rigid bodyinteraction. A common method to include a rigid body into an SPH simulation is by discretizing
the rigid body into a set of boundary points which can interact with other SPH fluid points byexchanging momentum using either a normal SPH interaction scheme for hydrodynamics, or amodified interaction scheme to suit the problem. In [9], the interaction between rigid bodies andfluids is done by using a Lennard-Jones potential repulsive boundary force, called the Monaghan
boundary force (MBF), to prevent penetration of fluid points into a rigid body. However, since
MBF is designed to prevent penetration of a point with the maximum velocity, the repulsive force
is thought to be too strong for slow moving points. In [10] an impulse-based boundary force (IBF)was proposed. Recently, a more robust IBF was introduced in [11] by using a sequential impulse
to solve the frictional contact problem with many contact points.
Since the IBF method relies on the normal at the surface of rigid bodies, the calculation fora more complex-shaped rigid body is not easy. Also, IBF is not purely based on hydrodynamic
forces. A more versatile interaction between rigid bodies and fluid points was introduced in [12],which also solves a neighbor deficiency problem near the boundary by the inclusion of boundary
points in the density calculation.
In this work, a rigid body is also discretized into a set of boundary points. An interaction
between fluid points and boundary points is done by using a boundary force based on purelyhydrodynamic forces. To prevent particles to penetrate into the rigid body, we give boundary
points a constant density function which is equal to a reference density, leaving the repulsive force
solely determined by fluid points near the boundary itself.

2

Governing equations

2.1

Fluid dynamics

The motion of the fluid can be observed in two different frames of reference: an Eulerian frame
which is fixed in space so that the flow of fluid is considered as a flux of physical quantities, anda Lagrangian frame which treats the fluid as material points, with each point carrying its physicalquantities as it moves. As an SPH method uses Lagrangian frame intensively, we will discuss thegoverning equation of the fluid in a Lagrangian description. The motion of fluid is governed byfollowing conservations laws which are the Euler’s equations of the fluid dynamics [13]:1. Conservation of mass:Dρ= −ρdiv(u)
Dt

(2.1)

2. Conservation of momentum (for an inviscid fluid):Du1= − ∇p + b
Dt
ρ

(2.2)

Reza Rendian S EPTIAWAN

53

Here ρ, p, and u are density, pressure, and velocity field, respectively, and b is a body force perDunit mass. Dt
is the so-called substantial derivative operator defined as
Df
∂f
=
+u·∇f
Dt
∂t
for any field function f (x,t) and velocity field u.

2.2

Rigid body dynamics

Following [14], let us consider a rigid body with a mass M independent of time specifically given
by
Z

M=

R

ρ(x) dx,

where ρ is a density, and R is a configuration of the rigid body in a principal frame. The center ofmass in the inertial frame is denoted by X(t), and the velocity is U(t) = dtd X(t). The linear motionof the rigid body is described by its linear momentum
G(t) = MU(t).

(2.3)

Let F i (t) = F(xi ,t) be a force that acts on a point xi on the rigid body at a given time t. The total
force acting on the rigid body equals to the rate of change of linear momentum
F(t) =dd(MU(t)) = M U(t) = MA(t),dtdtwhere F(t) = ∑i F i (t) is the resultant force acting on the rigid body, and A(t) is the linear acceleration of the rigid body.
Let J be the principal moment of inertia tensor. J is a diagonal matrix whose components are
given byZ

J1 := J11 =J2 := J22 =J3 := J33 =

ZR
ZR
R

ρ(x)(x22 + x32 ) dx,
ρ(x)(x12 + x32 ) dx,
ρ(x)(x12 + x22 ) dx.

Let R(t) be the rotation matrix of the rigid body. The angular velocity vector ω(t) of the rigid
body is defined in a way that for any fixed vector v ∈ R3 we haved(R(t)v) = W(t)R(t)v = ω(t) × R(t)v,
dt

(2.4)

T
T
where W(t) = dR(t)dt R (t). Notice that R(t) is an orthogonal matrix, that is, R(t)R (t) = I, where
I is the identity matrix. Hence,0= dR(t) TddRT (t)R(t)RT (t) =R (t) + R(t)= W(t) + WT (t).
dt
dt
dt

Therefore W(t) is a skew-symmetric matrix and can be written as0−ω3 (t) ω2 (t)W(t) =  ω3 (t)0−ω1 (t)  .
−ω2 (t) ω1 (t)
0

(2.5)

54

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Figure 1: Sideview of the illustration of the frame of system’s domain.

The vector ω(t) = (ω1 (t), ω2 (t), ω3 (t)) satisfies (2.4), since W(t)w = ω(t) × w for any vectorw ∈ R3 .
The angular momentum of the rigid body is defined as
H(t) = R(t)JRT (t)ω(t).
The total moment of force acting on the rigid body equals to the rate of change of its angularmomentum. Let M i (t) = (xi − X(t)) × F i (t) be the moment of force F i at a point xi on a rigid bodywith a center of mass X(t) at a given time t. Let K(t) = ∑i M i (t) be the total moment of force
acting on a rigid body at a given time t. Now we haveK(t) =d
H(t) = Ṙ(t)JRT (t)ω(t) + R(t)JṘT (t)ω(t) + R(t)JRT (t)ω̇(t),dtwhere Ṙ(t) = dtd R(t) and ω̇(t) = dtd ω(t). This can be rewritten as
˙K̂(t) = ω̂(t) × Jω̂(t) + Jω̂(t),

(2.6)


˙
where K̂(t) = RT (t)K(t), ω̂(t) = RT ω(t), and ω̂(t)= dtd R(t)T ω(t) = RT (t)ω̇(t). (2.6) is theEuler’s equation of rigid body dynamics [14].

2.3

ODE model of the surfer

Let us first introduce the geometry of the domain. In this work we use the frame of reference whichmoves together with the wave and we assume that the wave moves with a constant velocity. The
illustration of the domain can be seen in Figure 1. The goal of our surfing problem is to control
the position of the surfing board to be at the desired position Z̃. In our surfing problem, the only
parameter that we can control is the inclination angle of the board. We introduce the following
ODE control for the inclination angle:θ̇ (t) = a(Z(t) − Z̃) + b(V (t) − Ṽ ) + c(θ (t) − θ̃ ),(2.7)

where Z(t) and V (t) are the third component of position X(t) and linear velocity of the surfing
board, respectively, and θ (t) is the inclination angle of a surfing board. Z̃, Ṽ , and θ̃ are desired
position, desired velocity, and desired angle, respectively. Constants a, b, and c are to be fixed
below. The desired position Z̃ is given. Because of our choice of frame of reference, in a stablecondition, the velocity of the surfing board is zero relative to both the reference frame and theReza Rendian S EPTIAWAN

55

ocean wave, we take Ṽ = 0. θ̃ is a desired inclination angle that helps stabilize the board. Up to
now we do not have any information about θ̃ . Later we run simulations with various values of θ̃to find the value that gives the best result.
To find suitable parameters a, b, and c, we consider a simplified linearized ODE model with
the ODE control (2.7). Then we do a stability analysis of this ODE model and choose parametersa, b, and c that can stabilize the system.
From the definition, we have Ż(t) = V (t). We assume that the acceleration of the rigid bodydepends on an external body force (in this case, gravity) and a drag force. Since we do not know
the drag force in our system, by assuming that the system is close to the stationary point, we canlinearize the acceleration of the system as a function of the position, velocity, and the inclination
angle of the board. Let ξ (t) = (Z(t),V (t), θ (t)) be the unknown of the system. We consider asimplified linearized ODE model in the following form: Ż(t) = V (t),(2.8)
V̇ (t) = −µθ (t) − µvV (t) − µz Z(t) − µ0 ,θ̇ (t) = a(Z(t) − Z̃) + bV (t) + c(θ (t) − θ̃ ),where µ, µv , µz , and µ0 are constants related to the drag and gravity forces. Now, let us assumeµv = µz = 0 for simplicity. Later in the results section we will see that this assumption is notjustified, but it does not seem to influence the stability.
We can write (2.8) in a matrix form as




0
0 1 0

ξ˙ (t) =  0 0 −µ  ξ (t) + 
−µ0
a b c
−aZ̃ − cθ̃=: Aξ (t) + ζ .
The stationary point is stable if all eigenvalues of matrix A have negative real parts.
The characteristic equation of A isdet(A − λ I) = −λ 3 + cλ 2 − bµλ − aµ = 0.

(2.9)

By using Vieta’s formulas, we see that for all roots of (2.9) to be negative, the value of a and b
must be positive, and c must be negative.
Let us set the roots of (2.9) to be λ1 = −3, λ2 = −4, and λ3 = −5. This yields the equation−λ 3 − 12λ 2 − 47λ − 60 = 0.

(2.10)

47
Comparing this with (2.9), we get a = 60
µ , b = µ , and c = −12.
The inclination angle from the ODE control is passed to the simulation as a target inclination
angle for controlling the movement of the rigid body. We model the controlling action of the surfer
on a surfing board by using a two-points model, imitating the locations where the surfer places theirfeet. The location of those points in our model, each is denoted as rc1 and rc2 , respectively, canbe seen in Figure 2. They are located at the tips of the elongated axis of the rigid body. The total
force is constant, but the individual forces are controlled by a ratio T (t) as

Fc1 (t) = T (t)W

(2.11)

Fc2 (t) = (1 − T (t))W,(2.12)56

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Figure 2: The locations of contact points.

where W is the weight of the surfer, Fc1 (t) and Fc2 (t) are the forces given at each point rc1 and rc2 ,respectively, at a given time t, with the total force being constant as they represent the weight ofthe surfer. The ratio is controlled by the linear function
T (t) = 0.5 − σ (θ̂ (t) − min(θ (t), θm )),(2.13)

where θ (t) is the angle from the ODE control in Section 2.3, θm is an allowed maximum inclination angle, and θ̂ (t) is the observed angle at a given time t. σ is a given constant. Here weuse a clipped function for the inclination angle since the surfer does not want the front part oftheir surfing board to be immersed into the water. Note that we choose the inclination angle to be
negative when the board is inclined upward, and vice versa. We do not impose T (t) ∈ [0, 1], whichis reasonable if the surfer is using straps on their feet.

3

SPH method

3.1

Basic idea of SPH method

The basic idea of SPH comes from a convolution between any field function f (x) with a sufficiently smooth mollifier ψ(x),( f ∗ ψ) (x) :=Z

f (y)ψ(x − y) dy.

Rn

If f ∈ L1 (Rn ) and ψ ∈ Cck (Rn ), then ( f ∗ ψ) ∈ Cck (Rn ), which means the convolution product is assmooth as the smoothest of the convolved functions, see [15] for more details.
Let us define
1 xψh (x) := n ψ
.
h
h
R

R

Clearly Rn ψh (x) dx = Rn ψ(x) dx, but the “mass” of ψh is more concentrated to the origin as
R
h → 0. If Rn ψ(x) dx = 1, the family {ψh }h>0 is called an approximate identity. If f ∈ Ck (Rn ) forsome 1 ≤ k < ∞, then we have f ∗ ψh → f uniformly as h → 0,
f (x) ≈

ZRnf (y)ψh (x − y) dy,

(3.1)

57

Reza Rendian S EPTIAWAN

and by following the differentiation of a convolution in [15], the approximation for the derivativeof f is∂ α f (x) ≈ZRnf (y)∂ α ψh (x − y) dy,|α| ≤ k,(3.2)

for ψ ∈ Cck (Rn ). Both (3.1) and (3.2) are SPH approximations for a field function and its derivativein an integral form.
We can choose a smoothing function ψ. In this work we are using a compactly-supported
piecewise cubic kernel function [16]
(2 − |x|)3 − 4(1 − |x|)3 , 0 ≤ |x| < 1αn 
ψ(x) =(3.3)(2 − |x|)3 , 1 ≤ |x| < 2
6 0, 2 ≤ |x|
15
3
where αn is 1, 7π, or 2π
for n = 1, 2, 3 respectively. Note that ψ ∈ Cc2 (Rn ).

3.2

Discretization of the fluid equations using SPH

To discretize the fluid equations, we keep track of the physical quantities at discrete points thatmove with the fluid. Let us assume that we have N points with positions ri = ri (t) which arereasonably uniformly distributed. We approximate (3.1) and (3.2) by a discrete integration,NZ

f (x) ≈

Rnf (y)ψh (x − y) dy ≈ ∑ f (ri )ψh (x − ri )V (Ei ),

(3.4)

i=1

and
∂ α f (x) ≈

Z
Rn

N

f (y)∂ α ψh (x − y) dy ≈ ∑ f (ri )∂ α ψh (x − ri )V (Ei ),

(3.5)

i=1

where V (Ei ) is a volume of the set Ei around a point ri and we assume that the integrand is aconstant or close to being linear if Ei is symmetric and ri is at its center of mass. A common way
to approximate V (Ei ) is to use a mass and density at point ri ,Vi (Ei ) ≈

mi
,
ρi

(3.6)

where mi and ρi are mass and density at a point ri , respectively. By using (3.4) for the density field
function and by (3.6), we haveNNρ(x) ≈ ∑ ρi ψh (x − ri )V (Ei ) ≈ ∑ mi ψh (x − ri ).
i=1

i=1

Then we have an approximation for V (Ei ) asV (Ei ) ≈mi.
N
∑ j=1 m j ψh (ri − r j )(3.7)

Another way to update the density is to use an SPH discretized approximation of (2.1) at apoint ri at a given time t, that is,Nmjdρi= −ρi ∑
u j · ∇ψh (ri − r j ),dt
j=1 ρ j

(3.8)

58

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

where mi is a mass of a point ri , ρi (t) and ui (t) are density and velocity of a point ri at a given timet, respectively. The more common form for (3.8) is its anti-symmetrized version of (3.8) from [16]Nmjdρi(ui − u j ) · ∇ψh (ri − r j ).
= ρi ∑dt
j=1 ρ jIn this work we use (3.9) to update the density of points.
The SPH discretized approximation for (2.2) at a point ri at a given time t is


pj
dui
1 N
= − ∑ m j ∇ψh (ri − r j ) + bi ,dt
ρi j=1ρj(3.9)(3.10)

where pi (t) is a pressure at a point ri at a given time t, and bi is a body force per unit mass at a pointri . Remember that in this work we assume the fluid to be an inviscid fluid. The anti-symmetrizedversion [16] of (3.10) is
!
!
N
pj
pi
dui
= − ∑ mj+∇ψh (ri − r j ) + bi .
(3.11)dt
ρi2 ρ 2jj=1
In this work we use (3.11) to update the acceleration of points.
The pressure at a point ri is assumed to be a function of density by using Tait’s relation as
follows [17, 18, 19]:
 γ

ρi
c2 ρ0
−1 ,
(3.12)pi =γ
ρ0where ρ0 is a reference density of fluid, c is a speed of sound in a fluid, and γ = 7 for water-likefluid.
As is common in the SPH literature, we use the leapfrog time integrator scheme to evolve
physical quantities of material points as follows τ τ  duiui t += ui t −+(t)τ,(3.13)
22  dtτri (t + τ) = ri (t) + ui t +τ,(3.14)2 τ  duτi+ui (t + τ) = ui t +(t) ,(3.15)
2
dt
2
τ
dui
τ
ui (− ) = ui (0) −(0) ,(3.16)
2dt
2where τ is a timestep.

3.3

Rigid body dynamics

Since J is a diagonal matrix, we can write (2.6) componentwise asJα ω̂˙ α (t) = K̂(t)α + ω̂β (t)ω̂γ (t) Jβ − Jγ ,(3.17)

where (α, β , γ) = (1, 2, 3) , (2, 3, 1) and (3, 1, 2). Following [20], in the spirit of the leapfrog time
integrator, the angular velocity vector ω̂(t) can be updated using the fixed point iteration τ τ 
 τ

(m+1)
(m)
(m)
ω̂α
t+
= ω̂α t −+K̂α (t) + ω̂β (t)ω̂γ (t) Jβ − Jγ ,(3.18)22Jα59

Reza Rendian S EPTIAWAN

where τ is a timestep, and (m) is an iteration step. Initially, we can set
 τ
 τ
(0)
ω̂α t += ω̂α t −.
22
We also choose
(m)

(m)

ω̂β (t)ω̂γ (t) = τ τ 
1  τ  τ(m)(m)
ω̂β t −
ω̂γ t −
+ ω̂β
ω̂γ
.
t+
t+
2
2
2
2
2

Later, we update the rotation matrix by using (2.4) and in the sense of leapfrog integrator asd  τ
R t+
dt 2τ  τ= R(t) + τW t +
R t+
,
2
2

. After some algebraic operations,where W is a matrix from (2.5) and R t + τ2 = (R(t)+R(t+τ))2
we haveτ  τ −1 τ  τ R(t + τ) = I − W t +I+ W t +R(t)222 2 τ
=Θ t+R(t),(3.19)
2
R(t + τ) = R(t) + τ

where I is the identity matrix, and τ τ  τ −1 τ  τ Θ t+:= I − W t +I+ W t +
,
2
2
2
2
2
or more explicitly




2
2
 τ  I 1 − τ4 ω 2 t + τ2 − τW t + τ2 + τ2 (ω ⊗ ω) t + τ2
Θ t+
=
,

2
2
1 + τ ω2 t + τ
4where ω 2 t + 2


τ

3.4

2


= (ω · ω) t + τ2 .

Rigid body discretization and coupling with SPH

The rigid body is discretized into Nb points with a regular distance between them. In this work we
set the distance between rigid body points to be equal to h. A position of each point is denoted byri (t) = (ri,1 (t), ri,2 (t), ri,3 (t)) at a given time t. To make the rigid body dynamics calculation easier,we configure the calculation of the rigid body to be in a principal frame. Let J̌ be the moment ofinertia tensor of the rigid body. By assuming the density of the rigid body is uniform and equals
to ρb , components of J̌ areNb
2
2
Jˇ11 = ρb ∑ ri,2+ ri,3dx,

Jˇ12 = Jˇ21 = −ρb ∑ ri,1 ri,2 dx,22Jˇ22 = ρb ∑ ri,1+ ri,3dx,

Jˇ13 = Jˇ31 = −ρb ∑ ri,1 ri,3 dx,

22Jˇ33 = ρb ∑ ri,1+ ri,2dx,

Jˇ23 = Jˇ32 = −ρb ∑ ri,2 ri,3 dx.

i=1
Nb

i=1
Nb

i=1

Nb
i=1
Nb

i=1
Nb
i=1

(3.20)

60

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

We calculate the principal moment of inertia tensor J of the rigid body by doing an eigendecomposition of the moment of inertia tensor J̌. Eigenvalues of J̌ are diagonal components of J, while
their corresponding eigenvectors construct an orthogonal rotation matrix R(0).
J̌ = ΛJΛ−1 ,R(0) = Λ.

(3.21)

Similarly to fluid points, we keep track of the physical quantities associated with rigid body
points as well. Rigid body points also have a density function which we set to be equal to areference density ρ0 . Mass and volume of rigid body points are set to match our SPH configuration
described in (3.6) and (3.7). Note that the density of rigid body points is a quantity completely
different from the density of the rigid body ρb .
After the discretization process, interactions between fluid points and rigid body points arehandled in the exact same way as with interactions between fluid points. Fluid points see the rigid
body points as exactly the same objects as other fluid points and do not discriminate interactions
between them, on both the density update step and the acceleration update step.
As mentioned before in the introduction, in this work we consider a purely hydrodynamics
force for interactions between rigid body points and fluid points. Based on (3.11), since the densityof a rigid body point is equal to the density reference yielding zero pressure for a rigid body point,
the total force applied to rigid body point ri is!
!
Npj
fi = −mi cV(3.22)∑ m j ρ 2 ∇ψh (ri − r j ) + bi ,
j
j=1
3

h
. The constant cV is needed since there is a discrepancy on an SPH volume
where cV = V (Ei)
calculation V (Ei ) compared to the volume of the cube with length of the edge h occupied by eachrigid body point. Note that the sum in (3.22) ranges only over fluid points. That is, we do notconsider interactions between rigid body points.
To calculate a linear movement of the rigid body, from (2.3), we have!
!
Nb
1
A(t) =(3.23)∑ fi (t) + Fc1 + Fc2 ,M
i=1

where M is the mass of the rigid body, A(t) is a linear acceleration of the rigid body, fi (t) is thetotal force applied to a rigid body point ri , Fc1 and Fc2 are forces from the ODE control given at acontact point 1 and contact point 2, respectively. Updating a linear velocity and a position of thecenter of mass of the rigid body can be done using the leapfrog time integrator scheme in a similar
fashion described in (3.13)–(3.16).
Next, to calculate a rotation movement of the rigid body, we need to sum total moments offorce acting on the rigid body,!
NbK(t) =∑ (ri (t) − X(t)) × fi (t) + (rc1 (t) − X(t)) × Fc1 + (rc2 (t) − X(t)) × Fc2 ,(3.24)

i=1

where K(t) and X(t) are the total moment of force acting on the rigid body and the position of
the rigid body at a given time t, respectively, rc1 and rc2 are a position of a contact point 1 andcontact point 2, respectively. We use Euler’s equation of rigid body dynamics (2.6) and the iterative
scheme (3.18) to update the angular velocity of the rigid body. Then we update the rotation matrixof the rigid body by using (3.19). The position of the rigid body point is updated based on its
relative position to the center of mass of the rigid body asri (t + τ) = X(t + τ) + R(t + τ)RT (t) (ri (t) − X(t)) ,(3.25)

Reza Rendian S EPTIAWAN

61

where τ is the timestep. The velocity of a rigid body point is updated based on a linear velocity
and an angular velocity of the rigid body τ τ
ui (t + τ) = U t ++ A(t) + ω(t + τ) × (ri (t + τ) − X(t + τ)) ,(3.26)
22
where ui (t) is a velocity of a rigid body point ri at a given time t, U(t) and ω(t) are respectively alinear velocity and an angular velocity of the rigid body at a given time t.

3.5

Algorithm of the simulation

The simulation is started by an initialization of all parameters needed for the simulation and initialconfiguration for all SPH points as discretized representations of fluid and rigid body. The main
loop of the simulation is started with updating the density, pressure, and hydrodynamics forcesbased acceleration of all SPH points. The forces acting on rigid body points are used to updatea linear movement of the rigid body. Then the ODE control updates the contact forces needed
to control the rigid body. Then all forces including the contact forces are used to calculate therotational movement of the rigid body. Next, we update the velocity and position of the rigid body
points based on the evolution of the rigid body. Finally, we update the velocity and position of thefluid points, and go back to the starting point of the main loop. The full algorithm of the simulation
can be seen in Figure 3.

4

Simulation results and discussions

In this work, we use the frame of reference which moves together with the wave, assuming thewave moves with constant velocity V f = 2.5. The inclination angle of the domain (φ in Figure 1)π
is set to be φ = 18
, causing the gravity in our simulation frame to be slanted to −z-axis. Gravity
acts as an external body force per unit mass is set to be (0, −9.81 cosφ , −9.81 sinφ ). The size ofthe domain is set to be 1 × 0.6 × 1.6 in x, y, z axis, respectively. The system has a periodic boundary
condition in x-axis, non-zero Dirichlet boundary on some parts of the left boundary by using ghost
points (rendered with light blue color in Figure 4(a)) for −0.3 ≤ y < −0.12. The bottom boundaryis also set to be a non-zero Dirichlet boundary by using a non-moving boundary points (renderedwith dark blue color in Figure 4(a)) for −0.8 ≤ z < 0.56, and free boundary on the right boundary,
top boundary, the left boundary for −0.12 ≤ y < 0.3, and the bottom boundary for 0.56 ≤ z ≤ 0.8.
Note that the origin is located at the center of the domain. The depth of the fluid is 0.18. The sizeof the rigid body is 0.2 × 0.06 × 0.8. The rigid body is represented by red points in Figure 4(a).
Initially, the center of mass of rigid body is positioned at (0, −0.04, −0.27).
The density reference ρ0 is set to be ρ0 = 1000, while the density of the rigid body is ρb = 100.
The fluid is initialized to have an initial velocity ui (0) = V f = 2.5 toward +z-axis, and initial
density to be ρi (0) = ρ0 for all fluid points ri . We use the piecewise cubic kernel (3.3) as themollifier function and initialize the points on a regular grid with a distance h along each axis,
V (Ei ) and mi for all points are V (Ei ) = 8 × 10−6 and mi = 0.008. We set the parameter of kernel
function h = 0.02. Time step size is set to be τ = 0.0005 with speed of sound chosen to be c = 20.
W from (2.11) and (2.12) is set to W = 10. We choose σ and θm of (2.13) to be σ = 10 andθm = −0.05. The positions of contact point 1 and contact point 2 are rc1 = (0., −0.03, −0.4) and
rc2 = (0., −0.03, 0.4), respectively, relative to the position of center of mass of the rigid body.
Reynold’s number for our configuration, using the real value of the dynamic viscosity of the62

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Start simulation.
Initialize parameters and points configuration, including set t = 0, τ, and tmax .
Calculate principal inertia tensor J androtation matrix R(0) by (3.20)–(3.21).

t ≤ tmax ?
yesnoUpdate density ρi by (3.9).

Finish.

Update pressure pi by (3.12).
i
Update acceleration dudt and force f i by (3.11) and (3.22).

Update linear acceleration A of rigid body by(3.23), update linear velocity V and position X of
rigid body by a scheme similar with (3.13)–(3.16).
Back tomain loop.

Solve ODE (2.8), update contact
force Fc1 and Fc2 by (2.11)–(2.13).
Calculate total moment of forces K by (3.24), update angularvelocity ω by (3.18), and rotation matrix R by (3.19).
Update rigid body points’ position riby (3.25) and velocity ui by (3.26).

t = t +τUpdate velocity ui and position riof fluid points by (3.13)–(3.16).

Figure 3: Algorithm of the simulation.

water from [21], isRe =ρ0V f L 1000 × 2.5 × 0.2== 5 × 105 ,ν
0.001

where ν is the real dynamic viscosity of the water at 20 degree Celsius, and we take the characteristic length L to be the width of the rigid body. It is clear that the viscosity force is negligible
compared to the inertial force of the fluid. Hence, in this work we can neglect the viscosity andassume the fluid to be inviscid.
The free boundary condition is implemented by changing the type of any fluid points leaving
the domain into a ghost point whose velocity does not change with time and whose density always
equals to the reference density ρ0 , but still interacts with other points. If a ghost point reenters thedomain, it is marked as a normal fluid point again. But if it leaves the domain farther than h, thepoint is removed from the simulation.
Before we run the actual simulation, we run the “relaxation” process to stabilize the flow of
the water up to t = 1.5. The initial condition after relaxation can be seen on Figure 4(b).

63Reza Rendian S EPTIAWAN

(a)

(b)

Figure 4: (a) The initial configuration of the system, and (b) the initial condition after relaxation process.

(a)

(b)

(c)

(d)

Figure 5: Snapshots of the surfing board simulation without an ODE controller at (a) t = 1.00 s, (b) t =2.00 s, (c) t = 2.50 s, and (d) t = 2.75 s.

First, we try to run a simulation without an ODE control. We can see in Figure 5 that the board
cannot maintain its position and drifts away with the flow of the fluid.
To find the best θ̃ and µ, for each Z̃ ∈ {−0.6, −0.5, −0.4, −0.3, −0.2} we try each combination of θ̃ ∈ {−0.05, −0.06, −0.07, −0.08, −0.09, −0.1} and µ ∈ {1, 2, 5, 10, 20, 50}. For eachsimulation, we take the data of the position of the center of mass of the rigid body for each time
step to assess the quality of the corresponding parameters configuration. Let Z(µ, θ̃ , Z̃, mτ) be
the third component of the position of the center of mass of the rigid body at a time step mτ fora simulation with a parameters configuration µ, θ̃ , and Z̃, with m ∈ {0, 1, . . . , ϒ}. To assess thechosen parameters configuration, we take an average of the difference between Z(µ, θ̃ , Z̃, mτ) and64

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Figure 6: The graphs of the average positional error for different parameters θ̃ , µ, and Z̃.

a corresponding Z̃ for the whole simulation time,εZ (µ, θ̃ , Z̃) :=∑ϒm=0 Z(µ, θ̃ , Z̃, mτ) − Z̃
.
ϒ+1

(4.1)

Let us call ε¯Z the average positional error. The graphs of ε¯Z can be seen on Figure 6. Then we takean average of the average positional error for each Z̃,εZ (µ, θ̃ ) :=∑Z̃∈Z̃set εZ (µ, θ̃ , Z̃),#Z̃setZ̃set = {−0.6, −0.5, −0.4, −0.3, −0.2} ,(4.2)

and find cumulative errors for each θ̃ and µ. The average of average positional error and its
cumulative errors for each θ̃ and µ can be seen in Table 1.
From Table 1 we can see that θ̃ = −0.07 and µ = 5 give the smallest cumulative errors for θ̃and µ, respectively. Now let us see more in detail the simulation results for θ̃ = −0.07 in Figure 7.
As we can see in Figure 7, smaller values of µ give more oscillations to the position compared
to larger µ values. But as the larger value of µ dampens the amplitude and frequency of oscillation
to the position, it also shifts the stable position of surfing board. That problem also occurs when
we solve the ODE (2.8) directly by using an ODE solver which we can see in Figure 8.
The direct ODE solver is used to observe the behavior of the solution under different parameters
choices. The parameters a and b of the ODE control in (2.8) depend on µ. Note again that we
do not know the actual value µ of our system. Hence, we carry out trials and errors by varying
the value of µ in our simulation, yielding a and b that differ from the correct ones. By using the

65

Reza Rendian S EPTIAWAN

Table 1: Table of the average of average positional error εZ and its cumulative errors for each θ̃ and µ.

µ\θ̃

−0.05

−0.06

−0.07

−0.08

−0.09

−0.10

C.E.

1
2
5
10
20
50

5.74e-02
5.30e-02
5.02e-02
5.97e-02
8.47e-02
1.28e-01

5.81e-02
5.13e-02
5.29e-02
5.13e-02
7.04e-02
1.08e-01

5.33e-02
5.08e-02
4.62e-02
4.66e-02
5.07e-02
7.15e-02

5.80e-02
4.67e-02
4.76e-02
4.65e-02
5.07e-02
7.00e-02

5.98e-02
5.48e-02
4.53e-02
5.08e-02
6.92e-02
9.97e-02

5.44e-02
5.94e-02
4.89e-02
6.29e-02
9.33e-02
1.72e-01

3.41e-01
3.16e-01
2.91e-01
3.18e-01
4.19e-01
6.49e-01

C.E.

4.33e-01

3.92e-01

3.19e-01

3.20e-01

3.80e-01

4.91e-01

Figure 7: The z-axis-component position of surfing board for θ̃ = −0.07.

direct ODE solver we want to see the effect of the choice of µ in simulation that does not matchthe actual µ of the system. In the direct ODE solver we choose µactual = 20 and use various µ,effectively varying the values of a and b. We choose Z̃ = −0.5. From Figure 8 we can see thatas the value of µ increases, the oscillation decreases in both amplitude and frequency. But theincrease of µ also shifts the stable position farther from the desired position Z̃ = −0.5, justifying
the same behavior that occurs in our simulation. Adding an additional ODE control for θ̃ might
be the solution for this problem. θ̃ acts as a “target” inclination angle in our current ODE control.
By controlling θ̃ , it is possible to disturb the stability of the system when the board is not locatedat the desired position, forcing the board to nudge slowly to the desired position.
Notice that in the direct ODE solver, the value of µ affects the oscillation of θ and V instead

66

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Figure 8: Solution of the ODE system for different µ.

of θ and Z in our SPH simulations. This happens because our SPH simulation cannot translate thechange of the inclination angle into the change of the velocity fast enough. The delayed response
in velocity propagates the oscillation in the inclination angle to the position in SPH simulations.
By tweaking σ from (2.13) it is possible to transfer the change of the inclination angle θ from the
ODE control to the inclination angle of the board in the SPH simulation θ̂ faster, which makes afaster change in the velocity and helps the surfing board stabilizes faster.
The oscillation in Figure 7 also occurs since we set µv and µz from (2.8) to be zero, ignoring
the dependency of acceleration to the position and velocity. In fact, this is not correct since thedrag force depends not only on the shape of the interface between the fluid and the rigid body,
but also on the relative velocity between both of them. Initially we assume that the flow has aconstant velocity through the whole domain on z-axis of our frame. But, in reality, the velocity of
the flow depends on the position, as the gravity slows the velocity of the flow as a consequence ofour choice of the frame.
Snapshots of a simulation with the ODE control can be seen on Figure 9.

5

Summary

An ODE control was successfully implemented into a coupled inviscid fluid-rigid body SPH simulation in an attempt to control the movement of the surfing board. For our system, the best valuesfor θ̃ and µ are θ̃ = −0.07 and µ = 5. Although µ = 5 gave the best result, it does not mean theactual µ of our system is equal to 5. As shown in the result of the direct ODE solver, as the µ67

Reza Rendian S EPTIAWAN

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: Snapshots of the surfing board simulation with an ODE controller with parameters Z̃ = −0.2,µ = 5, and θ̃ = −0.07 at (a) t = 0.20 s, (b) t = 0.90 s, (c) t = 2.50 s, (d) t = 5.00 s, (e) t = 8.00 s, and (f)t = 10.00 s.

increases, the stable position shifts farther from the desired position. This problem can be solved
by adding an additional ODE control for θ̃ that can nudge the board toward the desired position.
By tweaking values of µz , µv and σ , it is possible to make the system control the surfing board
more accurately and swiftly. In this work we neglected the effect of viscosity since the inertial
force is dominating the viscous force.

Acknowledgement
This work is a part of doctoral research of the author. The author would like to thank his Ph.D.
supervisor Professor Seiro Omata and his colleague Professor Norbert Pozar for their invaluable
supports, advices, and moral supports. This work cannot be done without helps from both ofthem. The author would like to thank the Ministry of Education, Culture, Sports, Science, and68

An ODE control system of a rigid body on an ocean wave for a surfer simulation in the SPH method

Technology (MEXT) for their financial support.

References[1] Lucy, L. B.: A numerical approach to the testing of the fission hypothesis, Astro. J. 82 no. 12 (1977), 1013–1024.
[2] Gingold, R. A. and Monaghan, J. J.: Smoothed particle hydrodynamics: theory and application to non-sphericalstars, Mon. Not. Roy. Astron. Soc. 181 (1977), 375–389.
[3] Monaghan J. J.: Shock simulation by the particle method SPH, J. Compu. Phys. 52 (1983), 374–389.
[4] Monaghan, J. J.: SPH and Riemann solvers, J. Compu. Phys. 136 (1997), 298–307.
[5] Clear, P. W. and Monaghan, J. J.: Conduction modelling using smoothed particle hydrodynamics, J. Compu.
Phys. 148 (1999), 227–264.
[6] Colagrossi, A. and Landrini, M.: Numerical simulation of interfacial flows by smoothed particle hydrodynamics,
J. Compu. Phys. 191 (2003), 448–475.
[7] Lo, E. Y. M. and Shao S.: Simulation of near-shore solitary wave mechanics by an incompressible SPH method,
Appl. Ocean Res. 24 (2002), 275–286.
[8] Dalrymple, R. A. and Rogers, B. D.: Numerical modelling of water waves with SPH method, Coast. Eng. 53(2006), 141–147.
[9] Monaghan, J. J., Kos, A. and Issa, N.: Fluid motion generated by impact, J. Waterw. Port Coast. 129 (2003),
250–259.
[10] Oh, S. , Kim, Y. and Roh, B.: Impulse-based rigid body interaction in SPH, Comp. Anim. Virtual Worlds 20(2009), 215–224.
[11] Li, Y. and Asai, M.: Fluid-rigid body interaction simulation based on a stabilized ISPH method incorporated
with the impulse-based rigid body dynamics, Trans. Jap. Soc. Comput. Eng. Sci. 2018 Issue 2 (2018), Article No.
20182010.
[12] Akinci, N., Ihmsen, M., Akinci, G., Solenthaler, B. and Teschner, M.: Versatile rigid-fluid coupling for incompressible SPH, ACM Trans. Graph. 31 Issue 4 (2012), Article No. 62.
[13] Wendt, John F.: Computational Fluid Dynamics: An Introduction, Springer, 2009.
[14] Rao, Anil V.: Dynamics of Particles and Rigid Bodies: A Systematic Approach, Cambridge, 2006.
[15] Folland, G. B.: Real analysis, 2nd ed., Pure and Applied Mathematics, John Wiley & Sons, Inc., New York, 1999.
[16] Monaghan, J. J. and Lattanzio, J. C.: A refined particle method for astrophysical problem, Astron. Astrophys.
149 (1985), 135–143.
[17] Batchelor, G. K.: An Introduction to Fluid Dynamics, Cambridge, 1967.
[18] Monaghan, J. J.: Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994), 399–406.
[19] Monaghan, J. J.: Smoothed particle hydrodynamics, Rep. Prog. Phys. 68 (2005), 1703–1759.
[20] Omelyan, Igor P. : Algorithm for numerical integration of the rigid-body equations of motion, Phys. Rev. E 58
no.1 (1998), 1169–1172.
[21] Kestin, J., Sokolov, M. and Wakeham, W. A. : Viscosity of liquid water in the range −8o C to 150o C, J. Phys.
Chem. Ref. Data 7 no.3 (1978), 941–948.