Changeset 7549


Ignore:
Timestamp:
Oct 16, 2009, 12:04:30 PM (13 years ago)
Author:
jgriffin
Message:

corrected kinetic energy term

File:
1 edited

Legend:

Unmodified
Added
Removed
  • misc/tools/event_selection/wave_energy/wave_energy.py

    r7547 r7549  
    11"""
    2 Program to integrate the momentum quantities from a boundary time series.
    3 Used to develop different methods of estimating offshore tsunami hazard,
    4 rather than just wave height
     2Program to calcualte the energy density in offshore tsunami time series.
     3Energy is integrated over the time of the each wave crest.
     4Also calculates time integrated momentum quantities from the boundary time series.
     5Used to investiage different methods of estimating offshore tsunami hazard,
     6rather than just wave height.
    57
    6 Creator: Jonathan Griffin
     8Energy calculations from Johnson 1997. 'A modern introduction to the mathematical
     9theory of water waves' p23.
     10
     11Creator: Jonathan Griffin, Geoscience Australia
    712Created: 12 October 2009
    813"""
     
    1318import csv
    1419import glob
    15 ##import matplotlib
    16 ##matplotlib.use('Agg')
    17 ##import pylab
    18 ##from pylab import *
    1920import numpy
    2021
     22# Change path here!!!
    2123filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries'
    2224index = 2872
     
    2729filenamelist = []
    2830
     31
     32# Find all relevant gauges files and put into list
    2933for root, dirs, files in os.walk(filepath):
    30 ##    print 'root', root
    31 ##    print 'dirs', dirs
    32 ##    print 'files', files
     34
    3335    for filenames in files:
    3436        if filenames == filebasename:
     
    3638            filepathlist.append(filepathname)
    3739            filenamelist.append(filenames)
     40           
     41# Initialise lists for storing maximum values   
     42max_energy_list = []
     43max_energy_inst_list = []
     44counter = 0
    3845
    39 max_energy_list = []
    40 max_momentum_list = []
    41 max_energy_inst_list = []
    42 
    43 counter = 0
     46###########################################################
     47# Loop over gauges files
     48###########################################################   
    4449for filename in filepathlist:
    4550    time = []
     
    5762        y_momentum.append(float(line[3]))
    5863        # Calculate momentum norm
    59         momentum.append(numpy.sqrt(numpy.power(float(line[2]),2)+numpy.power(float(line[3]),2)))
    60 
     64       
    6165    time_diff = time[1]-time[0]
    6266
    6367    sign = 0
    64     max_momentum = 0
    65     momentum_sum = 0
    6668    energy_sum = 0
    6769    max_energy = 0
    6870    max_energy_inst = 0
     71   
     72    ###############################################################
     73    # Loop through time series and add values (integrate over time)
     74    # Sum set to zero when stage height changes sign
     75    ###############################################################
     76    for i in range(len(time)):
    6977
    70     for i in range(len(time)):
    71        
    72         # Kinetic energy
    73         KE = (1./2.)*(numpy.power(x_momentum[i],2)+numpy.power(y_momentum[i],2))
    74         # Potential energy
    75         PE = 9.81*(stage[i])
     78        ###############################################################
     79        # Energy calculations - note that quantities from gauge file are
     80        # already depth integrated.
     81        ###############################################################
     82        # Calculate velocities from 'momentum' terms (Note not true momentum,
     83        # rather defined as u*h. So we devide by the depth to get u).
     84        # roh = density = 1000 kg.m-3
     85        x_vel = x_momentum[i]/(100+stage[i])
     86        y_vel = y_momentum[i]/(100+stage[i])
     87        # Kinetic energy KE = 1/2 roh h u^2
     88        KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i])
     89        # Potential energy PE = roh g w (w = stage)
     90        PE = 9.81*(stage[i])*1000
    7691        energy = KE+PE
    7792
     93        # Maximum instantaneous energy density
    7894        if energy > max_energy_inst:
    7995            max_energy_inst = energy
     96
     97        ###############################################################
     98        # If sign of wave changes (goes from trough to crest or vice-versa)
     99        # multiply energy sum by timestep and reset sum value to zero
     100        ###############################################################
    80101       
    81     ##    print 'max_energy',max_energy
    82     ##    print 'energy_sum', energy_sum
    83102        if stage[i] > 0 and sign != 1:
    84103            sign = 1
    85             temp_max_momentum = momentum_sum*time_diff
    86104            temp_max_energy = energy_sum*time_diff
    87             if temp_max_momentum > max_momentum:
    88                 max_momentum = temp_max_momentum
    89105            if temp_max_energy > max_energy:
    90106                max_energy = temp_max_energy
     107               
    91108            energy_sum = energy
    92             momentum_sum = momentum[i]
    93 
    94109            start_time = time[i]
    95110           
    96111        elif stage[i] < 0 and sign != -1:
    97112            sign = -1
    98             temp_max_momentum = momentum_sum*time_diff
    99113            temp_max_energy = energy_sum*time_diff
    100             if temp_max_momentum > max_momentum:
    101                 max_momentum = temp_max_momentum
    102114            if temp_max_energy > max_energy:
    103115                max_energy = temp_max_energy
     116               
    104117            energy_sum = energy
    105             momentum_sum = momentum[i]
    106            # print momentum_sum
    107118            start_time = time[i]
    108119           
    109120        elif stage[i] == 0 and sign != 0:
    110121            sign = 0
    111             momentum_sum = momentum[i]
    112122            energy_sum = energy
    113 
    114123            start_time = time[i]
    115124           
    116125        else:
    117             momentum_sum += momentum[i]
    118126            energy_sum += energy
    119 
    120     ##    print 'KE',KE
    121     ##    print 'PE',PE
    122127           
     128    # Add maximum values to list       
    123129    max_energy_inst_list.append(max_energy_inst)
    124130    max_energy_list.append( max_energy)
    125     max_momentum_list.append(max_momentum)
    126 ##    print 'max_momentum', max_momentum
    127 ##    print 'max_energy', max_energy
     131
    128132    counter += 1
     133
     134   
    129135counter = 0
    130136total_max_energy = 0
    131 total_max_momentum = 0
    132137total_max_inst_energy = 0
    133138
     139# Print results for each file and calculate maximum value across all events
    134140for filename in filepathlist:
    135141    print filename
    136     print 'max energy', max_energy_list[counter]
    137     print 'max instantatneous energy', max_energy_inst_list[counter]
    138     print 'max momentum', max_momentum_list[counter]
     142    print 'max crest-integrated energy:', max_energy_list[counter],'J.s/m^2'
     143    print 'max instantatneous energy:', max_energy_inst_list[counter],'J/m^2'
    139144
    140145    if max_energy_list[counter] > total_max_energy:
     
    144149        total_max_inst_energy = max_energy_inst_list[counter]
    145150        max_energy_inst_index = counter
    146     if max_momentum_list[counter] > total_max_momentum:
    147         total_max_momentum = max_momentum_list[counter]
    148         max_momentum_index = counter
     151
    149152    counter+=1
    150153
    151 print '\n File with maximum energy integrated over a wave crest is: ', filepathlist[max_energy_index]
    152 print 'Energy is ', max_energy_list[max_energy_index], ' J.s/m^2'
    153 print '\n File with maximum instantaneous energy is: ', filepathlist[max_energy_inst_index]
    154 print 'Energy is ', max_energy_inst_list[max_energy_inst_index], ' J/m^2'
    155 print '\n File with maximum momentum is: ', filepathlist[max_momentum_index]
    156 print 'Momentum is ', max_energy_list[max_momentum_index], 'kg.m/s'
     154print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index]
     155print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2'
     156print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index]
     157print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
Note: See TracChangeset for help on using the changeset viewer.