
















Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Material Type: Project; Class: Fnd of Computatnl Sci; Subject: Computational Sci& Informatics; University: George Mason University; Term: Unknown 1989;
Typology: Study Guides, Projects, Research
1 / 24
This page cannot be seen from the preview
Don't miss anything!

















CSI 701 Final Paper
Structural damage study for Tsunami waves
Deepthi Achuthavarier Kristi Arsenault Ragini Dravida Musaddeque Hossein Dan Nagle
Submitted on: 05/02/
tsunami's energy flux, which is dependent on both its wave speed and wave height, remains nearly constant. Consequently, as the tsunami's speed diminishes as it travels into shallower water, its height grows. Because of this shoaling effect, a tsunami, imperceptible at sea, may grow to be several meters or more in height near the coast. When it finally reaches the coast, a tsunami may appear as a rapidly rising or falling tide, a series of breaking waves, or even a bore.
What happens when a tsunami encounters land?
As a tsunami approaches shore, it begins to slow and grow in height. Just like other water waves, tsunamis begin to lose energy as they rush onshore - part of the wave energy is reflected offshore, while the shoreward-propagating wave energy is dissipated through bottom friction and turbulence. Despite these losses, tsunamis still reach the coast with tremendous amounts of energy. Capable of inundating, or flooding, hundreds of meters inland past the typical high-water level, the fast-moving water associated with the inundating tsunami can crush homes and other coastal structures. Tsunamis may reach a maximum vertical height onshore above sea level, often called a run-up height, of 10, 20, and even 30 meters.
Numerical Modeling of Tsunami waves
. There are three distinct phases in the life of a tsunami. These are - tsunami generation, propagation and run-up. Numerical modeling of a tsunami involves simulating these three phases. Wave generation is often by an underwater earthquake or landslide. The sea level rise created as a result of such phenomena can be computed by simple analytical formulae.
Wave propagation is the energy propagation in deep ocean. The speed of tsunami waves
depends on the depth of the ocean ( c = gH , where g =acceleration due to gravity and
H =depth of ocean), hence waves undergo accelerations and decelerations in passing over an ocean bottom of varying depth. In deep ocean, tsunami waves travel at speeds of 500 to 1000 km/hr and their wavelengths range between 500 and 650 km.
Wave run-up in tsunami modeling involves the simulation of wave height as it approaches shallow waters. The National Oceanic and Atmospheric Administration (NOAA) defines wave run-up as the difference between the elevation of maximum tsunami penetration (inundation line) and the sea-level at the time of tsunami attack.
Figure 2 : Tsunami wave run-up and terminology associated with it Source: NOAA, www.noaa.gov
1.1 Objectives of this study
In this study, we first attempt to simulate a one-dimensional Tsunami wave propagation in deep ocean using linear shallow water theory and a computationally efficient numerical algorithm. Consequently, we not only try to simulate a wave run-up by a similar model but also include bottom friction and non-linear advection terms. Finally, we attempt to quantify the impact force exerted onto a coastal structure and the deflection thereby produced in the event of a Tsunami.
2 Methods
In this section, the theory used in modeling tsunami wave propagation, run-up and deflection of a coastal structure are described. The numerical methods used in each simulation are also explained.
2.1 Wave propagation
Tsunami waves in deep ocean can often be approximated as horizontally propagating linear shallow water gravity waves (Holton, 1992). The term shallow water arises from the fact that the wavelength is much greater than the depth of the medium for these types of waves. Shallow water waves can exist only if the fluid has either a free surface or an internal density discontinuity; the former is satisfied in the case of a tsunami wave. In shallow water waves, the restoring force is gravity and therefore the individual water particle oscillations are perpendicular to the propagation of the wave. Assuming that the motion is two dimensional in the x , z plane, equation of motion along the x direction can be written in a simplified form as shown below.
x
g t
u ∂
∂ η (1)
Figure: 3. Space staggered grid for the wave propagation problem (adapted from Kowalik and Murty, 1993)
The space index runs from j = jleft to j = jright. The velocity computations start at j = jleft +1 and proceeds to the point j = jright. Sea level computation starts at the left boundary ( j = jleft ) and ends at one grid point before the right boundary ( j = jright - 1 ).
At the left end of the channel a tsunami wave is generated as a sinusoidal signal.
= t T
Sin p
wave towards the right end of the channel is modeled. Once the wave is completely formed which happens when t = Tp , an open boundary condition is used at the left boundary. At the right boundary, open boundary condition is used throughout the simulation period. Open boundary condition is formulated in such a way that the signal is allowed to propagate outward without any reflection. Kowalik and Murty (1993) implemented a practical approach to the open boundary condition for shallow water waves, the basic idea of which stems from the radiating boundary condition proposed by Reid and Bodine (1968). According to this, for a one-dimensional shallow water wave, the velocity at the left boundary can be specified as follows.
S V S V S V S V S V S V S V
Velocity grid point Right boundary ( jright )
Sea level grid point Left boundary ( jleft )
g
Similarly, at the right boundary, velocity is computed using the sea level height value.
g
The space and time steps are chosen to satisfy the Courant-Friedrichs-Lewy condition which is given in equation (9).
h gH ≤ (9)
2.3 Wave run-up
To represent the run-up of the tsunami wave in a finer scale modeled approach, we use a one-dimensional model (varying in the x-direction and in time) originally developed by Kowalik and Murty (1993). The model calculates the run-up by first initiating the wave as a defined sinusoidal curve, similar to the wave set up for tsunami propagation program used also for this project. The governing equations for the model involve equations of motion and continuity along with the nonlinear advective term. Following the shallow water theory we constructed a simple tsunami wave run-up model. Horizontal momentum equation now incorporates advection and friction.
ruu x
g x
u u t
u − ∂
Here r is the bottom friction coefficient and all other terms are defined in previous sections. Equation of continuity used is the same as in equation 3.
2.4 Numerical scheme
The wave run-up problem is modeled using the following numerical scheme. The non- linear advective terms as well as the equation of continuity are approximated using the upwind/downwind scheme.
m j
m j m j
m j
m j
m j
m j
m j
m j
m j D D
ru u h
g h
u u un h
u u up T
u u
1
1 1 1
1
− + −
Equation of continuity is modeled as
K = the stiffness matrix of the system u&& (^) = the acceleration of the system
u = the displacement of the system
Equation 13 is similar in appearance to equation of motion of a single degree of freedom system without and damping (i.e. without a “dashpot”), with one important difference however. The M and K are matrices representing the mass and stiffness of the structure. For a multi story frame analysis, the mass of each floor (windows, doors, walls and fixtures, and all other live loads) is assumed to be concentrated at the center of each floor slab. The columns, besides functioning as floor framing supports are also assumed to provide the restoring force countering a floor deflection. Therefore the stiffness of the system is assigned to the columns. Therefore equation 13 represents N simultaneous equations with time derivatives of the displacement vector u. Formulation of the mass and the stiffness in this manner is referred to as the “lumped mass model”.
The formulation of the mass and stiffness matrices require a definition of a coordinate system. In accordance with the “lumped mass model”, displacement of the floor mass centered in floor slab takes place at the floor level only. Therefore each degree of freedom represents a possible lateral translation at the floor level. Hence, the number of degrees of freedom of the system is the same as the number of mass points of the system. This simplifies the matrix formulation. The mass matrix is then a diagonal matrix containing the mass of each floor as a diagonal element, and the stiffness matrix represents the stiffness of each mass of the system. Both matrices are square, with the mass matrix being diagonal as mentioned previously and the stiffness matrix having diagonal and non diagonal terms.
The lumped-mass method of mass and stiffness formulation results in equation 13 being a system of coupled equations also, since the stiffness of the system is interdependent on the stiffness of each of the stiffness elements. Application of normal coordinate transformation to the mass and stiffness matrices results in uncoupled equations, with each uncoupled equation representing a single-degree of freedom system.
The transformation matrix, also called the “mode shape” matrix is a matrix of the eigen vectors of the system in free vibration without any damping. The mode shape matrix represents the possible lateral displacement of system in free vibration. Therefore, for a mass matrix M, stiffness matrix K , the mode shape matrix, Φ is obtained by solving the eigen value problem:
M x = λ K x (14)
where, x = the eigen vectors of the system in free vibration λ = the eigen value of the system in free vibration
On applying the transformation as discussed above, the following transformed equation for the system is obtained:
M y&& +Ky=P * (15)
where, M* = ΦT^ M Φ C* = ΦT^ C Φ P* = ΦT^ P Φ y&&^ and y = acceleration and displacement respectively, associated with the transformed
equations.
As a result of the above transformation, the mass, stiffness and force matrix become diagonalized and equation of motion of the nth^ degree of freedom becomes:
M (^) n &y& +Kny=pn (16)
where P (^) n represents the component of the total force acting on the nth^ degree of freedom. Since all these parameters can be computed or measured, the system of uncoupled equations can be solved for the displacement using a computer. Once the displacement has been determined, it has to be transformed into the original system:
u = Φ y (17)
When the force varies with time, the system of equation will yield a displacement for each force increment. Consequently, the displacement at each degree of freedom will be the sum of each incremental displacement transformed to the original system:
n
j
n 1
= Φ (18)
Where an exact solution is not available, numerical methods must be employed to compute the displacement vector, y. The numerical schemes used in the project code to determine the displacement vector is discussed in section 2.6.
2.5.2 Computation of Tsunami Wave Force
The program generates the wave amplitude with computed coordinate points in two- dimension along the length of the wave. However since the wave shape is not a known function, the force of the wave acting on the structure will have to be computed along discrete points along the wave.
Since all points are referenced from the same elevation, the coordinates can be used to divide the approaching wave into rectangular slices along the length of the wave for force evaluation. If we assume that the length and height of each wave slice is d and h respectively, then the mass of each slice is
Wm = dhρ (19)
(y y hy ) h
4 y&& (^) n + 1 = −&y&n+ 2 n+ 1 − n− &n (23)
(y y) h
2 y& (^) n + 1 = −yn+ n+ 1 − n (24)
The equation of motion for the multi-degree of freedom system can also be modified accordingly:
My& &^ n + 1 +Kyn+ 1 =pn+ (^1) (25)
Displacement between time n and n+1 can be determined from the following:
y y ) h
4 y h
4 M K)y p M( h
4 ( 2 + (^) n+ 1 = n+ 1 + 2 n+ & (^) n+&&n (26)
Since the structure is assumed to be at rest prior to impact of the first wave slice, initial velocity and displacement is assumed to be zero. Upon impact of the first wave slice, the acceleration at each degree of freedom can be computed from equation 22. With the initial conditions determined, the displacement is determined from equation 26 and applied to equations 23 and 24 to determine the acceleration and velocity at time interval n+1. Using this methodology, the reaction of the structure as it is impacted by each wave slice can be determined. The time step “h” is same as that used for computing the wave run up.
The deflection y determined from equation 26 at each time interval at each degree of freedom can be summed and transformed into the deflection vector u at each degree of freedom in original coordinate system using equation 18.
Once the deflections have been determined, the shear force s in the columns of stiffness k is
s = u k (27)
This method can be implemented using an iterative loop to determine the displacement vector u over the entire wave.
3. Results
3.1. Wave propagation
In the first simulation in this section, we used the linear shallow water model developed in section 2.2 to propagate a tsunami wave along a channel of 1000km length. The depth of the channel is set as 4000m, which is the average depth of Pacific Ocean. Initial amplitude of the wave is given as 0.3 m and wave period 10 min. Using Fischer’s
algorithm, numerical integration is carried out for a space step of 2km and time step of 5sec. Figure 4 shows results of this simulation as snapshots of the wave as it propagates from the left end of the channel. We find the wave retaining its initial amplitude throughout the propagation, as no bottom friction dissipation is included in the model. The wave takes about one and a half hours to reach the right end of the channel with an approximate velocity of 769km/hr. This value is generally in agreement with the observed tsunami speed in deep ocean, which ranges from 300 km/hr to 900km/hr.
In the second simulation, we incorporated a sloping ocean bottom in the propagation model. The idea of this exercise was to understand how well a simple linear model reproduces tsunami wave amplification. In this simulation, the length of the channel is taken as 500km and the ocean bottom is sloping at angle of 0.74 degrees from 200km from the left boundary. The ocean depth is decreasing from 4000m to 150m (figure 5). Figure 6 shows results of this simulation as snapshots of wave as it propagates along an up sloping channel. We find the velocity of the wave is approximately 500km/hr as it travels along the sloping bottom. The velocity of the wave is found to be less than that obtained from the previous simulation which is in agreement with the general tsunami behavior. We also found that the wave height increases from 1m to 2m and wavelength decreases from about 120km to 40km as the wave travels upslope. A tsunami wave simulation found on the web page, http://www.math.uio.no/avdb/gitec/planok/ produces similar results and shows that the wave height doubles as it moves upslope.
In the third simulation we incorporated the advective term in horizontal momentum equation to understand the influence of this term in wave propagation. We ran the model for the same initial and boundary conditions but one with and the other without horizontal advection. We found that in deep ocean propagation, the simple linear model is able to produce results as accurate as the non-linear model (figure 8). This helps us conclude that in deep ocean propagation, non-linear terms can be omitted safely.
3.2 Wave run-up
The original run-up model was updated and modified for this project, though the original parameters were used for code testing and for the testing of the structural model component. Figure 8 below shows the amplitude of an idealized tsunami wave (represented by the red line) as a function of distance from the left boundary as it runs up (stage 1-3) and runs down (stage 4).These results are generated by the updated version original tsunami run-up channel model (Kowalik and Murty, 1993) and similar to their results confirming the proper working of the code. The modified code was used as a control simulation, using the original parameters as designated by Kowalik in his original version. As indicated in the figure 8, this simulated tsunami runs up the land to a height of more than 4 meters before retreating. The model actually produces a series of waves, a realistic representation of typical tsunami events, but only a single wave is shown here.
structure deflections were determined from the output of the validated tsunami code, and the inter-story shears were computed. The code provides the maximum deflection and forces produced in a stand alone structure that is absorbing the complete impact force of a tsunami run-up wave. However, this simplification is not realistic as it does not take into account structure wave interaction between other structures that may be adjacent to the structure being examined such as other buildings, wave breakers and seawalls, and other costal features that may reduce the magnitude of the impact. Study and incorporation of any influence on the impact force due wave-structure interaction between the runup and structures adjacent to the structure of interest was beyond the scope this project. Additionally, the wave is assumed to transfer all its energy to the structure. However, this is not the case in reality. Upon impact, a run-up wave becomes a breaking wave. Part of the breaking wave transmits its energy to the impacted structure. The remaining part of the breaking wave interferes with the subsequent run-up waves, reducing the impact somewhat. Incorporation of influence due to wave interference on the impact force was also beyond the scope of the project.
Furthermore, the plan shape of the structure has not been taken into account. The plan shape of the structure will influence where and how the wave impacts the structure. This is influence the overall deflection and the inter-storey forces. Normally, accepted building codes provide coefficients to accommodate for any plan irregularities with which the design force is modified to account for the influence of the plan irregularity on the impact force. Development of such coefficients was beyond the scope of this project.
Based on the above discussion, suggested future work is as follows:
improve readability and ease of use by other users. The main program was converted into a subroutine (now called runup_channel.f90) and is now called by a driver program (runup_driver.f90) that reads in user inputs for the model from input file named, tsunami.crd. The tsunami case study used in this project utilized real surveyed measurement values to set up and run simulations with. To generate ocean bottom conditions for the models, a Fortran program was developed to read in a topographic dataset for a basin domain. The program then locates the depths of where an earthquake or tsunami-originating point and the location of where the tsunami ran up onto a land point, both specified by the user. From these depths, the slope and distance between the points are calculated and provided as the output. A text input file is provided for the user to enter the latitude and longitude for both the tsunami rupture and runup geographic locations. The code was tested out on a real set of subsetted bathymetric datasets and validated against topographic maps
6. Citations - Holton, J.R., (1992), an introduction to dynamic meteorology, Academic Press, San Diego. - Kowalik, Z and Murty, T.S (1993), Numerical modeling of ocean dynamics, World scientific, Singapore. - Mercado, A. and W. McCann, 1998. Numerical simulation of the 1918 Puerto Rico tsunami, J. Natural Hazards, Vol. 18, 57-76. - Sielecki, A and Wurtele, M.G.,1970. The numerical integration of non linear shallow water equations with sloping boundaries. J. Computational Physics, 6, 219- - http://www.geophys.washington.edu/tsunami/general/physics/physics.html - Smith, W.H.F. and D.T. Sandwell, 1997, Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings, Science 277 (5334), p.1956-1962. - Ghali, A and M. Neville Structural Analysis: A Unified classical and Matrix approach. - Humar, J.L., Dynamics of structures
Figure 5 : Ocean depth as a function of distance
Figure 6 : Propagation of tsunami wave along an up sloping channel
0 100 200 300 400 500
0
2
0
2
0
2
0
2
0
2
0
2
t = 10min
amplitude of wave (m)
distance from offshore boundary (km)
t = 20min
t = 30min
t = 40mins
t = 50min
t = 60min
Figure 7 : variation of maximum wave height with slope angle
Figure 8 : Wave propagation with and without advection term for two points tin time
400 420 440 460 480 500 520
-2.
-1.
-1.
-0.
Wave propagation with and without advective terms
t = 60 min
with advection without advection
amplitude of wave( (m)
distance from offshore boundary (km)
1.0 1.2 1.4 1.6 1.8 2.0 2.
Increase in wave amplitude with slope
maximum wave amplitude (m)
sea bottom slope (degree)