wiat.TM.Radiation

1 # coding: latin-1 2 # Copyright 2008 Antoine Lefebvre 3 # 4 # This file is part of "The Wind Instrument Acoustic Toolkit". 5 # 6 # "The Wind Instrument Acoustic Toolkit" is free software: 7 # you can redistribute it and/or modify it under the terms of 8 # the GNU General Public License as published by the 9 # Free Software Foundation, either version 3 of the License, or 10 # (at your option) any later version. 11 # 12 # "The Wind Instrument Acoustic Toolkit" is distributed in the hope 13 # that it will be useful, but WITHOUT ANY WARRANTY; without even 14 # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 15 # See the GNU General Public License for more details. 16 # 17 # You should have received a copy of the GNU General Public License 18 # along with "The Wind Instrument Acoustic Toolkit". 19 # If not, see <http://www.gnu.org/licenses/>. 20 # 21 # All rights reserved. 22 23 """ Radiation.py - Classes defining the radiation impedance models 24 25 """ 26 27 __all__ = ['UnflangedOpenEnd','FlangedOpenEnd','IdealOpenEnd','ClosedEnd', 28 'SaxophoneBell'] 29 30 import os 31 import math 32 33 import numpy 34 from numpy import log, sqrt 35 36 import scipy.integrate 37 from scipy.integrate import quad 38 39 import scipy.special 40 from scipy.special import jn_zeros, j1 41 42 import scipy.interpolate 43 from scipy.interpolate import UnivariateSpline 44 45 import wiat 46 from wiat.util import air_properties, transform_R2Z

47 48 # in this file we can add radiation models for standard horns 49 # possibly pre-calculated from BEM or MM models. 50 # 51 # Data should then be loaded from a file and impedance interpolated 52 # for tha values we need 53 54 -class UnflangedOpenEnd():

55 """ Unflanged open end radiation impedance models 56 57 References: 58 59 1. J. Dalmont and C.J. Nederveen, 60 "Radiation impedance of tubes with different flanges: 61 numerical and experimental investigations ," 62 Journal of Sound and Vibration, vol. 244, pp. 505-534, 2001. 63 2. R. Caussé, J. Kergomard, and X. Lurton, 64 "Input impedance of brass musical instruments - 65 Comparaison between experimental and numerical models," 66 J. Acoust. Soc. Am., vol. 75, pp. 241-254, 1984. 67 3. H. Levine and J. Schwinger, 68 On the radiation of sound from an unflanged circular pipe. 69 Phys. Rev., 73(4), pp. 383-406, 1948. 70 71 TODO: implement the exact solution from Levine and Schwinger 72 73 """ 74

75 - def __init__(self,d,model='Dalmont_Nederveen',m=1.):

76 """ Radiation from an unflanged open pipe (Levine and schwinger) 77 78 @param d: diameter of the open end 79 @param m: multiplication factor (to take care of sperical wave, 80 approximation from Causse) 81 """ 82 self.a = d/2. 83 self.m = m 84 self.model = model 85 self.Z = {'Dalmont_Nederveen': self.Z_Dalmont_Nederveen, 86 'Causse': self.Z_Causse, 87 'Levine_Schwinger': self.Z_Levine_Schwinger}

88

89 - def __call__(self,f,T):

90 return self.Z[self.model](f,T)

91

92 - def TM(self,f,T):

93 return numpy.array([[self.Z[self.model](f,T)], 94 [numpy.ones(numpy.size(f))]],dtype=complex)

95

96 - def Z_Dalmont_Nederveen(self,f,T):

97 """ Unflanged pipe radiation impedance approximation (ka < 3.5) from 98 Dalmont and Nederveen 99 100 R = -|R0|exp(-2j*ka*delta), where delta is the frequency dependant 101 length correction 102 103 See reference [1] page 509 104 """ 105 air = air_properties(T) 106 ka = 2*numpy.pi*f*self.a/air['c'] 107 R0 = (1 + 0.2*ka - 0.084*ka**2)/(1 + 0.2*ka + (0.5 - 0.084)*ka**2) 108 delta = 0.6113 * ((1+0.044*ka**2)/(1+0.19*ka**2) - 109 0.02*(numpy.sin(2*ka))**2 ) 110 R = - R0 * numpy.exp( -2 * 1j * ka * delta) 111 Zl = transform_R2Z(R) 112 return self.m*Zl

