# Nicholson Canyon Simulation

### From Gerris

Revision as of 05:19, 9 August 2011EmilyMLane (Talk | contribs) (→Velocity of Sliding Block) ← Previous diff |
Current revisionEmilyMLane (Talk | contribs) (→Velocity of Sliding Block) |
||

Line 17: |
Line 17: | ||

In our simulations | In our simulations | ||

\theta=9.09^\circ = \alpha_1 b=1.7km w=1,7km l=0.272 km | \theta=9.09^\circ = \alpha_1 b=1.7km w=1,7km l=0.272 km | ||

- | In the example Watts et al 2000 used they took C_n=0, C_d=1.7 and C_n=0.8, there does not seem to be much in the literature about the correct values to use for these variables in realistic geophysical scale problems | + | We take the density of the block to be 2 g/cm^3 = 2000 kg/m^3, thus m_b = \rho_b V. In the example Watts et al 2000 used they took C_n=0, C_d=1.7 and C_n=0.8, there does not seem to be much in the literature about the correct values to use for these variables in realistic geophysical scale problems. |

+ | |||

+ | Using these values gives | ||

+ | a_0=0.5535m/s^2 | ||

+ | u_t=39.6 m/s | ||

+ | |||

+ | Because Gerris cannot handle the sliding block impinging on the edge of the domain we need to ensure that the block stops before it reaches the bottom of the canyon. We do this by picking a point where the block begins decelerating and calculate the linear deceleration necessary for the block to stop before it impinges on the far side of the canyon. By playing around with the distance in which the block stops (and thus the deceleration rate) we can also see how this affects the size of the wave generated. | ||

== Two-Dimensional Simulations == | == Two-Dimensional Simulations == |

## Current revision

## Contents |

## Problem Definition

These simulations are loosely based on the geometry of the Nicholson Canyon in Cook Strait. It assumes a simple sliding block with a prescribed velocity. This work is similar to that of the benchmark cases of Enet *et al.* and Liu *et al.* with the difference being that while the two benchmark cases involved laboratory scales this simulation involves realistic geophysical scales.

### Velocity of Sliding Block

Watts et al 2000 uses <math> a_0 = \frac{g(m_b-m_0)(\sin{\theta} - C_n \cos{\theta})}{m_b+C_m m_0} </math> and <math> u_t= \sqrt{\frac{2g(m_b-m_0)(\sin{\theta} - C_n \cos{\theta})}{C_d \rho_0 w b \sin{\theta}} </math> In our simulations \theta=9.09^\circ = \alpha_1 b=1.7km w=1,7km l=0.272 km We take the density of the block to be 2 g/cm^3 = 2000 kg/m^3, thus m_b = \rho_b V. In the example Watts et al 2000 used they took C_n=0, C_d=1.7 and C_n=0.8, there does not seem to be much in the literature about the correct values to use for these variables in realistic geophysical scale problems.

Using these values gives a_0=0.5535m/s^2 u_t=39.6 m/s

Because Gerris cannot handle the sliding block impinging on the edge of the domain we need to ensure that the block stops before it reaches the bottom of the canyon. We do this by picking a point where the block begins decelerating and calculate the linear deceleration necessary for the block to stop before it impinges on the far side of the canyon. By playing around with the distance in which the block stops (and thus the deceleration rate) we can also see how this affects the size of the wave generated.

## Two-Dimensional Simulations

These simulations are run as a 2D vertical slice (i.e. assuming that the canyon and landslide are infinitely long and that there is no movement in the along-canyon direction

### Generating the Landslide Body Geometry

First we need to create a .gts object that is the sliding block. One difficulty is that we have to ensure that the block is very close to the edges of the domain while still all lying inside the domain (otherwise those points lying outside the domain get 'left behind'). In the two dimensional simulation we still need a three dimensional block to ensure that it properly intersects the domain. This is made by first creating points on the surface of the block using an .awk file and then triangulating them using the gts delaunay triangulation

# Generating points using .awk file

awk -f nicholson.awk > points.tmp

# Delaunay triangulation and transformation using gts tools

( wc points.tmp | awk '{print $1" 0 0"}'; cat points.tmp ) | gtsdelaunay -d | gtstransform --rx=-90 > nichol_test.gts

# Shifting points to the domain boundary

sed -e "s/ -0.0168/0/g" < nichol_test.gts > nichol_testout.gts

The .awk file used to generate the points is

BEGIN {

a=0.272

b=1.7

alpha=atan2(a,b);

h=a*cos(alpha);

l1=b*cos(alpha);

l2=a*sin(alpha);

for (k=0; k <= 100; k++) {

for (i=-1; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -

0.85+1.7*k/100 , h*i/100) }

for (i=1; i <= 101; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100,

-0.85+1.7*k/100, h*(1-i/100)) }

}

dx=1.0e-4

k=0.

j=-1.

for (i=-1; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -0.

85+1.7*k/100-dx,h*j/100) }

for (i=1; i <= 101; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100, -

0.85+1.7*k/100-dx,h*j/100) }

k=100.

for (i=-1; i <= 100; i++) { printf("%5.4f %5.4f %5.4f\n",l1*i/100, -0.

85+1.7*k/100+dx,h*j/100) }

for (i=1; i <= 101; i++) { printf("%5.4f %5.4f %5.4f\n",l1+l2*i/100, -

0.85+1.7*k/100+dx,h*j/100) }

}

### Numerical "wave gauge"

The height of the water is measured using the same technique as the Liu_et_al_Benchmark. Probes are placed two and one kilometre before the edge of the canyon, at the edge of the canyon, halfway along the sliding block, at the original edge of the sliding block, above the deepest point the the canyon, at the far edge of the canyon and one and two kilometres beyond the far edge of the canyon. These are labelled probe1 through probe9 respectively.

The .awk file used to generate the probe positions (position.awk) is

BEGIN {

angle = atan2(4,25)

for (y = -0.5; y <= 0.5; y+= 0.001) {

x = -2/cos(angle)

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe1x"

x = -1/cos(angle)

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe2x"

x = 0

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe3x"

x = 0.8608

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe4x"

x = 1.7216

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe5x"

x = 2.5318

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe6x"

x = 3.9496

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe7x"

x = 3.9496+1/cos(angle)

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe8x"

x = 3.9496+2/cos(angle)

print cos(angle)*x - sin(angle)*y, sin(angle)*x + cos(angle)*y, 0 > "probe9x"

}

}

The file that calculates the water level from the raw output (distance.awk) is

BEGIN {

FS=" |:"

sum=0

}

{

if (NR > 2 && ( t0 != $1 || Y0 > $3 ) ) {

print $1,0.5-sum;

fflush (stdout);

sum=0;

}

if ($1 == "#") {

# get column indices of the relevant fields (X,Y)

for (i = 1; i <= NF; i++) {

if ($i == "T")

iT = $(i-1);

}

}

else {

if ($1 == t0 && Y0 <= $3) {

T = $iT

sum = sum + 0.001*(T+T0)/2.;

}

}

T0=$iT; X0=$2; Y0=$3; t0=$1

}

Note that this is very similar to that used in the Liu_et_al_benchmark differing only in small details (lines 6 and 7).