# source:production/sydney_2006/smf_volume_conservation.py@2673

Last change on this file since 2673 was 2673, checked in by sexton, 18 years ago

Script to determine slump parameters to ensure volume conservation

File size: 2.3 KB
Line
1"""
2Determining volume conservation for SMF slump formulation.
3Jane Sexton, GA 2006
4"""
5
6from math import sqrt, atan, degrees, exp, cosh, pi, radians
7from Numeric import ones
8from pylab import *
9import scipy
10from scipy.special import erf
11
12ion()
13hold(False)
14
15
16x = arange(-20,25,0.1)
17y = arange(-10,10,0.1)
18X,Y = meshgrid(x,y)
19
20# Grilli and Watts 2005 example
21dx = 2.0
22x0 = 10.0
23y0 = 0.0
24w = 2.0
25kappa = 3
26g = 9.8
27d = 0.259
28T = 0.052
30lam0 = 5.0 #sqrt(g*d)
31
33
34    numerator = exp(-((x-x0)/lam0)**2) - kappad*exp(-((x-dx-x0)/lam0)**2)
35    denominator = cosh(kappa*(y-y0)/(w+lam0))**2
36    return -numerator/denominator
37
38eta1 = func(X,0,0.83)
39eta2 = func(X,0,1)
40
42#im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
43#plot(X/lam0,eta)
44#axis([0, 5, -2, 2])
45
46figure(1)
47plot(X/lam0,eta1,'b.-',X/lam0,eta2,'g.-')
48title('Surface elevation of SMF: kappa = 0.83 and 1')
49xlabel('x/lambda_0')
50ylabel('eta(x,y=0)')
51axis([0, 5, -2, 2])
52
53def set_integral_zero(x,dx):
54    numerator = erf( (x-x0)/sqrt(lam0) )
55    denominator = erf( (x-dx-x0)/sqrt(lam0) )
56    return numerator/denominator
57
58def set_integral_zero_2(x,x0):
59    numerator = erf( (x-x0)/sqrt(lam0) )
61    #print 'xg', xg
62    denominator = erf( (x-xg-2*x0)/sqrt(lam0) )
63    return numerator/denominator
64
65n = 100
67dx_vec1 = []
68step = 0.001
69dx = -step #+ 0.016
70for i in range(n):
71    dx += step
73    dx_vec1.append(dx)
76
77figure(2)
79title('Volume conservation')
80xlabel('delta x')
81ylabel('kappa_d')
82savefig('Volume_conservation')
83#axis([0,0.1,0,80])
84
85n = 60
87x0_vec2 = []
88step = 0.1
89x0 = -step
90for i in range(n):
91    x0 += step
93    x0_vec2.append(x0)