| 1 | """ |
|---|
| 2 | Determining volume conservation for SMF slump formulation. |
|---|
| 3 | Jane Sexton, GA 2006 |
|---|
| 4 | """ |
|---|
| 5 | |
|---|
| 6 | from math import sqrt, atan, degrees, exp, cosh, pi, radians |
|---|
| 7 | from Numeric import ones |
|---|
| 8 | from pylab import * |
|---|
| 9 | #import scipy |
|---|
| 10 | #from scipy.special import erf |
|---|
| 11 | |
|---|
| 12 | ion() |
|---|
| 13 | hold(False) |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | x = arange(-20,25,0.1) |
|---|
| 17 | y = arange(-10,10,0.1) |
|---|
| 18 | X,Y = meshgrid(x,y) |
|---|
| 19 | |
|---|
| 20 | # Grilli and Watts 2005 example |
|---|
| 21 | dx = 2.0 |
|---|
| 22 | x0 = 10.0 |
|---|
| 23 | y0 = 0.0 |
|---|
| 24 | w = 2.0 |
|---|
| 25 | kappa = 3 |
|---|
| 26 | g = 9.8 |
|---|
| 27 | d = 0.259 |
|---|
| 28 | T = 0.052 |
|---|
| 29 | theta = 15.0 #acos(radians(-d/T)) |
|---|
| 30 | lam0 = 5.0 #sqrt(g*d) |
|---|
| 31 | |
|---|
| 32 | def func(x,y,kappad): |
|---|
| 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 | |
|---|
| 38 | eta1 = func(X,0,0.83) |
|---|
| 39 | eta2 = func(X,0,1) |
|---|
| 40 | |
|---|
| 41 | #import sys; sys.exit() |
|---|
| 42 | #eta = func(X,Y,kappad) |
|---|
| 43 | #im1 = imshow(eta, cmap=cm.jet, interpolation='nearest') |
|---|
| 44 | #plot(X/lam0,eta) |
|---|
| 45 | #axis([0, 5, -2, 2]) |
|---|
| 46 | |
|---|
| 47 | figure(1) |
|---|
| 48 | plot(X/lam0,eta1,'b.-',X/lam0,eta2,'g.-') |
|---|
| 49 | title('Surface elevation of SMF: kappa = 0.83 and 1') |
|---|
| 50 | xlabel('x/lambda_0') |
|---|
| 51 | ylabel('eta(x,y=0)') |
|---|
| 52 | axis([0, 5, -2, 2]) |
|---|
| 53 | |
|---|
| 54 | def set_integral_zero(x,dx): |
|---|
| 55 | numerator = erf( (x-x0)/sqrt(lam0) ) |
|---|
| 56 | denominator = erf( (x-dx-x0)/sqrt(lam0) ) |
|---|
| 57 | return numerator/denominator |
|---|
| 58 | |
|---|
| 59 | def set_integral_zero_2(x,x0): |
|---|
| 60 | numerator = erf( (x-x0)/sqrt(lam0) ) |
|---|
| 61 | xg = d/tan(radians(theta)) + T/radians(sin(theta)) |
|---|
| 62 | #print 'xg', xg |
|---|
| 63 | denominator = erf( (x-xg-2*x0)/sqrt(lam0) ) |
|---|
| 64 | return numerator/denominator |
|---|
| 65 | |
|---|
| 66 | n = 100 |
|---|
| 67 | kappad_vec1 = [] |
|---|
| 68 | dx_vec1 = [] |
|---|
| 69 | step = 0.001 |
|---|
| 70 | dx = -step #+ 0.016 |
|---|
| 71 | for i in range(n): |
|---|
| 72 | dx += step |
|---|
| 73 | kappad1 = float(max(set_integral_zero(x,dx))) |
|---|
| 74 | dx_vec1.append(dx) |
|---|
| 75 | kappad_vec1.append(kappad1) |
|---|
| 76 | #print 'dx: %0.5f kappad: %f' % (dx, kappad1) |
|---|
| 77 | |
|---|
| 78 | figure(2) |
|---|
| 79 | plot(dx_vec1,kappad_vec,'.-') |
|---|
| 80 | title('Volume conservation') |
|---|
| 81 | xlabel('delta x') |
|---|
| 82 | ylabel('kappa_d') |
|---|
| 83 | savefig('Volume_conservation') |
|---|
| 84 | #axis([0,0.1,0,80]) |
|---|
| 85 | |
|---|
| 86 | n = 60 |
|---|
| 87 | kappad_vec2 = [] |
|---|
| 88 | x0_vec2 = [] |
|---|
| 89 | step = 0.1 |
|---|
| 90 | x0 = -step |
|---|
| 91 | for i in range(n): |
|---|
| 92 | x0 += step |
|---|
| 93 | kappad2 = float(max(set_integral_zero_2(x,x0))) |
|---|
| 94 | x0_vec2.append(x0) |
|---|
| 95 | kappad_vec2.append(kappad2) |
|---|
| 96 | print 'dx: %0.5f kappad: %f' % (x0, kappad2) |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | figure(3) |
|---|
| 100 | plot(x0_vec2,kappad_vec2,'.-') |
|---|
| 101 | title('Volume conservation - using xg') |
|---|
| 102 | xlabel('x0') |
|---|
| 103 | ylabel('kappa_d') |
|---|
| 104 | savefig('Volume_conservation_xg') |
|---|
| 105 | |
|---|
| 106 | |
|---|
