Chapter 1 Introduction
1.3 Thesis Organization
The following chapters are organized as follows. In Chapter 2, we will introduce two mainstream directions of fluid simulation in computer graphics, one is small volume liquid simulation and the other is ocean simulation. In Chapter 3, we introduce our proposed tsunami model and how we simulate the propagation of tsunami. Then we introduce how to generate light reflection on the ocean surface. The implementation and results will be demonstrated in Chapter 4. Finally, Conclusion and future works will be given in Chapter 5.
Chapter 2
Related Works
2.1 Fluid Simulation
Currently, the simulation of complex water effects is often based on three-dimensional Navier-Stokes [1] equations. Through 3D Navier-Stokes equations, we can simulate churn, splash and sprinkle in the 3D space in detail. With advanced physically-based rendering system, we can simulate a really vivid scene. But the drawback of methods mentioned above is computational-consuming. It is hard to simulate in real-time. Foster and Fedkiw [2] made significant contributions to the simulation and control of 3D fluid simulations through the introduction of a hybrid liquid volume model combining implicit surfaces and mass-less marker particles.
Considering the simulation of wave motions of large water volume by using 3D Navier-Stokes equations, in order to reduce the computational cost, hybrid models are popular. Losasso et al. [7] introduced a way that coarsens away from the water surface with an octree data structure. But this doesn’t make full use of the highly accurate MAC grids and numerical dissipation will increase. Besides large octree cell can not represent bottom topography. There are also methods based on height fields such as 2D shallow water equations [8]. Though this approach can capture effects of complex bottom topography rather well, it does not support overturning or other 3D behaviors.
Irving and Guendelman [5] proposed a water volume model which makes a combination of uniform grids and tall cells. In uniform grids near the surface, the
advection is calculated by 3D Navier-Stokes equations; in tall cells, advection equations are simplified. Pressures and velocities are interpolated within the cells.
This hybrid method can keep details of the water surface and reduce the computation cost under water surface.
2.2 Ocean Simulation
The classic reference on ocean waves rendering is the work of [3]. Fournier and Reeves presented a simple model for the surface of the ocean. It is based on the Gerstner model where particles of water describe circular or elliptical stationary orbits and the foam generated by the breaker is modeled by particle systems. The wave shape is modeled as follows,
x = x 0 + r ∗ sin( κ x 0 − ω t ) (1a) z = z 0 − r ∗ cos( κ x 0 − ω t ) (1b) where t is the time, z is the vertical axis and is the particle location at rest, r is a distance from the center of a circle of radius
) , ( x 0 z 0
κ
1 rolling over a line.
Figure 2 shows the trochoid of Gerstner wave model. The surface is modeled as a parametric surface and can be rendered by traditional rendering methods. By combining more sine, cosine functions, the wave surface will be less regular.
Figure 2.1: The trochoid of Gerstner wave model
Gonzato and Bertrand [4] proposed a complete geometrical model of the ocean waves accounting for refraction, diffraction, reflection, transmission and multi-wave
trains. Then they rendered the scene by using a particular illumination model and displacement mapping texture to deal with multi-wave trains. It makes physically-realistic result but not in real-time.
However oceanographers do not take Gerstner wave as a realistic model of the ocean. Instead, they introduce statistical models which combine experiment observations. In statistical model, the wave height is formulated as a function , where x is the horizontal position and t is the time. Jason L. Mitchell [6] present a multi-band Fourier domain approach to synthesize the ocean surface and rendering the waves on GPU entirely. The height of the water at location x at time t is :
) , ( tx f
=
∑
k
eikx
t k h t
x
h ~( , )*
) ,
( (2)
where k is a 2D vector with components (kx,ky), kx =2πn/Lx,ky =2πm/Ly, and n and m are integers with bounds –N/2 < n < N/2 and –M/2 < m < M/2. Then IFFT (Inverse Fast Fourier Transform) will generate a height field at discrete point
, and we can finally construct ocean water.
) / , /
(nL N mL M
x= x y
Yaohua Hu et al. [12] presented a system that can render calm ocean wavers with sophisticated lighting effects at a high frame rate. They construct the wave geometry as a displacement map with surface detailed by a dynamic bump map.
Xudong Yang et al. [11] presented a multi-resolution mesh model of the ocean surface based on a straightforward terrain level of detail scheme, Tiled Quad-tree. By this scheme, the ocean surface can be extended limitless and accelerated by GPU.
Chapter 3
Tsunami Wave Model
In the following sections, we will describe how we construct the tsunami wave model and render the complete ocean scene. In section 3.1, we describe the basic equations for our simulation system. Then we demonstrate how we obtain the linear system from equations we described above in section 3.2. In section 3.3, we introduce a numerical method to solve the linear system. Ocean surface construction is described in section 3.4. Finally, we show how we render the ocean scene in section 3.5.
3.1 Shallow Water Equation
Two-dimensional depth-averaged equation is generally called shallow water equation (SWE) which is integrated from 3-D Navier-Stokes equations. The depth is integrated from the bottom to the top surface. The SWE describes the evolution of a hydrostatic homogeneous, incompressible flow on the surface of the sphere. This equation is accurate when the ratio of the vertical scale to the ratio of the horizontal scale is small. SWE can be used in atmospheric and oceanic modeling, but are much simple than its primitive equations ( Navier-Stokes equations ).
3.1.1 Derivation of Shallow Water Equation
Assuming a column of fluid of height h, and a base area A. Then the total mass of fluid at the base over the area A is
m=hAρ, (3) and it corresponds to a force in the direction of gravity such as
mg
f = . (4) Then we can find out that the pressure at the base area is therefore
p= f /A=hAρg/A=hρg. (5) Hence we can rewrite the pressure term in (5) into
∇p=∇hρg. (6) In large-scale motion of the fluid, “shallow” means that there is only little
variation in y direction, so the vertical acceleration can be ignored. Now the fluid is said to be in hydrostatic equilibrium, and it implies that the Newton’s second law of motion can be simplified as
U U g h U Fe
where is the external forces, , and the changes in the horizontal plane are nonzero. So we can reduce (7) to the following two equations:
Fe U =[u,v,w]T
The mass preservation in terms of the height is found by make an integration on z, so surface, hence
w0 wh
w0 =0, (12)
dt
wh = dh. (13)
Then insert (12), (13) into (11), we can find that
y
↓
Then the SWE can be obtained by using , and assuming inviscid fluid, that is
3.2 Implicit Semi - Lagrangian Time Integration Method
3.2.1 Semi - Lagrangian Method
The semi-Lagrangian method can be viewed as a combination of Eulerian method and Lagrangian method. In the Eulerian method, we describe one kind of characteristics of the fluid at a fixed point in the world space by a function E(x,y,z,t), and we observe the world evolving around the fixed point. In the Lagrangian method, we follow a particle moving along its trajectory, so the physical factors which at certain point X in the world space at time T are only affected by their initial positions and the time variable T.
In the semi-Lagrangian scheme, at every time step the grid points of the numerical mesh are representing the arrival points of backward trajectories at the future time steps. The point reached during this backward tracking defines where a small space was at the beginning of the time step. The particle is subjected to different physical forces during its traveling. All predictive variables on this departure point are then found by interpolating the values on the grids at the previous time steps, as
shown in Figure 3.1.
Figure 3.1: The Lagrangina derivation is approximated along particle trajectories
The advantage of semi-Lagrangian method is that it allows the use of large time step without limiting the stability and it is superior in performance. The limitations for stability are that trajectories of fluid particles should not cross and not overtake another. Hence, the choice of time step in the semi-Lagrangian method is only limited by numerical accuracy.
3.2.2 Applying semi-Lagrangian method to SWE
To show how the SWE can be integrated by the semi-Lagrangian method, we first consider equations (16) and then we replace the substantial derivates in (16) with ordinary derivates
x Then, by using forward differences, we can estimate the ordinary derivates by tracking a particle at grid point (i, j),
t
t at grid point (i, j), based on an assumption that these current values are also good approximated results in the previous steps. The departure point of hypothetical particle at (i, j) can be estimated as follows :
by bilinear interpolation,
(
a)(
b)
u( ⎣ ⎦ ⎣ ⎦
x y) (
a b)
u( ⎡ ⎤ ⎣ ⎦
x y)
abu( ⎡ ⎤ ⎡ ⎤
x y) (
a)
bu( ⎣ ⎦ ⎡ ⎤
x y)
u~ = 1− 1− ~ , ~ + 1− ~ , ~ + ~ , ~ + 1− ~ , ~ , (20)
where a=x~ −
⎣ ⎦
x~ , b=~ −y⎣ ⎦
~y . Then v~ and h~are also similar. Now we discretize (17) and combine with (18), then we obtain,
x
Then we calculate the derivative of (21a) with respect to x and (21b) with respect to y, we obtain,
)
after rearranging the terms,
))
Then we apply first order and second order central differences for the spatial derivates, we can get the following equation,
)) system as a linear system for
)
The algorithm is shown in Figure 3.2:
Figure 3.2: Algorithm for the SWE
3.3 Conjugated Gradient Method
In mathematics, the conjugated gradient method is an algorithm for the numerical solution of particular system of linear equations, namely those whose matrix is symmetric and positive definite. The method proceeds by generating vector sequences of iterates, residuals corresponding to iterates, and searching directions used in updating iterates and residuals. It is an algorithm for finding the nearest local minimum of a function of n variables. The idea of this algorithm is that the convergence to the solution could be accelerated if we minimize a function P over the hyperplane that contains all previous search directions, instead of minimizing P over just the line points down gradient. There are several attractive features of conjugated gradient method：First, under the premise that we regard numerical round-off, we can obtain acceptable accuracy of convergence within N iterations. Second, applying different preconditioning schemes, the rate of convergence can be improved greatly.
Third, users do not have to choose parameters experimentally.
Conjugated gradient method can be viewed as an improvement of steepest descent method [1]. Now suppose we want to solve the following system of linear
equations,
b
AX = , (27) where A is an n by n symmetric and positive definite matrix and A is real. We denote the unique solution of this system is . Without loss of generality,
is the initial guess for conjugate direction. Then forms a basis of
Dk Let demonstrates an inner product of U and V. To solve equation (27), we find the minimum value of this following quadric equation (29),
) This suggests that we can take the first basis vector to be the gradient of at , which means We note that is the negative gradient of , so the gradient decent method would take as the moving direction. But here we insist that the direction need to be conjugate to each other, so we choose the direction closest to under the conjugate constraint. Then it gives the following equation (31),
rk q( X) and the final resulting algorithm m is shown in Figure 3.3.
1
Figure 3.3: The algorithm of conjugate gradient method.
3.4 Ocean Surface Construction
Generally speaking, except using 3D-Navier Stokes equations to simulate ocean surface, most methods for constructing ocean surface are two dimensional methods using height map representing the variation of ocean surface. Though these 2D methods are hard to simulate the splashing or churning behavior of fluid in detail, they can reduce lots of computation in simulating large scale scenes.
We choose a method to produce fractal surfaces by using Perlin noise function [9]. With laying sets of different amplitudes and frequencies of noise function, we can
generate more fractal surface. Equation (32) demonstrates the output height of a fractal surface.
∑
= i perlin noise x xParameter α represents the amplitude. The higher value of α will give a rougher surface. In the other hand the lower value of α will give a smoother surface.
Then we can pre-generate three or more fractal surface height maps with different sets of parameters. We can synthesize these fractal surface maps into a final height map by an interpolation by dividing the phase of cosine function equally.
Through this interpolation method, the wave motion may rise and fall much more irregular in comparison with linear interpolation between fractal surface maps.
Equation (33) demonstrates the interpolation method for three fractal surfaces.
3
3.5 Rendering Ocean Surface
In order to render physically looking ocean surface, we introduce more physically optic reflection model. First of all, we assume that the ocean surface is perfect specular reflector and we know about the relation of reflection and transmission clearly. But under certain circumstances, we can not consider ocean surface as a perfect specular reflector. When the camera is at a large distance to the ocean surface, the reflected sunlight from the waves appears to be spread out and diffuse. Generally speaking, the behavior of light proceeding respected to a surface are split into two components. First, rays of light are reflected above the ocean surface.
Second, rays of light are transmitted though the surface. Here we will discuss how the two components work together.
3.5.1 Reflection
It is well known that in a perfect specular reflection, the incident ray and the reflected ray have the same angle respect to the surface normal. We may build up the following equation to express the reflected ray. In the beginning, let us define some notations: the wave height and horizontal position x, so the point r on the ocean surface in three-dimensional space can be expressed as follows,
) yv is a pointing-up unit vector. At the point r, the normal of the surface can be
calculated by its surface slope,
) then the surface normal is expressed as,
)
We denote the direction of incident ray of light is , according to the angle of incidence equals to the angle of reflection. We can obtain the reflected ray of light,
Ni
3.5.2 Transmission
The expression for transmission of light is not as simple as reflection of light.
We have to consider two things. First, the direction of transmitted ray is only related with the surface normal and the direction of incident ray. Second, the transmission between two different mediums obeys Snell’s Law.
Assuming that the direction of incident ray is in a medium with index of refraction , and the transmitted ray is in the medium with index of , so we can
Ni
ni nt
calculate the angle between incident ray and surface normal,
sinθi = 1−(Ni ⋅Ns)2 = Ni×Ns , (38) and the angle of the transmitted ray will be
s t
t = N ×N
θ
sin . (39) According to the Snell’s law, we can calculate the direction of refracted ray.
3.5.2 Fresnel Reflectivity and Transtivity
When light proceeds through media, parts of it will transmit from one medium to another. The other parts of it will be reflected. So theoretically no light is lost during the proceeding. In other words, reflectivity R and transtivity T are related by this following constraint,
=1 + T
R . (40) However the derivation of the expression for R and T is based on the electromagnetic theory of dielectrics. We do not attempt to introduce where it came from. We consult some charts in [10] which referenced some papers in surface wave optics, we directly use the following Fresnel term in our algorithm, which stands for R :
Chapter 4
Implementation and Results
We implement our system on a PC with 1.6GHz AMD Sempron Processor, 2.0 GB of system RAM, and a GeFroce 6200 PCI-X with 128MB of video memory. In average, it takes 0.29s to generate each frame. The conjugated gradient method solver for simulating the tsunami surface movement takes most of the computing time.
Though we only store non-zero parts of the n-by-n symmetric and positive-definite matrix in the conjugated gradient method, it indeed reduces the complexity from to ( n depends on the number of vertex). Solving conjugated gradient method is still time-consuming. Then we implement environment mapping in our shader program on GPU to make the surface reflection more physically real.
) (n2
O O(n)
In the following two sections, we will demonstrate how we simulate the propagation of tsunami wave in our system. In Section 4.1, we show the difference of water surface when different ocean topographies are applied for simulation. In Section 4.2, we show the simulations with different initial tsunami amplitudes. And the boat appearing in the scene is about 60 meters in length.
4.1 Different ocean topographies
Figure 4.1 is a height map which we download from a weather website [14]. It is an ocean topography in the direction of northeast of Taiwan. The brighter pixels mean the closer to the sea level.
Figure 4.1 The height map of the underwater terrain in the northeast of Taiwan In Figures 4.2 - 4.6, we only demonstrate the tsunami surface propagating on the ocean surface by setting Figure 4.1 as the ocean topography at different frames.
Figure 4.2 Frame No. 99 (Ocean256.raw)
Figure 4.3 Frame No. 199 (Ocean256.raw)
Figure 4.4 Frame No. 299 (Ocean256.raw)
Figure 4.5 Frame No. 399 (Ocean256.raw)
Figure 4.6 Frame No. 499 (Ocean256.raw) Figure 4.7 A height map edited by using image processing tool [13].
Figure 4.7 The edited height map
In Figures 4.8-4.12, we demonstrate the corresponding movement of surface.
Figure 4.8 Frame No. 99 (DEM256.raw)
Figure 4.9 Frame No. 199 (DEM256.raw)
Figure 4.10 Frame No. 299 (DEM256.raw)
Figure 4.11 Frame No. 399 (DEM256.raw)
Figure 4.12 Frame No. 499 (DEM256.raw)
Through these simulations, we can find out that the generated tsunami surface with real ocean topography is much rougher than the one with edited height map.
The reason is the depth of the ocean determines the velocity of wave. So the initial ocean wave propagating velocities on the surface with real ocean topography are changing. The depth varies so abruptly in the edited height map that the initial ocean wave propagating velocities differ a lot in that area. Then we can observe the surface generating much higher waves at the position.
The following Figures 4.13a-f and Figures 4.14a-f are the simulation of tsunami waves with applying ocean surface.
Figure 4.13a Frame No. 90 Figure 4.13b Frame No. 160
Figure 4.13c Frame No. 260 Figure 4.13d Frame No. 360
Figure 4.13e Frame No. 500 Figure 4.13f Frame No. 560
Figure 4.14a Frame No. 1 Figure 4.14b Frame No. 40
Figure 4.14c Frame No. 90 Figure 4.14d Frame No. 160
Figure 4.14e Frame No. 260 Figure 4.14e Frame No. 360
4.2 Different initial tsunami amplitudes
This section we describe the differences by using different initial amplitudes of tsunami waves.
Figures 4.15a-c show that the initial amplitude of tsunami wave is 15 meters, and Figures 4.16a-c show that the initial amplitude of tsunami wave is 1 meter. We can see that, the simulated tsunami surface is different apparently with different initial amplitude values. We also can notice that different amplitudes result in different wave velocities. So at the same frame, the wave crests are at different positions. Higher amplitude results in higher speed.
Figure 4.15a Frame No. 150 (amplitude 15m) Figure 4.16a Frame No. 150 (amplitude 1m)
Figure 4.15b Frame No. 300 Figure 4.16b Frame No. 300
Figure 4.15c Frame No. 540 Figure 4.16c Frame No. 540
Chapter 5
Conclusion and Future Works
We have presented a simulation system for tsunami wave propagation. By combining the shallow water equations which are often used in ocean engineering to simulate the movement of tsunami waves and the rendering techniques developed in computer graphics to render the 3D scene. We can obtain more physically-correct tsunami wave propagation. The shallow water equation is simplified by integrating the three-dimensional Navier-Stokes equations in z direction. It can simulate the wave movement correctly when the ratio of depth of water to the wave length is rather small, just like tsunami waves. Then we apply a more stable and efficient numerical method, called Semi-Lagrangian method, on SWE and it allows that we can maintain the stability of the simulation system with larger time steps. Conjugated gradient
We have presented a simulation system for tsunami wave propagation. By combining the shallow water equations which are often used in ocean engineering to simulate the movement of tsunami waves and the rendering techniques developed in computer graphics to render the 3D scene. We can obtain more physically-correct tsunami wave propagation. The shallow water equation is simplified by integrating the three-dimensional Navier-Stokes equations in z direction. It can simulate the wave movement correctly when the ratio of depth of water to the wave length is rather small, just like tsunami waves. Then we apply a more stable and efficient numerical method, called Semi-Lagrangian method, on SWE and it allows that we can maintain the stability of the simulation system with larger time steps. Conjugated gradient