Newer
Older
Mads M. Pedersen
committed
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
import numpy as np
from py_wake.site.distance import StraightDistance
from numpy import newaxis as na
from py_wake.utils.grid_interpolator import GridInterpolator
class ISONoiseModel:
def __init__(self, src_x, src_y, src_h, freqs, sound_power_level, elevation_function=None):
"""ISONoise model based on
DSF/ISO/DIS 9613-2
Acoustics – Attenuation of sound during propagation – Part 2:
Engineering method for the prediction of sound pressure levels outdoors
DS/ISO 9613-1:1993
Akustik. Måling og beskrivelse af ekstern støj. Lydudbredelsesdæmpning udendørs. Del 1:
Metode til beregning af luftabsorption
The code is a vectorized version of the implementaion made by Camilla Marie Nyborg <cmny@dtu.dk>
The model models the sound pressure level from a number of sound sources at a number of receivers taking into
account
- spherical geometrical spreading
- ground reflection/absorption
- atmospheric absorption (DS/ISO 9613-1:1993)
Parameters
----------
src_x : array_like
x coordinate of sound sources
src_y : array_like
y coordinate of sound sources
src_h : array_like or float
height of sound sources (typically wind turbine hub height)
freqs : array_like
Frequencies of sound_power_level
sound_power_level : array_like
Emitted sound power level, dim=(n_sources, n_freq)
"""
if elevation_function:
src_h += elevation_function(src_x, src_y)
self.elevation_function = elevation_function
self.src_x, self.src_y = src_x, src_y
self.src_h = np.zeros_like(src_x) + src_h
self.n_src = len(src_x)
self.distance = StraightDistance()
self.freqs = freqs
self.sound_power_level = sound_power_level
def ground_eff(self, ground_distance_ij, distance_ij, ground_type):
# Ground effects ISO
freqs_init = np.array([63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 8000.0])
src_h_ij = self.src_h[:, na]
rec_h_ij = self.distance.dst_h_j[na]
hei_check_ij = 30.0 * (src_h_ij + rec_h_ij)
q_ij = np.where(ground_distance_ij <= hei_check_ij, 0, 1 - hei_check_ij / ground_distance_ij)
Gs = Gr = Gm = ground_type
# for f=63Hz Am =-3*q, for other frequencies Am = -3 * q * (1-Gm)
fak = np.where(freqs_init == 63.0, 1, (1 - Gm))
Am_fij = -3.0 * q_ij[na] * fak[:, na, na]
# 125 Hz
a_source_ij = 1.5 + 3.0 * np.exp(-0.12 * (src_h_ij - 5.0)**2) * (1.0 - np.exp(-ground_distance_ij / 50.0)) + \
5.7 * np.exp(-0.09 * src_h_ij**2) * \
(1.0 - np.exp(-2.8 * (10.0**(-6.0)) * (ground_distance_ij**2)))
a_rec_ij = 1.5 + 3.0 * np.exp(-0.12 * (rec_h_ij - 5.0)**2) * (1.0 - np.exp(-ground_distance_ij / 50.0)) + \
5.7 * np.exp(-0.09 * rec_h_ij**2) * (1.0 - np.exp(-2.8 * (10.0**(-6.0)) * (ground_distance_ij**2)))
# 250 Hz
b_source_ij = 1.5 + 8.6 * np.exp(-0.09 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
b_rec_ij = 1.5 + 8.6 * np.exp(-0.09 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
# 500 Hz
c_source_ij = 1.5 + 14.0 * np.exp(-0.46 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
c_rec_ij = 1.5 + 14.0 * np.exp(-0.46 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
# 1000 Hz
d_source_ij = 1.5 + 5.0 * np.exp(-0.9 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
d_rec_ij = 1.5 + 5.0 * np.exp(-0.9 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
zeros = np.zeros_like(ground_distance_ij)
As_fij = np.array([zeros - 1.5, # 63Hz
-1.5 + Gs * a_source_ij, # 125Hz
-1.5 + Gs * b_source_ij, # 250Hz
-1.5 + Gs * c_source_ij, # 500Hz
-1.5 + Gs * d_source_ij, # 1000Hz
zeros - 1.5 * (1.0 - Gs), # 2000Hz
zeros - 1.5 * (1.0 - Gs) # 8000Hz
])
Ar_fij = np.array([zeros - 1.5, # 63 Hz
-1.5 + Gr * a_rec_ij, # 125 Hz
-1.5 + Gr * b_rec_ij, # 250 Hz
-1.5 + Gr * c_rec_ij, # 500 Hz
-1.5 + Gr * d_rec_ij, # 1000 Hz
zeros - 1.5 * (1.0 - Gr), # 2000 Hz
zeros - 1.5 * (1.0 - Gr) # 8000 Hz
])
# Interpolation to other frequencies are actually not following the standard completely.
ip = GridInterpolator([freqs_init], np.moveaxis([As_fij, Ar_fij, Am_fij], 0, -1))
# interpolate to freq and sum up As+Ar+Am
Agr_ijf = np.moveaxis(ip(np.array([self.freqs]).T).sum(-1), 0, -1)
# Area of sphere: 4pi R^2
# 10*log_10(4pi) ~ 11
# 10 * log_10(R^2) = 20 * log_10(R) =
# reference area = 1.0 m^2
Adiv_ij = 20.0 * np.log10(distance_ij / 1.0) + 11.0 # The geometrical spreading (divergence)
ISO_ground_ijf = - Adiv_ij[:, :, na] - Agr_ijf
return ISO_ground_ijf
def atmab(self, distance_ij, T0, RH0):
# Atmospheric absorption
T_0 = 293.15
T_01 = 273.16
T = T0 + 273.15
p_s0 = 1.0
psat = p_s0 * 10.0**(-6.8346 * (T_01 / T)**1.261 + 4.6151)
h = p_s0 * (RH0) * (psat / p_s0)
F_rO = 1.0 / p_s0 * (24.0 + 4.04 * (10**4.0) * h * (0.02 + h) / (0.391 + h))
F_rN = 1.0 / p_s0 * (T_0 / T)**(0.5) * (9.0 + 280 * h *
np.exp(-4.17 * ((T_0 / T)**(1.0 / 3.0) - 1.0)))
alpha_ps = self.freqs**2.0 / p_s0 * (1.84 * (10**(-11)) * (T / T_0)**(0.5) + (T / T_0)**(-5.0 / 2.0) *
(0.01275 * np.exp(-2239.1 / T) / (F_rO + self.freqs**2.0 / F_rO) +
0.1068 * np.exp(-3352 / T) / (F_rN + self.freqs**2.0 / F_rN)))
ISO_alpha = alpha_ps[na, na] * 20.0 / np.log(10.0) * distance_ij[:, :, na]
return ISO_alpha
def transmission_loss(self, rec_x, rec_y, rec_h, ground_type, Temp, RHum):
# transmission loss = ground effects + atmospheric absorption
rec_h = np.zeros_like(rec_x) + rec_h
if self.elevation_function:
rec_h += self.elevation_function(rec_x, rec_y)
self.distance.setup(self.src_x, self.src_y, self.src_h, [rec_x, rec_y, rec_h])
ground_distance_ij = np.sqrt(self.distance.dx_ij**2 + self.distance.dy_ij**2)
distance_ij = np.sqrt(self.distance.dx_ij**2 + self.distance.dy_ij**2 + self.distance.dh_ij**2)
atm_abs_ijf = self.atmab(distance_ij, T0=Temp, RH0=RHum) # The atmospheric absorption term
ground_eff_ijf = self.ground_eff(ground_distance_ij, distance_ij, ground_type)
return ground_eff_ijf - atm_abs_ijf # Delta_SPL
def __call__(self, rec_x, rec_y, rec_h, Temp, RHum, ground_type):
"""Calculate the sound pressure level at a list of reveicers
Parameters
----------
rec_x : array_like
x coordinate, [m], of receivers
rec_y : array_like
y coordinate, [m], of receivers
rec_h : array_like or float
height, [m], of receivers (typically 2m)
Temp : float
Temperature (deg celcius)
RHum : float
Relative humidity (%)
ground_type : float
Factor describing the amount of ground absorption, 0=hard reflective ground,
e.g. paving, water, ice and concrete. 1= soft porous absorbing ground,
e.g. grass, trees, other vegetation and farming land
Returns
-------
total_spl : array_like
Total sound pressure level at each receiver
spl : array_like
Sound pressure level of freqs, dim=(n_receiver, n_freq)
"""
# Computing total transmission loss
Delta_SPL_ijf = self.transmission_loss(rec_x, rec_y, rec_h, ground_type, Temp, RHum)
sound_power_ijxxf = self.sound_power_level[:, na]
I, J, F = Delta_SPL_ijf.shape
shape = [I, J] + [1] * (len(sound_power_ijxxf.shape) - len(Delta_SPL_ijf.shape)) + [F]
Delta_SPL_ijxxf = Delta_SPL_ijf.reshape(shape)
# Add negative transmission loss to emitted sound and sum over sources to get sound pressure level
spl_jxxf = 10.0 * np.log10(np.sum(10.0**((self.sound_power_level[:, na] + Delta_SPL_ijxxf) / 10.0), axis=0))
# Sum over frequencies to get total sound pressure level
total_spl_jxx = 10.0 * np.log10(np.sum(10.0**(spl_jxxf / 10.0), axis=-1))
return total_spl_jxx, spl_jxxf
def main():
if __name__ == '__main__':
import matplotlib.pyplot as plt
from py_wake.deficit_models.gaussian import ZongGaussian
from py_wake.flow_map import XYGrid
from py_wake.turbulence_models.crespo import CrespoHernandez
from py_wake.site._site import UniformSite
from py_wake.examples.data.swt_dd_142_4100_noise.swt_dd_142_4100 import SWT_DD_142_4100
from py_wake.utils.layouts import rectangle
from py_wake.utils.plotting import setup_plot
wt = SWT_DD_142_4100()
wfm = ZongGaussian(UniformSite(), wt, turbulenceModel=CrespoHernandez())
x, y = rectangle(5, 5, 5 * wt.diameter())
sim_res = wfm(x, y, wd=270, ws=8, mode=0)
nm = sim_res.noise_model()
ax1, ax2, ax3 = plt.subplots(1, 3, figsize=(16, 4))[1]
sim_res.flow_map().plot_wake_map(ax=ax1)
ax1.plot([x[0]], [1000], '.', label='Receiver 1')
ax1.plot([x[-1]], [1000], '.', label='Receiver 2')
ax1.legend()
total_sp_jlk, spl_jlkf = nm(rec_x=[x[0], x[-1]], rec_y=[1000, 1000], rec_h=2, Temp=20, RHum=80, ground_type=0.0)
ax2.plot(nm.freqs, spl_jlkf[0, 0, 0], label='Receiver 1')
ax2.plot(nm.freqs, spl_jlkf[1, 0, 0], label='Receiver 2')
setup_plot(xlabel='Frequency [Hz]', ylabel='Sound pressure level [dB]', ax=ax2)
plt.tight_layout()
nmap = sim_res.noise_map(grid=XYGrid(x=np.linspace(-1000, 5000, 100), y=np.linspace(-1000, 1000, 50), h=2))
nmap['Total sound pressure level'].squeeze().plot(ax=ax3)
wt.plot(x, y, ax=ax3)
plt.show()
main()