| 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 | #eta = func(X,Y,kappad) |
|---|
| 42 | #im1 = imshow(eta, cmap=cm.jet, interpolation='nearest') |
|---|
| 43 | #plot(X/lam0,eta) |
|---|
| 44 | #axis([0, 5, -2, 2]) |
|---|
| 45 | |
|---|
| 46 | figure(1) |
|---|
| 47 | plot(X/lam0,eta1,'b.-',X/lam0,eta2,'g.-') |
|---|
| 48 | title('Surface elevation of SMF: kappa = 0.83 and 1') |
|---|
| 49 | xlabel('x/lambda_0') |
|---|
| 50 | ylabel('eta(x,y=0)') |
|---|
| 51 | axis([0, 5, -2, 2]) |
|---|
| 52 | |
|---|
| 53 | def 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 | |
|---|
| 58 | def set_integral_zero_2(x,x0): |
|---|
| 59 | numerator = erf( (x-x0)/sqrt(lam0) ) |
|---|
| 60 | xg = d/tan(radians(theta)) + T/radians(sin(theta)) |
|---|
| 61 | #print 'xg', xg |
|---|
| 62 | denominator = erf( (x-xg-2*x0)/sqrt(lam0) ) |
|---|
| 63 | return numerator/denominator |
|---|
| 64 | |
|---|
| 65 | n = 100 |
|---|
| 66 | kappad_vec1 = [] |
|---|
| 67 | dx_vec1 = [] |
|---|
| 68 | step = 0.001 |
|---|
| 69 | dx = -step #+ 0.016 |
|---|
| 70 | for i in range(n): |
|---|
| 71 | dx += step |
|---|
| 72 | kappad1 = float(max(set_integral_zero(x,dx))) |
|---|
| 73 | dx_vec1.append(dx) |
|---|
| 74 | kappad_vec1.append(kappad1) |
|---|
| 75 | #print 'dx: %0.5f kappad: %f' % (dx, kappad1) |
|---|
| 76 | |
|---|
| 77 | figure(2) |
|---|
| 78 | plot(dx_vec1,kappad_vec,'.-') |
|---|
| 79 | title('Volume conservation') |
|---|
| 80 | xlabel('delta x') |
|---|
| 81 | ylabel('kappa_d') |
|---|
| 82 | savefig('Volume_conservation') |
|---|
| 83 | #axis([0,0.1,0,80]) |
|---|
| 84 | |
|---|
| 85 | n = 60 |
|---|
| 86 | kappad_vec2 = [] |
|---|
| 87 | x0_vec2 = [] |
|---|
| 88 | step = 0.1 |
|---|
| 89 | x0 = -step |
|---|
| 90 | for i in range(n): |
|---|
| 91 | x0 += step |
|---|
| 92 | kappad2 = float(max(set_integral_zero_2(x,x0))) |
|---|
| 93 | x0_vec2.append(x0) |
|---|
| 94 | kappad_vec2.append(kappad2) |
|---|
| 95 | print 'dx: %0.5f kappad: %f' % (x0, kappad2) |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | figure(3) |
|---|
| 99 | plot(x0_vec2,kappad_vec2,'.-') |
|---|
| 100 | title('Volume conservation - using xg') |
|---|
| 101 | xlabel('x0') |
|---|
| 102 | ylabel('kappa_d') |
|---|
| 103 | savefig('Volume_conservation_xg') |
|---|
| 104 | |
|---|
| 105 | |
|---|