113

114 - def Z_Causse(self,f,T):

115 """ Unflanged pipe radiation impedance approximation (ka < 1.5) 116 from Causse & al. 117 See [2] 118 """ 119 air = air_properties(T) 120 ka = 2*numpy.pi*f*self.a/air['c'] 121 Zl = 0.25*(ka)**2 + 0.0127*(ka)**4 + 0.082*(ka)**4 * numpy.log(ka) -\ 122 0.023*(ka)**6 + 1j * (0.6113*(ka) - 0.036*(ka)**3 + \ 123 0.034*(ka)**3 * numpy.log(ka) - 0.0187*(ka)**5) 124 return Zl

125

126 - def Z_Levine_Schwinger(self,f,T):

127 air = air_properties(T) 128 ka = 2*numpy.pi*f*self.a/air['c'] 129 130 #z = ka 131 132 #int 2: Eq. (VI.4) of Levine and Schwinger 133 s = numpy.empty((len(ka))) 134 for i,z in enumerate(ka): 135 s[i] = scipy.integrate.quad( 136 lambda x,z: log(1/(2*scipy.special.i1(x)*\ 137 scipy.special.k1(x)))/(x*scipy.sqrt(x**2 + z**2)), 138 0,20,(z,))

139

140 #int 3: Eq. (VI.5) of Levine and Schwinger 141 # lambda x,z: numpy.arctan( scipy.special.k1(x)/(numpy.pi*scipy.special.i1(x)) ) * \ 142 # (1 - z/(numpy.sqrt(z**2 + x**2))) / x 143 #int 4: Eq. (VI.4) of Levine and Schwinger 144 # lambda x,z: nump.log(numpy.pi*numpy.abs(scipy.special.j1(x))*\ 145 # numpy.sqrt(scipy.special.j1(x)**2 + scipy.special.y1(x)**2))/\ 146 # (x*numpy.sqrt(z**2 - x**2)); 147 148 -class FlangedOpenEnd():

149 - def __init__(self,d,model='Dalmont_Nederveen'):

150 151 self.a = d/2 152 self.model = model 153 self.Z = {'Dalmont_Nederveen': self.Z_Dalmont_Nederveen}

154

155 - def __call__(self,f,T):

156 return self.Z[self.model](f,T)

157

158 - def TM(self,f,T):

159 return numpy.array([[self.Z[self.model](f,T)],[numpy.ones(numpy.size(f))]],dtype=complex)

160

161 - def Z_Dalmont_Nederveen(self,f,T):

162 """ Flanged pipe radiation impedance approximation (ka < 3.5) 163 from Dalmont and Nederveen 164 165 R = -|R0|exp(-2j*ka*delta), where delta is the frequency dependant 166 length correction 167 168 Z = (1+R)/(1-R) 169 See [1] page 509 170 """ 171 air = air_properties(T) 172 k = 2*numpy.pi*f/air['c'] 173 ka = k * self.a 174 Rinf = (1 + 0.323*ka - 0.077*ka**2)/(1 + 0.323*ka + (1 - 0.077)*ka**2) 175 #delta = 0.8216 * (1 + (0.77*ka)**2/(1+0.77*ka))**-1 # Dalmont 176 177 delta = (0.82159-0.49*ka**2)/(1 -0.46*ka**3) # Norris 178 179 R = - Rinf * numpy.exp( -2 * 1j * ka * delta) 180 Zl = transform_R2Z(R) 181 return Zl

182 183 @staticmethod

184 - def Norris_Sheng_dbeta(s, x, i):

185 return 2*s*j1(s)**2 / (sqrt(x**2-s**2) * (s**2-jn_zeros(1,i)**2))

186

