|
|
|
|
Team Number: |
098 |
|
School Name: |
Santa Fe Preparatory School |
|
Area of Science: |
Geophysics |
|
Project Title: |
Fire on the Mountain |
|
Project Abstract: |
|
|
Interim Report: |
|
|
Final Report: |
http://mode.lanl.k12.nm.us/97.98/finalreports/098/finalreport.html |
EXECUTIVE SUMMARY
The purpose of our project was to create a finite difference program that would predict the areas of devastation from a volcanic eruption. We are using basic conservation equations of mass, momentum, and later energy to figure out the location of these areas and the magnitude of effect. The program will be based on a three dimensional matrix of finite difference objects encapsulating these three conservation equations. Looking at time intervals, we are going to display our results in a comprehensive manor.
When considering actual sites variables such as the viscosity of the flow, slope topography, vegetation, temperature, and the wind will have to be taken into consideration. The dynamics of the eruption is also a consideration in our predictions. Though for simplification’s sake and because of the time limit some of these more minor variables will have to be eliminated from our program.
The program will be achieved computationally by iterating a process in which every element is altered by interacting with its adjacent cells in accordance with proven techniques of numerical analysis. The iterations will represent the behavior of the flow in brief intervals of time. When repeated by a super-computer the simulation will be a mathematically accurate approximation.
INTRODUCTION
The devastation that a volcano can create is important to predict not only as a scientific multi-phase problem but also because predictions could save lives, property, and other mean necessary for human survival. There are many variables in a volcanic eruption wind, caldera size and shape, temperature of the volcanic material, the compositions of the material, the viscosity, so on and so forth. Not all these variables could be accounted for in our meager program, so we chose the ones that were prevalent.
We used a finite difference scheme to help us go from partial difference equations to finite difference equations. This way we can look at cubic cells instead of every individual point in our defined (x, y, z) mesh. The place of the eruption is defined with two boundaries. The first most obvious one is where our graphical representation ends. On the top and sides the cells along this boundary flow over. They material and velocity going in the positive z and x direction as well as the matter moving in the negative x direction pours out of the mesh that we define to have a finite height, width, and depth. The other boundary condition is the caldera itself. Along the solid sides of the cone the falling and spilling mater that comes from the vent faces a curve of cells whose input can only be zero. Thus the mater comes to rest and piles up.
* Out of the vent is where the material receives its initial condition. For now the conditions we have programmed in are density and velocity. Later, we will add energy. Using explicit differentiation both of these physical laws have been broken down into (x, y, z) coordinates where mass is conserved and vectors of velocity. Gravity and the source information are taken into account with these equations.
This is a simplified form of an eruption. Temperature, forces outside of the cone aside from gravity, the terrain of the cone, and other typical world conditions are not taken into consideration. However, as I stated before we only have so much knowledge and there is only so much time. We hope that the decisions we have made about which equations to leave in and which to leave out and the ones left in will prove prudent.
CODE
//Program: idealfluid.cpp
//Author: Devin Green
//Date: 11/30/1997
//Description: flow program. 2d for now
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <math.h>
#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <except.h>
void initializeg () ;
int main () {
volatile int c, i, j, x;//used for referencing
//volatile int k, s, a, b; //LIC
const int n = 40; //array size
float stepXY, stepT, alpha;
int topography[n], noise[n][n];//noise is for LIC graphics
double **p, **u, **v, **p2, **u2, **v2;
try {
p = new double*[n];
u = new double*[n];
v = new double*[n];// w[n][n],
p2 = new double*[n];
u2 = new double*[n];
v2 = new double*[n];//, w2[n][n];
for (c = 0; c < n; c++) p[c] = new double[n];
for (c = 0; c < n; c++) u[c] = new double[n];
for (c = 0; c < n; c++) v[c] = new double[n];
for (c = 0; c < n; c++) p2[c] = new double[n];
for (c = 0; c < n; c++) u2[c] = new double[n];
for (c = 0; c < n; c++) v2[c] = new double[n];
//density and velocity arrays
}
catch (xalloc) { // ENTER THIS BLOCK ONLY IF xalloc IS THROWN.
// YOU COULD REQUEST OTHER ACTIONS BEFORE TERMINATING
cout << "Could not allocate. Bye ...";
exit(-1);
}
initializeg();
cout << "Enter the XY step in meters: "; cin >> stepXY;
cout << "Enter the time step in seconds: "; cin >> stepT;
alpha = stepT / 2 / stepXY;
const float g = 9.8 * stepT; // use SI units
for (x = 1; x <= 15; x++) {
setrgbpalette(x, x*4, 0, 0);
setpalette(x, x);
}
setrgbpalette(0, 0, 0, 50);
setpalette(0,0);
setrgbpalette(1, 0, 50, 0);
setpalette(1,1);
randomize();
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
noise[i][j] = random(15) + 1;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++) {
p[i][j] = 0; p2[i][j] = 0;
u[i][j] = 0; u2[i][j] = 0;
v[i][j] = 0; v2[i][j] = 0;
}
// for (x = 0; x < n; x++)
// topography[x] = 0;
// for (x = 0; x < n; x++)
// topography[x] = n/2 - abs(x - n/2); //10 + floor(10*sin(x*3.14/10));
for (x = 0; x < n; x++)
topography[x] = 3;
for (i = 15; i < 25; i++)
for (j = 4; j < 14; j++)
p[i][j] = 1;
for (x = 0; x < n; x++)
line(x, n, x, n-topography[x]);
for (i = 0; i < n; i++)
for (j = 4; j < n; j++) {
if (p[i][j] == 0) c = 0;
else { /* this will be LIC
k = sqrt(pow(u[n*i + j], 2) + pow(v[n*i + j], 2));
a = i; b = j; s = 0;
for (c = 0; c < k; c++) {
s = noise[a][b];
Yet to do here
}
c = s/k; */
c = p[i][j];
c = c % 14 +2;
}
putpixel(i, n-j, c);
}
getch();
for (x = 0; !kbhit(); x++) {
for (i = 1; i < n-1; i++)
for (j = topography[i] + 1; j < n-1; j++) {
p2[i][j] = p[i][j] - alpha * ( u[i][j] * (p[i+1][j] - p[i-1][j])
+ v[i][j] * (p[i][j+1] - p[i][j-1])
+ p[i][j] * (u[i+1][j] - u[i-1][j])
+ p[i][j] * (v[i][j+1] - v[i][j-1]) );
u2[i][j] = u[i][j] - alpha * ( u[i][j] * (u[i+1][j] - u[i-1][j])
+ v[i][j] * (u[i][j+1] - u[i][j-1]) );
v2[i][j] = v[i][j] - alpha * ( u[i][j] * (v[i+1][j] - v[i-1][j])
+ v[i][j] * (v[i][j+1] - v[i][j-1]) );
}
for (i = 1; i < n-1; i++)
for (j = topography[i] + 1; j < n-1; j++) {
p[i][j] = p2[i][j] - alpha * ( u2[i][j] * (p2[i+1][j] - p2[i-1][j])
+ v2[i][j] * (p2[i][j+1] - p2[i][j-1])
+ p2[i][j] * (u2[i+1][j] - u2[i-1][j])
+ p2[i][j] * (v2[i][j+1] - v2[i][j-1]) );
u[i][j] = u2[i][j] - alpha * ( u2[i][j] * (u2[i+1][j] - u2[i-1][j])
+ v2[i][j] * (u2[i][j+1] - u2[i][j-1]) );
v[i][j] = v2[i][j] - alpha * ( u2[i][j] * (v2[i+1][j] - v2[i-1][j])
+ v2[i][j] * (v2[i][j+1] - v2[i][j-1]) );
}
//if (!(x%10)) { ;//rungraphic program
for (i = 0; i < n; i++)
for (j = topography[i]; j < n; j++) {
if (p[i][j] == 0) c = 0;
else { /* this will be LIC
k = sqrt(pow(u[n*i + j], 2) + pow(v[n*i + j], 2));
a = i; b = j; s = 0;
for (c = 0; c < k; c++) {
s = noise[a][b];
Yet to do here
}
c = s/k; */
c = p[i][j];
c = c % 14 +1;
}
putpixel(i, n-j, c);
}
//}// end of graphics
// settextjustify(n+20, n/2);
// cout << (x+.5);
}// end main loop
return 0;
}
inline void initializeg () {
int gdriver = VGA, gmode = VGAHI, errorcode;
/* initialize graphics and local variables */
initgraph(&gdriver, &gmode, "");
/* read result of initialization */
errorcode = graphresult();
if (errorcode != grOk) { /* an error occurred */
printf("Grap2hics error: %s\n", grapherrormsg(errorcode));
printf("Press any key to halt:");
getch();
exit(1); /* terminate with an error code */
}
}
Another Version
//Program: idealfluid1darray2.cpp
//Author: Devin Green
//Date: 11/15/1997
//Description: flow program 2d for now
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <math.h>
#include <graphics.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
void initializeg () ;
int main () {
volatile int c, i, j, x; //used for referencing
//volatile int k, s, a, b; //LIC
const int n = 100, n2 = 10000; //array size
float stepXY, stepT, alpha;
int topography[n], noise[n][n];//noise is for LIC graphics
//density and velocity arrays
float *temp, *p, *p2, *u, *u2, *v, *v2;//, *w, *w2;
p = new float[n2]; p2 = new float[n2];
u = new float[n2]; u2 = new float[n2];
v = new float[n2]; v2 = new float[n2];
// w = new float[n2]; w2 = new float[n2];
cout << "Enter the XY step in meters: "; cin >> stepXY;
cout << "Enter the time step in seconds: "; cin >> stepT;
alpha = stepT / 2 / stepXY;
const float g = 9.8 * stepT; // use SI units
initializeg();
for (x = 2; x <= 15; x++) {
setrgbpalette(x, x*4, 0, 0);
setpalette(x, x);
}
setrgbpalette(0, 0, 12, 40);
setpalette(0,0);
setrgbpalette(1, 0, 50, 0);
setpalette(1,1);
randomize();
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
noise[i][j] = random(15) + 1;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++) {
p[n*i + j] = 0; p2[n*i + j] = 0;
u[n*i + j] = 0; u2[n*i + j] = 0;
v[n*i + j] = 0; v2[n*i + j] = 0;
}
for (x = 0; x < n; x++)
topography[x] = n/2 - abs(x - n/2); //10 + floor(10*sin(x*3.14/10));
for (i = 40; i < 60; i++)
for (j = 70; j < n; j++) {
p[n*i + j] = 5;
}
for (x = 0; x < n; x++)
line(x, n, x, n-topography[x]);
for (x = 0; !kbhit(); x++) {
for (i = 1; i < n-1; i++)
for (j = topography[i] + 1; j < n-1 ; j++) {
p2[n*i + j] = p[n*i + j] - alpha*( u[n*i + j] * (p[n*(i+1) + j] - p[n*(i-1) + j])
+ v[n*i + j] * (p[n*i + j+1] - p[n*i + j-1])
+ p[n*i + j] * (u[n*(i+1) + j] - u[n*(i-1) + j])
+ p[n*i + j] * (v[n*i + j+1] - v[n*i + j-1]));
u2[n*i + j] = u[n*i + j] - alpha * ( u[n*i + j] * (u[n*(i+1) + j] - u[n*(i-1) + j]) );
//+ v[n*i + j] * (u[n*i + j+1] - u[n*i + j-1]));
v2[n*i + j] = -g + v[n*i + j] - alpha * ( u[n*i + j] * (v[n*(i+1) + j] - v[n*(i-1) + j])
+ v[n*i + j] * (v[n*i + j+1] - v[n*i + j-1]));
}
for (i = 1; i < n-1; i++) {
j = topography[i];
p2[n*i + j] = p[n*i + j] - alpha*( u[n*i + j] * (p[n*(i+1) + j] - p[n*(i-1) + j])
//+ v[n*i + j] * (p[n*i + j+1] - p[n*i + j-1])
+ p[n*i + j] * (u[n*(i+1) + j] - u[n*(i-1) + j]) );
//+ p[n*i + j] * (v[n*i + j+1] - v[n*i + j-1]));
u2[n*i + j] = u[n*i + j] - alpha * ( u[n*i + j] * (u[n*(i+1) + j] - u[n*(i-1) + j])
+ v[n*i + j] * (u[n*i + j+1] - u[n*i + j-1]));
//v2[n*i + j] = -g + v[n*i + j] - alpha * ( u[n*i + j] * (v[n*(i+1) + j] - v[n*(i-1) + j])
// + v[n*i + j] * (v[n*i + j+1] - v[n*i + j-1]));
}
temp = p; p = p2; p2 = temp;
temp = u; u = u2; u2 = temp;
temp = v; v = v2; v2 = temp;
if (!(x%10)) { ;//rungraphic program
for (i = 0; i < n; i++)
for (j = topography[i]; j < n; j++) {
if (p[n*i + j] == 0) c = 0;
else { /* this will be LIC
k = sqrt(pow(u[n*i + j], 2) + pow(v[n*i + j], 2));
a = i; b = j; s = 0;
for (c = 0; c < k; c++) {
s = noise[a][b];
Yet to do here
}
c = s/k; */
c = p[n*i + j];
c = c % 14 +1;
}
putpixel(i, n-j, c);
}
}// end of graphics
}//end of mainloop
delete[] temp, p, p2, u, u2, v, v2;//, *w, *w2;
return 0;
}//end of main
RESULTS
We hope that the output of the code will be a numerically accurate representation of a volcanic eruption. We then will display these numbers in a VRML program. By picking a time interval the user of the program will be able to look at the devastation and the progression at the time interval they pick. VRML allows you to look at the volcano with a three-dimensional perspective.
CONCLUSION
Drawn when we have more formal results.
ACKNOWLEDGEMENTS
We would like to thank Rob Huasman for all the time and help he has given to us above and beyond the call of duty. A thanks also to Mr. Taylor our supporting teacher. To Ken Wohletz and Greg Valentine for their outstanding help with the scientific part of our project.
REFERENCES
Harbaugh, John H. and Graeme Bonham-Carter. Computer Simulation in Geology. Wiley-Interscience. New York.
Valentine, Greg A. and Kenneth H. Wohletz. Effects of Topography on Facies and Compositional Zonation in Caldera-Related Ignimbrites.
Valentine, Greg A. and Kenneth H. Wohletz. Numerical Models of Plinian Eruptions Columns and Pyroclastic Flows. Los Alamos National Laboratory.
Wood, Charles A. and Jurgen Kienle. Volcanoes of North America. Cambridge University Press.