Thesis for the Degree of Master of Science in Engineering Physics I N T E R A C T I V E S I M U L AT I O N O F H Y D R O D Y N A M I C S F O R A R B I T R A R I LY S H A P E D O B J E C T S Author: s a n d r a å l s t i g s a n d b e rg Supervisor: n i k l a s m e l i n Co-supervisor: c l au d e l ac o u r s i è r e Examiner: m a rt i n s e rv i n June 12, 2014
Sandra Ålstig Sandberg: Interactive simulation of hydrodynamics for arbitrarily shaped ob-jects, c June 12, 2014
A B S T R A C T
i n t e r ac t i v e s i m u l at i o n o f h y d ro dy n a m i c s f o r a r b i t r a r i ly s h a p e d o b j e c t s
Simulating hydrodynamics can require extensive calculations which becomes a problem when doing interactive simulations. This thesis investigates an efficient method for hydro-dynamic simulations with the effects of buoyancy, drag, lift and added mass that is imple-mented and tested with the help of AgX Dynamics using triangulated meshes.
For buoyancy, drag and lift a method of numerical integration over triangles was used to calculate the forces and torques acting on each triangle of a mesh. For added mass a large part of the calculations could be done before the simulation starts using a Boundary El-ement Method (BEM). The final value for the added mass was calculated each time step based on how the object was submerged.
The method of triangle integration produced results that were close to the analytical values with a certain mesh dependence. The results had an increasing accuracy when the mesh had a more exact representation of the object. The drag and lift coefficients could however be better adjusted. The added mass results also had a mesh dependence, but with accuracy increasing with number of triangles even for shapes that already had an exact representa-tion, e.g. a cube. For a fully submerged sphere with 4900 triangles the maximum error for the added mass was 0.65%.
The time required for precalculations using BEM had a rapid growth with increasing number of triangles due to the factorization of a dense matrix that has a complexity of O(n3). For the hydrodynamic calculations done each time step the time requirement increased linearly with number of triangles.
S A M M A N FAT T N I N G
i n t e r a k t i v s i m u l e r i n g av h y d ro dy n a m i k f ö r o b j e k t av g o d t y c k l i g f o r m
Att simulera hydrodynamik kan kräva omfattande uträkningar vilket blir ett problem när man ska göra interaktiva simuleringar. Det här examensarbetet undersöker en effektiv metod för att göra hydrodynamiska simuleringar som tar hänsyn till flytförmåga, motståndskrafter, lyftkraft och adderad massa som är implementerad och testad med hjälp av AgX Dynamics och triangulerade objekt.
För att beräkna flytförmåga, motståndskrafter och lyftkraft så användes en metod för nu-merisk integrering av trianglar som beräknar krafter och vridmoment på alla trianglar. För den adderade massan så gjordes en stor del av uträkningar i förväg med hjälp av en rand-värdesmetod (BEM). Sedan kunde den slutliga adderade massan beräknas i varje tidssteg beroende på vilka delar av objektet som låg under vattenytan.
Triangelintegreringen gav goda resultat som låg nära de analytiska värdena med ett visst trianguleringsberoende. Precisionen ökade när trianguleringen hade en mer exakt represen-tation av objektets form. Koefficienterna för motståndskrafterna och lyftkraften skulle dock kunna justeras bättre. Även resultaten för added mass hade ett trianguleringsberoende men där så ökade noggrannheten med antalet trianglar även för objekt som hade en exakt repre-sentation, till exempel en kub. För en sfär med 4900 trianglar som var helt under vattenytan så blev det största felet för added mass 0.65%.
Tiden som behövdes för att använda BEM ökade kraftigt med antal trianglar på grund av faktoriseringen som görs på en tät matris som har komplexitet O(n3). För de uträkningar som gjordes varje tidsteg så ökade tidsåtgången linjärt med antal trianglar.
A C K N O W L E D G E M E N T
This project was done at Algoryx Simulation AB as a thesis for the degree of master of science in engineering physics. A lot of people helped me along the way and I would like to thank them all greatly.
To my supervisor Niklas Melin for all discussions on hydrodynamic effects and all the help and support during this project.
To Claude Lacoursière for having the patience to answer all my questions regarding the theory used in this project.
To everyone at Algoryx for helping me with everything from debugging to tests and for making me feel welcome.
C O N T E N T S
1 i n t ro d u c t i o n 1
1.1 Background . . . 1
1.2 Related work . . . 1
1.3 Scope and approach . . . 2
1.3.1 Method . . . 2
1.3.2 Trimesh . . . 2
1.3.3 Water . . . 2
1.4 Outline . . . 3
2 t h e o ry 5 2.1 Numerical integration over triangles . . . 5
2.2 Boundary Element Method (BEM) for Laplace’s equation . . . 6
2.2.1 Numerical integration . . . 7
2.3 Hydrodynamic effects . . . 9
2.3.1 Buoyancy . . . 9
2.3.2 Drag & Lift . . . 10
2.3.3 Added Mass . . . 12
3 s o f t wa r e i m p l e m e n tat i o n 15 3.1 Algorithm description . . . 15
3.2 Pressure field . . . 16
3.3 Constants and coefficients . . . 16
4 r e s u lt s 19 4.1 Buoyancy . . . 19
4.2 Pressure drag forces . . . 21
4.3 Viscous drag forces . . . 22
4.4 Lift forces . . . 23 4.5 Added mass . . . 24 4.6 Performance . . . 26 4.7 Demonstration . . . 27 5 d i s c u s s i o n 29 6 c o n c l u s i o n s 33 6.1 Further work . . . 33 r e f e r e n c e s 35
1
I N T R O D U C T I O N 1.1 b ac kg ro u n dSimulations are a way of observing the physical behavior of a system under certain condi-tions without the need to recreate it physically. By knowing the physics behind it systems can be simulated without any material cost.
In maritime context it is of great interest to know how objects fully or partially submerged in water react in different situations. Hydrodynamic effects are imperative in simulations to get correct physical behavior in these environments. In a static system only buoyancy is present, but in a dynamic system there are also effects that are dependent on velocity, e.g. lift and drag. When the object is accelerating in a fluid it is also accelerating part of the fluid it is submerged in. This effect can be accounted for by adding a virtual mass to the object that, in accordance to Newton’s second law, will give a damping in the acceleration or deceleration of the object. This effect, called added mass, is dependent on how the object is moving and what part of it is submerged and may therefore change during the simulation. The computational cost for calculating these effect can be high, which becomes a problem when doing interactive simulations. At the same time it is important that the calculations are physically correct. The aim of this thesis is to find a method for interactive hydrodynamics, implement it and to validate the results through simulations.
1.2 r e l at e d wo r k
To calculate added mass Berklite [2] used a method later known as a boundary element method (BEM). By taking into account the depth of the water and the effects of the water surface he calculated the potential flow around an object to get the added mass coefficients. Even though he presented good results there are many variables to consider in the calcu-lations. A simpler way of using this method was presented by Ghassemi and Yari [4]. By having fully submerged objects in infinite water the flow would not change even if the object was moving close to the surface or near a wall.
Another common way for calculating added mass is by using strip theory. Under the as-sumption that the object is a slender body, i.e. long and thin, it can be divided into small cross sections and the added mass can be calculated for these strips separately. The values are then added to give the added mass tensor. Hem Lata and Thiagarajan [6] compared this
2 i n t ro d u c t i o n
method with two boundary element methods and discovered that the slender body method is faster and produces good results. On the other hand the boundary element methods works on arbitrary objects while strip theory only works for slender bodies.
A method for simulation that used both added mass and drag combined was tested by Weiß-mann and Pinkall [10]. Even if they only used viscous drag the results visually corresponded well to the experimental movement of the tested objects. The tests, however, only included objects that were falling or sinking freely in a fluid that they were completely submerged in. It did not handle partially submerged objects and the resulting tensor had to be scaled differently for different objects to get good results.
1.3 s c o p e a n d a p p roac h
The aim of this project is to find a method that enables interactive simulation of an arbi-trarily shaped object submerged in water and to implement and validate this method. The hydrodynamic effects of buoyancy, drag, lift and added mass will be taken into account.
1.3.1 Method
For calculating the drag and lift effects a method of numerical integration over triangles is used. For calculating the added mass the flow for each triangle around the object are precalculated with BEM assuming that the object is fully submerged in a infinite fluid. The added mass contribution from each submerged triangle will then be added to get the total added mass for that time step.
1.3.2 Trimesh
A usual way to represent an arbitrarily shaped object is to use a triangulated mesh, a trimesh. This means that the objects represented by a finite number of points in space creating a surface of triangles. Since drag and lift is calculated over triangle integration and the calculations of the added mass also is done over a triangulated grid the simulations will be done using trimesh objects.
1.3.3 Water
There is no simulation of water in this project. The water level is represented by a height field with a grid creating a triangulated mesh. The height field can simulate waves on the surface but the water will be considered still. The height field sole purpose is to determine what
1.4 outline 3
part of the simulated object is beneath the surface and at what depth. In the simulations all surfaces underneath the water level will be considered to be in contact with water, as illustrated in Figure 1. It will not be affected by the movement of objects submerged in it. When integrating the buoyancy force the height of the water surface is considered to be constant.
water surface
Figure 1: A bucket in water.
1.4 o u t l i n e
The methods used in this project is described in Chapter 2 together with a presentation of the hydrodynamic effects of buoyancy, drag, lift and added mass. There is also be a brief description on how these methods are used to implement these effects.
In Chapter 3 the details about the simulations are presented. The general algorithm for the simulations and the constants and coefficients used are described.
The results are presented in Chapter 4 for the effects of drag, lift and added mass com-pared to analytical results. There are also performance tests for both the precalculations and the simulation when running.
A discussion about the results and other effects of the implementation are held in
Chap-ter 5 and a brief summary of the conclusions in ChapChap-ter 6 together with suggestions for
2
T H E O RY2.1 n u m e r i c a l i n t e g r at i o n ov e r t r i a n g l e s
The simulation object consists of a triangulated mesh that has to be integrated over to get the effects from buoyancy, drag and lift. Since the triangles can be of arbitrary shape a change of variables is introduced so that the integration is over the unit triangle, T , instead. A generic triangle, T0, is defined by three vertices, v1, v2 and v3. By defining two vectors
p = v2− v1 and q = v3− v1 all points r(x, y, z) inside the triangle can be expressed as
r= v1+sp+tq, where s > 0, t > 0 and s+t ≤ 1, see Figure 2. Since the integration is done in the plane of a triangle, the problem is reduced to 2D.
x y z v2 v3 v1 p q (a) s t 1 1 (b)
Figure 2: A triangle in xyz-coordinates (a) and in st-coordinates (b).
The change of variables theorem in two dimensions states that
Z Z T0 f(x, y)dx dy= Z Z T f(s, t) ∂(x, y) ∂(s, t) ds dt, (1)
when changing coordinates in the region of integration [8]. For this problem ∂(x,y) ∂(s,t) =
||p × q||=2A, where A is the area of the triangle [7]. Eq. (1) can now be written as
Z Z T0 f(x, y)dx dy=2A Z Z T f(s, t)ds dt. (2)
The results from the integrals for the functions f(s, t)over the unit triangle T that is needed for buoyancy, drag and lift are given in Table 1.
6 t h e o ry
Table 1: Results from integrating different functions over the unit triangle.
f(s, t) 1 s t st s2 t2 s2t st2 s3 t3 value 12 16 16 241 121 121 601 601 201 201
2.2 b o u n da ry e l e m e n t m e t h o d ( b e m ) f o r l a p l ac e ’ s e q u at i o n
To determine the added mass of an object the potential flow φ around that object is calcu-lated using BEM. This method transforms a three dimensional problem into a calculation over a surface. A flow that is incompressible, irrotational and inviscid can be expressed with Laplace’s equation,
∇2φ=0. (3)
A Greens’ function [8] G(x, y) to this equation is a function solving the equation
∇2φ(x) =δ(x − y), (4)
while having the properties G(x, y) =G(y, x),
∇2G=δ(x − y), (5)
where δ is the Dirac delta function. Inserting this into Green’s second identity leads to
Z V (φ∇2G − G∇2φ)dV = Z S φ∂G ∂n − G ∂φ ∂n dS, (6)
where V is a region in space, S is its boundary and n is the unit normal pointing outwards [8]. Using the properties in Eq. (5) and choosing the integrating variable to be y the above expression becomes Z V φ(y)δ(x − y)dV = Z S φ∂G ∂n − G ∂φ ∂n dS. (7)
With the sifting property of the δ-function[8] the potential can be expressed as φ(x) = Z S φ∂G ∂n − G ∂φ ∂n dS, (8)
for any point x. The surface of integration is discretized to linear segments which in this case are triangles. The triangles are considered to have constant potential φ. This yields
φ(x) = N X j=1 φj Z ∂G ∂ndSj − N X j=1 " ∂φ ∂n j Z G dSj # , (9)
2.2 boundary element method (bem) for laplace’s equation 7
where N is the number of triangles. Expanding this to an equation system involving all centroid point of the triangles this can be written in matrix form:
φ=Dφ − S ∂φ ∂n , (10) where φ= [φ1, φ2, ..., φN]T, Dij = Z ∂G ij ∂nj dSj, where i, j =1, 2, ..., N Sij = Z GijdSj, where i, j=1, 2, ..., N ∂φ ∂n = ∂φ ∂n 1 , ∂φ ∂n 2 , ..., ∂φ ∂n N T . (11)
Knowing ∂φ∂n the it is possible for solve for φ using Eq. (10).
The Green’s function for Laplace’s equation in free space, i.e. unbounded space, is [5] G(x, y) =− 1
4π||x − y||. (12)
Inserting this into the expressions for Dij and Sij Eq. (11) yields
Dij =− 1 4π Z ∂ ∂n 1 ||x − y|| dSj = 1 4π Z 1 ||x − y||3(x − y)· ˆn dSj, Sij =− 1 4π Z 1 ||x − y||dSj. (13)
The integrals of the element in D and S are calculated with numerical integration, except for the diagonal entries of D that can not be solved since the integration involves the singularity when ||x − y||=0 and the function to be integrated becomes 00. The entries of D can instead be expressed as [4] Dij = 1 2 when i=j, 1 4π R 1 ||x−y||3(x − y)· ˆn dSj when i 6=j. (14) 2.2.1 Numerical integration
The integrands in Eq. (13) are complicated functions that are not multinomials in s and t and therefore the method of integration described in section 2.1 is not used. They are
8 t h e o ry
however still expressed as functions of s and t which enables the use of Eq. (2), that is rewritten so the region of integration is mapped onto a rectangle.
2A Z Z T f(s, t)ds dt=2A 1 Z 0 1 Z 0 f(u(1 − t), t)(1 − t)du dt. (15)
Using the Gaussian integration formula [1], 1 Z 0 1 Z 0 f(u, t)du dt ' 1 16 16 X i=1 f(ui, ti), (16)
the integrand is evaluated at the following 16 (u,t)-coordinates:
(u1, t1) = 1 4 + 1 4√3, 1 4+ 1 4√3 , (u2, t2) = 1 4 + 1 4√3, 1 4− 1 4√3 , (u3, t3) = 1 4 − 1 4√3, 1 4+ 1 4√3 , (u4, t4) = 1 4 − 1 4√3, 1 4− 1 4√3 , (u5, t5) = 3 4 + 1 4√3, 1 4+ 1 4√3 , (u6, t6) = 3 4 + 1 4√3, 1 4− 1 4√3 , (u7, t7) = 3 4 − 1 4√3, 1 4+ 1 4√3 , (u8, t8) = 3 4 − 1 4√3, 1 4− 1 4√3 , (u9, t9) = 3 4 + 1 4√3, 3 4+ 1 4√3 , (u10, t10) = 3 4 + 1 4√3, 3 4− 1 4√3 , (u11, t11) = 3 4 − 1 4√3, 3 4+ 1 4√3 , (u12, t12) = 3 4 − 1 4√3, 3 4− 1 4√3 , (17)
2.3 hydrodynamic effects 9 (u13, t13) = 1 4 + 1 4√3, 3 4+ 1 4√3 , (u14, t14) = 1 4 + 1 4√3, 3 4− 1 4√3 , (u15, t15) = 1 4− 1 4√3, 3 4+ 1 4√3 , (u16, t16) = 1 4− 1 4√3, 3 4− 1 4√3 . (18) 2.3 h y d ro dy n a m i c e f f e c t s 2.3.1 Buoyancy
An object submerged in a fluid is affected by the pressure caused by the fluid surrounding it. When the fluid and the body is at rest the pressure acting on it, known as the hydro-static pressure, can be expressed as ps = ρgh where ρ is the density of the fluid, g is the
gravitational constant and h is the vertical distance to the surface. To find the force acting on the body the pressure is integrated over the surface of the object:
Fb =−
Z Z
psˆn dS=−ρg
Z Z
h ˆn dS, (19)
where ˆn is the normal pointing outwards [3]. Since pressure is increasing with depth the
object will experience a larger pressure from underneath than from above. This results in a net force pointing upwards called buoyancy.
The torque produced by a force F is given by
τ = (r − rcm)× F, (20)
where r is the position that the force is acting on and rcm is the position of the center of
mass [9].
To get the buoyancy force in Eq. (19) the pressure have to be expressed in terms of the position r. If ξ(r) is the height of the surface at a given position the vertical distance h between the point and the surface can be expressed as h = ξ − r · ˆz. If ξ is constant, i.e.
10 t h e o ry
the water level is flat, combining Eq. (19) with Eq. (2) and making the variable change
r=v1+sp+tq leads to Fb =−2Aρg ˆn Z Z ξ −(v1+sp+tq)· ˆz ds dtS =−2Aρg ˆn Z Z ξ − v1· ˆz − s(p · ˆz)− t(q · ˆz)ds dt =−2Aρg ˆn Z Z h0− sa − tb ds dt, (21)
where a=p · ˆz, b=q · ˆz and h0 =ξ − v1· ˆz is the vertical distance between the vertex v1 and the surface. From Table 1 the final expression for the buoyancy force becomes
Fb =−2Aρg h 0 2 − p · ˆz 6 − q · ˆz 6 ˆ n=−Aρg h0− (p+q)· ˆz 3 ˆ n. (22)
Combining Eq. (20) and Eq. (21) the torque produced by this force is τb = Z Z (r − rcm)× −ρg(ξ − r · ˆz)n dSˆ =−2Aρg Z Z (v1+sp+tq − rcm)×(h0− sa − tb)ˆn ds dt =−2Aρg Z Z (h0− sa − tb) [(v1+sp+tq − rcm)× ˆn] ds dt =−2Aρg Z Z (h0− sa − tb)(v1− rcm)× ˆn+ (h0s − s2a − stb)p × ˆn+ (h0t − sta − t2b)q × ˆn ds dt. (23)
The net torque acting on the object is then τb =−2Aρg h 0 2 − a+b 6 (v1− rcm)× ˆn+ h 0 6 − a 12− b 24 p × ˆn+ h 0 6 − a 24− b 12 q × ˆn . (24)
2.3.2 Drag & Lift
If the object is moving the system will no longer be static and there will be additional forces. The velocity of the object will give rise to drag and lift forces. The pressure drag, Fp, and the lift, Fl, is caused by the dynamical forces from the fluid and the viscous drag, Fv, is
caused by the frictional shear stress from the viscosity of the fluid. The expressions for these forces are given by
Fp =− Z Z 1 2ρ Cp(v · ˆn) 2sgn(v · ˆn)n dS,ˆ (25a) Fl =− Z Z 1 2ρ Cl(v · ˆtv) 2ˆn dS, (25b) Fv =− Z Z 1 2ρ Cv(v · ˆtv) 2ˆt vdS, (25c)
2.3 hydrodynamic effects 11
where v is the velocity at the objects surface, ρ is the density of the fluid, ˆtv is the direction
of the velocity tangent to the surface and Cp, Cl and Cv are the drag and lift coefficients [7].
For the case of pressure drag, Eq. (25a), there is a signum function that for some triangles will be constant and for others will have different signs on different parts of the triangle. The triangles with two different signs are clipped into three new triangles with a constant sign to avoid the integration over the signum function, see Figure 3. Clipping is also done at the water surface to make sure that only submerged triangles are integrated.
+ −
− −
+
Figure 3: Sign clipping of a triangle.
With the signum function removed from the equation the expressions in Eq. (25) all are of the form
Fi =−
Z Z 1
2ρ Ci(v · ˆk)
2ˆl dS, (26)
where i=p, l, v for the different kind of forces and ˆk and ˆl are replaced by the unit vectors ˆtv and/or ˆn. To get this in terms of s and t, v is written as v = vcm+ω ×(r − rcm) = vcm+ω ×(v1+sp+tq − rcm), where ω is the angular velocity of the object [9]. This
expression inserted in v · ˆk yields
v · ˆk= [vcm+ω ×(v1+sp+tq − rcm)]· ˆk
= [vcm+ω ×(v1− rcm)]· ˆk+s(ω × p)· ˆk+t(ω × q)· ˆk
=γ+sα+tβ,
(27)
where the constants γ, α and β are γ = [vcm+ω ×(v1− rcm)]· ˆk
α = (ω × p)· ˆk
β = (ω × q)· ˆk.
12 t h e o ry
Eq. (26) now becomes
Fi =−2A 1 2ρ Ci ˆlZ Z (γ+sα+tβ)2ds dt =−Aρ Ciˆl Z Z γ2+2αγs+2γβt+2αβst+α2s2+β2t2ds dt. (29)
With values from Table 1 the net force is
Fi =−Aρ Ci γ 2 2 + γα 3 + γβ 3 + αβ 12 + α2 12 + β2 12 ! ˆl. (30)
From Eq. (29) and Eq. (20) the torque acting on the object is
τi =−AρCi
Z Z
(v1+sp+tq − rcm)×(γ2+2αγs+2γβt+2αβst+α2s2+β2t2)ˆl ds dt.
(31) Determining the integral with Table 1 the net torque becomes
τi= γ2 2 + 2αγ 6 + 2γβ 6 + 2αβ 24 + α2 12 + β2 12 ! (v1− rcm)׈l + γ 2 6 + 2αγ 12 + 2γβ 24 + 2αβ 60 + α2 20 + β2 60 ! p ׈l + γ 2 6 + 2αγ 24 + 2γβ 12 + 2αβ 60 + α2 60 + β2 20 ! q ׈l. (32) 2.3.3 Added Mass
When the velocity of an object that is moving through a fluid is changed the velocity for part of the fluid surrounding it also changes. The effective mass of the object has therefore increased. The added mass depends on the shape of the object, what part of the object is submerged and the direction of motion [3]. The direction of motion not only includes translational motion in x, y and z-direction but also three rotational motions roll, pitch and yaw explained in Figure 4, which means that also the inertia tensor of the object is affected. The added mass becomes a 6 × 6 matrix A, where the rows m and the columns n represents the different degrees of freedom.
Assuming that the flow is irrotational, inviscid and incompressible the velocity potential φ fulfills Laplace’s equation in Eq. (3), where the velocity, u, of the fluid can be expressed as u = ∇φ [3]. When calculating the potential flow φ the object is assumed to be fully submerged in infinite water with no floor and no surface. However, when using the flow to
2.3 hydrodynamic effects 13 x y z roll yaw pitch
Figure 4: The six degrees of freedom.
calculate the added mass the object may be partially submerged.
By analyzing the kinetic energy of the accelerated fluid [4] the added mass coefficients amn
can be expressed as amn =−ρ Z S ∂φ ∂n m φndS, (33)
where m and n represents the degrees of freedom, ρ is the density of the fluid and n is the surface normal pointing outward. The potential φ is calculated using the BEM described in section 2.2. Since the object is assumed to be submerged in an infinite fluid, the free space Green’s function in Eq. (12) is used. That means that the only boundary in this problem is the one of the submerged object. The condition on this boundary is expressed as
∂φ
∂n =−v · n, (34)
where v is the velocity of the object. This means that the fluid can not flow through the surface of the submerged object. By separating the six degrees of freedom and considering a unit potential we get the following conditions:
∂φ ∂n m = ni, for m = 1,2,3 and i = 1,2,3 (r × n)i, for m = 4,5,6 and i = 1,2,3 (35)
where m is the degrees of freedom, i is the vector component and r is the position on the boundary relative the center of mass [4]. Using Eq. (10) for each degree of freedom separately the flow is solved for all triangles. Since φ and n is constant over the triangles Eq. (33) becomes mmn=ρ N X k=1 ∂φ ∂n k m φknAk, (36)
14 t h e o ry
where Ak is the area of the triangle. If the triangle is partially submerged only the
3
S O F T WA R E I M P L E M E N TAT I O NThe simulations are made with the help of the physics engine AgX Dynamics, a product developed by Algoryx Simulations AB. AgX does collision detection between the geometries and applies the appropriate response to handle collisions and contacts. The implementation made in this project provide a way to handle the case when there is an intersection between an object and water. AgX also do the time integration and the rendering, except for the visualization of the pressure field that is done using shaders. The main code in this project is written in C++ and the shaders are written in OpenGL Shading Language (glsl).
3.1 a l g o r i t h m d e s c r i p t i o n
Algorithm 1 The main simulation algorithm
1: Scene setup
2: while the simulation is running do
3: Collision detection
4: if there is an intersection between the object and water then
5: if the flow potential for the objects is not stored then
6: Calculate φ (Eq. (10)) and store it.
7: end if
8: Calculate pressure and depth at all vertices.
9: Send values of pressure and depth to shader to enable pressure field rendering. 10: for all triangles do
11: if the height of all three vertices are positive then
12: Set the submerged area to be zero and go to the next triangle.
13: end if
14: if the height of the vertices are of different signs then
15: Perform surface clipping.
16: for the new triangles that are underneath the water surface do
17: Add their area to the submerged area of the triangle. 18: Perform sign clipping of the dot product v · n. 19: Calculate hydrodynamic forces and torques. 20: Apply forces and torques to the object.
21: end for
16 s o f t wa r e i m p l e m e n tat i o n
22: else
23: The triangle is completely submerged. Calculate the area. 24: Perform sign clipping of the dot product v · n.
25: Calculate hydrodynamic forces and torques. 26: Apply forces and torques to the object.
27: end if
28: Add the triangles contribution to the added mass with the calculated area.
29: end for
30: Apply the added mass tensor to the object.
31: end if
32: Render the scene. 33: Time integration. 34: end while
3.2 p r e s s u r e f i e l d
To visualize the pressure field acting on the submerged object the pressure at each vertex is transformed to a color in a graphics shader. Vertices above the water level will have the original color of the object while vertices underneath the waterlevel will be given a color according to a color scale between pmin and pmax as given in Figure 5. The minimum
value is zero and the maximum value is chosen differently depending on what object is being simulated. The color of a point on a triangle will be determined by an interpolation between the colors of the three vertices. If the triangle is partially submerged the shader creates new vertices and new triangles so all triangles are either fully submerged or not submerged at all.
pmin pmax
Pa
Figure 5: The color gradient used to visualize the pressure field
3.3 c o n s ta n t s a n d c o e f f i c i e n t s
The constants and coefficients of the simulations are given in Table 2. The density is the tabulated value for water [9] and the the gravitational constant is a default value used in AgX Dynamics. The viscous drag coefficient is dependent on the Reynolds number, i.e. ratio between inertial and viscous forces, and can be neglected for high velocities. The Reynolds number is dependent on velocity, characteristic length of the object and the kinematic
vis-3.3 constants and coefficients 17
cosity of the fluid and therefore the viscous drag coefficient has a simple velocity dependence and is only present at low velocities.
Table 2: Values for constants and coefficients.
Constant/coefficient Notation Value Unit
Density of water ρ 1000 kg/m3
Gravitational constant g 9.80665 m/s2
Pressure drag coefficient Cp 0.6316
-Lift coefficient Cl 0.01
-Viscous drag coefficient Cv 0.1 (v<1)
-0.2-0.1v (1<v<2) 0 (v>2)
4
R E S U LT STo test the validity of the methods the following tests are made to analyze the effects of hydrodynamics and the performance of the implementation:
• Spheres and cubes of different tessellations are dropped from a height. From this simulation the buoyancy force is compared to analytical values and the angular velocity is analyzed.
• The pressure drag for different objects is calculated for a several objects to determine the value for the drag coefficient.
• A sphere held underneath the water level is given a torque to start a rotation. The angular velocity for an object with viscous drag is compared to one without viscous drag.
• A sphere is dragged horizontally by a given force. The height of an object with lift enabled is compared to one with lift disabled.
• The added mass for spheres and cubes of different tessellations are computed and compared to analytical values.
• The time requirement for the precalculations of the potential flow is measured for objects of different size.
• The time requirement for the runtime calculations for hydrodynamics is measured for objects of different size.
4.1 b u oya n c y
Two spheres with a radius of 1 m and a density of 550 kg/m3 is simulated with two different tessellations, one with 480 triangles and the other with 4900 triangles. The only hydrody-namic effect is the buoyancy while the other effects are disabled. In the beginning of the simulation the sphere is dropped from a height of 2 m and then falls towards the water surface. When hitting the surface a buoyancy force appears and is increasing with depth. The force acting on the object at each time step is plotted against the height and compared to the analytical value. As seen in Figure 6 the buoyancy force for the sphere with finer tessellation is closer to the analytical result with an error of 0.36 % than the sphere with a coarse tessellation that has an error of 3.49 %.
20 r e s u lt s
Figure 6: Buoyancy force for two spheres with 480 and 4900 triangles.
The same test is done on two cubes with side lengths of 2 m and the same density as for the spheres. The first cube has the simplest representation possible with 2 triangles per side and a total of 12 triangles while the other has a number of 128 triangles per side resulting in a total of 768 triangles. From the values in Figure 7 the errors are 1.85 · 10−14 % for the cube with 128 triangles and 5.21 · 10−7 % for the cube with 768 triangles.
(a) (b)
Figure 7: Two cubes with 12 triangles (a) and 768 triangles (b) compared to the analytical result.
To test the torque part of the buoyancy the angular velocities for a cube that is dropped from 2 m above the surface is plotted against time for two different tessellations. The angular velocity for both cubes varies over time and the more tessellated cube have a generally larger angular velocity than the coarse tessellated cube, as seen in Figure 8.
4.2 pressure drag forces 21
(a) (b)
Figure 8: Angular velocity of two cubes with 12 triangles (b) and 768 triangles (b).
This test is also done with two spheres of different tessellation. In Figure 9 is it shown that the angular velocity of a sphere is increasing with time and that the angular velocity is larger for a sphere with finer tessellation.
Figure 9: Angular velocity for two spheres with 480 and 4900 triangles.
4.2 p r e s s u r e d r ag f o rc e s
The drag is calculated for a range of velocities and compared to tabulated values to see what drag coefficient would give the same force. For the tabulated value both pressure and viscous
22 r e s u lt s
drag is included. In this project the viscous drag is very small so the calculated coefficient would correspond to the pressure drag coefficient. This test is done for two different models and with two different directions of motion. The four test setups can be seen in Figure 10.
(a) (b) (c) (d)
Figure 10: The four test setups.
In Table 3 it is shown that the pressure drag coefficient in the simulation Cpis not the same
coefficient as the tabulated coefficient Cd.
Table 3: Drag coefficients for different setups
Test Cube (a) Cube (b) Cylinder (c) Cylinder (d)
Calculated Cp 0.59 0.90 0.57 0.47
Tabulated Cd 1.07 0.81 0.68 0.85
The chosen Cp is the average value from the four calculated values in this test which is Cp =0.6316. When testing the drag force for this value the plots in Figure 11 is obtained,
where the force is calculated for different tessellations and compared to the tabulated value.
4.3 v i s c o u s d r ag f o rc e s
To test the viscous drag a sphere is placed 3 m underneath the water surface and a torque of 50,000 Nm is applied in the beginning of the simulation. To minimize the disturbance from the numerical errors described above the sphere is held in place by a constraint during the simulation. This is first done with a sphere with buoyancy, drag and lift enabled and is then repeated with the viscous drag disabled. For the sphere with the viscous drag enabled the rotation is decreasing more rapidly then for the sphere without viscous drag, as seen in Figure 12.
4.4 lift forces 23
(a) (b)
(c) (d)
Figure 11: Difference between tabulated values and simulated values for the drag force for a cube (a), a tilted cube (b), a standing cylinder (c) and a lying cylinder (d).
4.4 l i f t f o rc e s
To see the effects of lift a sphere is dragged through water with a force of 10000 N. This in done two times, first with lift enabled and then with no lift. The difference can be viewed in Figure 13, where the sphere effected by lift has a higher position then the one without lift.
24 r e s u lt s
Figure 12: Angular velocity over time for two spheres, one with viscous drag and one without.
Figure 13: Height above the water level for spheres with and without lift.
4.5 a d d e d m a s s
The added mass tensor A for a sphere with a radius of 1 m and 4900 triangles is calculated to be A= 2087.49 0.00 0.03 −0.37 0.00 0.02 0.00 2080.71 −0.00 0.00 0.78 0.00 −0.04 −0.00 2087.49 −0.02 −0.00 −0.37 −0.38 0.00 −0.02 0.02 0.00 0.00 0.00 0.77 −0.00 0.00 0.04 −0.00 0.02 0.00 −0.38 −0.00 −0.00 0.02 (37)
4.5 added mass 25
with BEM. The analytical values are amn=2094 for i=j= 1, 2, 3 and zero for the other
entries. For a cube with length 2 m the approximate value for the added mass is amn =5600
for the same entries. In Table 4 the added mass for a sphere and a cube are compared to the analytical values for the first three diagonal values. For both the sphere and the cube the finer tessellations produces results closer to the analytical value.
Table 4: Added mass for a sphere of radius 1 m.
Object Entry Analytical value BEM value Relative error (%)
Sphere 4900 triangles a11 2094 2087 -0.33 a22 2094 2081 -0.65 a33 2094 2087 -0.33 Sphere 480 triangles a11 2094 2059 -1.71 a22 2094 2059 -1.68 a33 2094 1993 -4.82 Cube 768 triangles a11 5600 5207 -7.0 a22 5600 5208 -7.0 a33 5600 5206 -7.0 Cube 12 triangles a11 5600 6490 15.9 a22 5600 6492 15.9 a33 5600 6491 15.9
26 r e s u lt s
4.6 p e r f o r m a n c e
To test the performance for the precalculations of the flow φ the time requirement is mea-sured for a range of different models with various number of triangles. The most time consuming part of these calculations, that are performed only once in the beginning of the simulation, is solving the matrix equation, as seen in Figure 14. The rest of the time is used setting up the matrices by calculating the integrals in Eq. (13).
Figure 14: Time required for precalculation of flow for different objects.
The time it takes for the hydrodynamics calculations is measured during a simulation of 1 s in 60 Hz. The test is repeated for various models that are held under water with a constraint. In Figure 15 the measured time is plotted against the number of triangles. The time for the added mass calculations are marked in the figure, the rest of the time is used for surface and sign clipping and to calculate the effects of buoyancy, drag and lift.
4.7 demonstration 27
Figure 15: Time required to perform hydrodynamic calculations during a simulation of 1 s in 60 Hz.
4.7 d e m o n s t r at i o n
To show how an object behaves in the simulation a teapot is dropped from a height. Shots from the scene is saved every 0.3 s in the simulation and can be seen in Figure 16. Above the water level the objects color is visible. Underneath the water level the color is determined by the pressure acting on the object.
28 r e s u lt s
(a) (b) (c)
(d) (e) (f)
(g) (h) (i)
(j) (k) (l)
5
D I S C U S S I O NIn the buoyancy test it is apparent that the force on the sphere is dependent on the tes-sellation while that of the cube is almost unaffected. This can be explained by the fact that the shape of the cube is represented in an exact way, while for the sphere the shape is an approximation that gets better with finer tessellation. The results for the method of triangle integration is therefore only dependent on how well the tessellated object represents the actual shape. The small errors for the cube is a numerical error in the calculations and increases therefore with number of triangles, but is still very small.
A sphere or a box with a flat face down, dropped from a height into water should not get any torque from the buoyancy and should therefore not gain any angular velocity. The errors in the torques are small numerical errors that increases with the number of triangles. For the coarse tessellated cube the error is small and mostly resembles a noise. For the cube with finer tessellation the angular velocity is larger but still seems to stay at a stable level. A cube that receives a torque and gains angular velocity will in the next time step be slightly tilted. This will produce a buoyancy torque in the other direction stabilizing the position of the cube. For the sphere however, there is no such stability and the angular velocity will increase every time a small torque is added. A sphere with coarse tessellation has sharp corners that gives it some stability which is an explanation to why the angular velocity is larger for the sphere with finer tessellation. In Figures 8 and 9 there are short periods of time where the angular velocity does not change. For the cube with a coarse tessellation these intervals is when the cube is not in contact with water at all. For the cube with finer tessellation there is also a constant angular velocity when the cube is fully submerged. This means that the error becomes a problem only when the object is partially submerged. This test is done with only buoyancy enabled. It is however possible that this error would be less of a problem when the other effects, such as viscous drag, also is present. It is also important to note that when the object is affected by torque from other effects or collisions the buoyancy torque might become negligible.
From the pressure drag test it is clear that the tabulated value for the drag coefficient is not the same as the coefficient in the simulation. For three of the test cases the value for the simulation is smaller, which could be explained by the fact that the simulation adds forces for both the front and the back of the object while the tabulated value is based on the frontal area of the object. For the case of the tilted cube however, the value for the simulation is larger than the tabulated value. The value used in the simulations for an ob-ject of arbitrary size is the average of the values for these test cases. This is a very crude approximation that clearly does not work for all types of objects. It is also important to
30 d i s c u s s i o n
note that the coefficients vary, not only with geometry, but also with the Reynolds number of the system, which is determined from the velocity, a characteristic length of the object and the kinematic viscosity of the fluid.
Both the viscous drag and the lift behaves as expected from the tests. The lift coefficient is however determined by subjective tests on how the object should behave in water. For both the lift and the viscous drag coefficients to be verified more test needs to be done. At high velocities there may be some instability caused by too strong forces and torques acting on the object.
The added mass for a sphere and a cube agrees well with the analytical value, especially when there is a fine tessellation of the object. In the case of added mass the result is depen-dent on the tessellation of an object and the error decreases when the number of triangles increases. The effect from assuming complete submersion when calculating the flow velocity potential and only count the contribution from the submerged part when calculating added mass in unknown. The flow has the strongest dependence from the neighbor elements so the errors would be mainly from near the water surface. When adding up the contributions this error should be small.
When doing performance tests on the precalculations for the flow potential the time con-sumption increases rapidly with number of triangles. Most of that time is required to solve the matrix equation since it requires a LU-factorization of a dense matrix which has a complexity of O(n3). The matrix solver used in this project could easily be replaced which could speed up the calculations. These calculations could also be done separately from the simulation and the data could be loaded from a file. Then the calculations could be done only when using the model for the first time instead of in the beginning of each simulation. The test for simulation time is also increasing with number of triangles but with a linear dependence. The time that is required to calculate the added mass tensor is very small compared to the other calculations thanks to the precalculations of the potential flow. For the other calculations it is possible to optimize the code to receive even better results. The code could also be parallelized which could decrease the running time considerably.
The demonstration in Figure 16 shows the pressure acting on the object when it is submerged. In (a) and (e) the teapot has similar submersion but in (a) there is a higher pressure at the bottom shown by the red color. This is because the teapot has a higher velocity from falling and is affected by pressure drag. In (e) the teapot has a lower velocity in the other direction and the pressure distribution seen here is mostly from the hydrostatic pressure.
d i s c u s s i o n 31
The water in the simulations is represented by a height field determining the level of the water surface. A limitation in using this representation of water is that it does not handle surfaces that is underneath the water level but sheltered from water by, for instance, the hull of a ship, as shown in Figure 1. All surfaces underneath the surface level will be considered to be in contact with water and will experience a pressure from it. It is therefore important to consider models carefully.
6
C O N C L U S I O N SThe method presented in this report works well for interactive simulations of arbitrarily shaped objects. The tests produced results close to the analytical values and the time con-sumption in each time step is reasonable with a linear dependece on the number of triangles. There is a tessellation dependence meaning that the accuracy increases with number of tri-angles of the object. For the numerical integration over tritri-angles this depence is only present when the tessellation affects how well the shape is represented. For a cube, that has an ex-act representation with only 12 triangles, a finer tessellation does not give a more accurate result. For BEM however, the accuracy will increase with finer tessellation for all objects.
6.1 f u rt h e r wo r k
The methods implemented in this project produced good results, but there is more to be done on the subject. The following points are suggestions for further work.
• The coefficients for lift and drag could be tested more to be better adjusted. There should also be tests checking the stability of these effects at high velocities.
• More simulations could be done testing the numerical error that appeared in the implementation of buoyancy. There should also be a way of handling this error if it becomes a problem in the simulations.
• In this project the added mass was only tested for fully submerged objects. More tests could be done with objects that are partially submerged.
• The code could be further optimized to give both shorter precalculations and to enable objects of even finer tessellation to be interactively simulated.
• Simulation of water that could include waves, vorticity and wakes.
R E F E R E N C E S
[1] Whye-Teong Ang. A beginner’s course in boundary element methods. Universal Pub-lishers, 2007.
[2] Ronald Betts Berklite. Added mass of submerged objects of arbitrary shape. Master’s thesis, Naval postgraduate school, Monterey, California, 1972.
[3] James W. Daily and Donald R. F. Harleman. Fluid Dynamics. Addison-Wesley, 1966. [4] Hassan Ghassemi and Ehsan Yari. The added mass coefficient computation of sphere, ellipsoid and marine propellers using boundary element method. Polish Maritime Re-search, 18(1):17–26, April 2011. ISSN 1233-2585.
[5] G. Steven Gipson. Boundary Element Fundamentals - Basic Concepts and Recent Development in the Poisson Equation. Computational Mechanics Publications, 1987. [6] W. Hem Lata and K. P. Thiagarajan. Comparison of added mass coefficients for a
floating tanker evaluated by conformal mapping and boundary element methods. In 16th Australasian Fluid Mechanics Conference (AFMC), pages 1388–1391, December
2007.
[7] Claude Lacoursière. Computing forces between fluids and rigid bodies.
[8] Jerrold E. Marsden and Anthony J. Tromba. Vector Calculus. W. H. Freeman and Company, 1988.
[9] Carl Nordling and Jonny Österman. Physics Handbook for Science and Engineering. Studentlitteratur, 2006.
[10] Steffen Weißmann and Ulrich Pinkall. Underwater Rigid Body Dynamics. ACM Trans. Graph., 31(4):104:1–104:7, July 2012. ISSN 0730-0301.