187 # def Z_Norris_Sheng(self,f,T,number_of_modes = 15): 188 # """ 189 # Norris, A. N., & Sheng, I. C. (1989). 190 # Acoustic radiation from a circular pipe with an infinite flange. 191 # Journal of Sound and Vibration, 135(1), 85-93. 192 # """ 193 # 194 # air = air_properties(T) 195 # k = 2*numpy.pi*f/air['c'] 196 # ka = k * self.a 197 # 198 # 199 # alpha_00 = ka**-1 - j1(2*ka)/ka**2 200 # 201 # M_00 = alpha_00 + 1.j/(ka) # eq. 9 202 # 203 # # eq. 14, R as a function of ka 204 # sum = 0. 205 # for i in range(1,number_of_modes): 206 # # eq. 13c 207 # beta_i_1 = 1j*quad(Norris_Sheng_dbeta, 0, x, args=(ka,i)) + \ 208 # quad(Norris_Sheng_dbeta,x,scipy.inf,args=(ka,i)) 209 # alpha_i0 = beta_i_1 # eq. 13b 210 # if (ka > jn_zeros(1,i)): 211 # xhi_i = numpy.sqrt(ka**2-jn_zeros(1,i)**2)/a 212 # else: 213 # xhi_i = 1j*numpy.sqrt(jn_zeros(1,i)**2-ka**2)/a 214 # 215 # B_i = A_i*xhi_i/(2*k) # eq. 12 216 # sum = sum + alpha_i0*B_i 217 # R = -1 + 2*(alpha_00 - sum)/M_00 218 219 220 221 -class FlangedPiston():

222 pass

223

224 -class UnflangedPiston():

225 pass

226

227 -class IdealOpenEnd():

228 - def __call__(self,f,T):

229 return numpy.zeros(numpy.size(f))

230

231 - def TM(self,f,T):

232 return numpy.array([[numpy.zeros(numpy.size(f))], 233 [numpy.ones(numpy.size(f))]],dtype=complex)

234

235 -class ClosedEnd():

236 - def __call__(self,f,T):

237 return numpy.inf*numpy.ones(numpy.size(f))

238

239 - def TM(self,f,T):

240 return numpy.array([[numpy.ones(numpy.size(f))], 241 [numpy.zeros(numpy.size(f))]],dtype=complex)

242

243 244 -class SaxophoneBell():

245 """ Radiation model of a saxophone bell. 246 247 It uses the multimodal calculation method to calculate the planar 248 wave input impedance of the bell. 249 """

250 - def __init__(self, d, theta):

251 """ 252 @param d: diameter of the bell at the junction with the body (small end) 253 the diameter of the large end is implicit. 254 @param theta: opening angle of the bell at the junction with the body. 255 This is usually the same at the body angle in order to have a continuity 256 in the tangent. 257 """ 258 self.a = d/2. 259 260 261 # the calculations are calculated as a function of ka and interpolated 262 # if there is a diameter change afterward 263 self.kas = numpy.arange(0.,5.,0.1) 264 265 # We use a representative dimension for the losses to be valid 266 self.d0 = 0.1 267 self.de = 0.15 268 269 # there is no Unflanged model yet... 270 self.Zr = wiat.MM.Radiation.RadiationFlanged(self.de) 271 272 self.bell = wiat.MM.Waveguides.QuarterCircleHorn( 273 self.d0,self.de,math.radians(theta),math.radians(80.)) 274 275 self.number_of_modes = 5 276 277 self.backup_filename = os.path.join(os.path.dirname(__file__), 278 'SaxophoneBell.mat')

279

280 - def __str__(self):

281 return "SaxophoneBell :: %s, length = %.4f"%\ 282 (str(self.bell), self.bell.L)

283 284

285 - def __call__(self,f,T):

286 air = wiat.util.air_properties(T) 287 288 # these are the value of ka for which we want the impedance to be 289 # returned 290 kas = 2*numpy.pi*f*self.a/air['c'] 291 292 293 # this will create the impedance vector for the pre-defined value of 294 # ka (self.kas) 295 Z = wiat.util.load_or_calculate(self.backup_filename, 296 wiat.MM.Solver.SolveImpedance, 297 geometry=self.bell, Zr=self.Zr, 298 # ks=2*self.kas/self.d0 299 freqs=air['c']*self.kas/(self.d0*numpy.pi), 300 air=air, 301 N=self.number_of_modes) 302 303 spl_real = UnivariateSpline(self.kas, numpy.real(Z),s=0.0000001) 304 spl_imag = UnivariateSpline(self.kas, numpy.imag(Z),s=0.0000001) 305 306 return spl_real(kas) - 1j*spl_imag(kas)

307 308

309 - def TM(self,f,T):

310 return numpy.array([[self(f,T)],[numpy.ones(numpy.size(f))]], 311 dtype=complex)

312

Source Code for Module wiat.TM.Radiation