/
gplume_stability.py
128 lines (100 loc) · 3.56 KB
/
gplume_stability.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import numpy as np
class receptorGrid:
def __init__(self,x,y,z):
self.x = x
self.y = y
self.z = z
self.yMesh,self.zMesh,self.xMesh = np.meshgrid(y,z,x)
class pointSource:
def __init__(self,x,y,z,rate,H):
self.x = x
self.y = y
self.z = z
self.rate = rate
self.H = H
self.sourceType = 'point'
class areaSource:
def __init__(self,x0,dx,nx,y0,dy,ny,z,rate,H):
self.x = np.linspace(x0,nx*dx,nx+1)
self.y = np.linspace(y0,ny*dy,ny+1)
self.z = z
self.sourceType = 'area'
self.yMesh,self.zMesh,self.xMesh = np.meshgrid(self.y,z,self.x)
self.H = H
self.rate = rate
self.dx = dx
self.dy = dy
class stabilityClass:
def __init__(self,letter):
self.letter = letter
if letter == 'A':
Iy = -1.104
Jy = 0.9878
Ky = -0.0076
Iz = 4.679
Jz = -1.7172
Kz = 0.2770
elif letter == 'B':
Iy = -1.634
Jy = 1.0350
Ky = -0.0096
Iz = -1.999
Jz = 0.8752
Kz = 0.0136
elif letter == 'C':
Iy = -2.054
Jy = 1.0231
Ky = -0.0076
Iz = -2.341
Jz = 0.9477
Kz = -0.0020
elif letter == 'D':
Iy = -2.555
Jy = 1.0423
Ky = -0.0087
Iz = -3.186
Jz = 1.1737
Kz = -0.0316
elif letter_ == 'E':
Iy = -2.754
Jy = 1.0106
Ky = -0.0064
Iz = -3.783
Jz = 1.3010
Kz = -0.0450
elif letter == 'F':
Iy = -3.143
Jy = 1.0148
Ky = -0.0070
Iz = -4.490
Jz = 1.4024
Kz = -0.0540
def sy(dist):
return np.exp(Iy + Jy*np.log(dist) + Ky*(np.log(dist)**2))
def sz(dist):
return np.exp(Iz + Jz*np.log(dist) + Kz*(np.log(dist)**2))
self.sz = sz
self.sy = sy
class gaussianPlume:
def __init__(self,source,grid,stability,U):
self.grid = grid
self.source = source
self.stability = stability
self.U = U
def calculateConcentration(self):
conc = np.zeros_like(self.grid.xMesh,dtype=float)
if self.source.sourceType == 'area':
for x in self.source.x:
for y in self.source.y:
a = self.source.rate*self.source.dx*self.source.dy / (2 * np.pi * self.U * self.stability.sy(self.grid.xMesh -x) * self.stability.sz(self.grid.xMesh - x))
b = np.exp(-(self.grid.yMesh - y)**2/(2*self.stability.sy(self.grid.xMesh -x)**2))
c = np.exp(-(self.grid.zMesh-self.source.H)**2/(2*self.stability.sz(self.grid.xMesh - x)**2)) + np.exp(-(self.grid.zMesh+self.source.H)**2/(2*self.stability.sz(self.grid.xMesh - x)**2))
conc += a*b*c
if self.source.sourceType == 'point':
x = self.source.x
y = self.source.y
a = self.source.rate / (2 * np.pi * self.U * self.stability.sy(self.grid.xMesh -x) * self.stability.sz(self.grid.xMesh - x))
b = np.exp(-(self.grid.yMesh - y)**2/(2*self.stability.sy(self.grid.xMesh -x)**2))
c = np.exp(-(self.grid.zMesh-self.source.H)**2/(2*self.stability.sz(self.grid.xMesh - x)**2)) + np.exp(-(self.grid.zMesh+self.source.H)**2/(2*self.stability.sz(self.grid.xMesh - x)**2))
conc += a*b*c
return conc