wiat.TM.Waveguides

1 # coding: latin-1 2 # 3 # Copyright 2008 Antoine Lefebvre 4 # 5 # This file is part of "The Wind Instrument Acoustic Toolkit". 6 # 7 # "The Wind Instrument Acoustic Toolkit" is free software: 8 # you can redistribute it and/or modify it under the terms of 9 # the GNU General Public License as published by the 10 # Free Software Foundation, either version 3 of the License, or 11 # (at your option) any later version. 12 # 13 # "The Wind Instrument Acoustic Toolkit" is distributed in the hope 14 # that it will be useful, but WITHOUT ANY WARRANTY; without even 15 # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 16 # See the GNU General Public License for more details. 17 # 18 # You should have received a copy of the GNU General Public License 19 # along with "The Wind Instrument Acoustic Toolkit". 20 # If not, see <http://www.gnu.org/licenses/>. 21 # 22 # All rights reserved. 23 24 """ Waveguides.py - Contains classes defining the transmission matrices for 25 simple geometries and functions to generate them. 26 27 """ 28 29 __all__ = ['Cylinder','Cone','CurvedCone', 30 'ConeApproximatedWithCylinders','EquationApproximatedWithCones', 31 'CreateCylindersAndConesFromList', 'Jump', 'ToneHole'] 32 33 from scipy import special 34 from scipy.integrate import quad 35 36 import numpy 37 from numpy import size, sqrt, cosh, sinh, tanh, ones 38 39 import Radiation 40 from wiat.util import air_properties 41

42 -def PropagationConstant(f,a,air,beta):

43 """ Calculate the complex propagation constant in an acoustic waveguide 44 45 See Allan D. Pierce p. 531-534 46 47 @param f: frequency array 48 @param a: tube radius 49 @param air: Air instance 50 @param beta: factor to take into account wall roughness 51 @return: S{Gamma} = S{alpha} + 1j * S{omega} / v_S{phi} 52 where S{alpha} is the attenuation, S{omega} is 2S{pi}f and 53 v_S{phi} is the phase velocity 54 55 """ 56 omega = 2*numpy.pi*f 57 58 # slightly modified version of eq. 10-5.8 59 alpha = beta * (1/a) * sqrt(omega*air['mu']/(2*air['rho']*air['c']**2)) *\ 60 (1 + (air['gamma']-1)/air['nu']) 61 62 # slightly modified version of eq. 10-5.10 63 G = 1j*(omega/air['c']) + (1+1j)*alpha 64 return G

65 66

67 -class Waveguides:

68 " Base class for waveguides "

69 - def __init__(self):

70 self.a1 = 0. 71 self.L = 0.

72

73 - def get_input_surface(self):

74 " Return the input surface of the section. " 75 return numpy.pi * self.a1 ** 2

76

77 - def get_length(self):

78 " return the length of the section. " 79 return self.L

80

81 - def checkx(self,x):

82 assert -1e-6 <= x <= self.L+1e-6,\ 83 "value of x outside the waveguide by %.4f"%(x-self.L)

84

85 -class Cylinder(Waveguides):

86 " Class representing a cylinder "

87 - def __init__(self,d,L):

88 """ Cylinder transmission matrix calculation class 89 90 d: diameter 91 L: length 92 """ 93 self.a1 = d/2. 94 self.L = L

95

96 - def __str__(self):

97 return "Cylinder(d=%.4f,L=%.4f)"%(2*self.a1,self.L)

98

99 - def radius(self,x):

100 self.checkx(x) 101 return self.a1

102

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

104 " Transmission matrix of the complete cylinder " 105 TM = self.subTM(f,T,0.,self.L) 106 return TM

107

108 - def subTM(self,f,T,xi,xf):

109 """ Transmission matrix of a section of the cylinder 110 f: array of frequencies 111 T: temperature 112 xi: beginning of the section 113 xf: end of the section 114 """ 115 assert 0. <= xi <= self.L,\ 116 "cylinder subsection out of range %f <= %f <= %f"%(0.,xi,self.L) 117 assert xi <= xf <= self.L+0.00001,\ 118 "cylinder subsection out of range %f <= %f <= %f"%(xi,xf,self.L) 119 L = xf - xi 120 air = air_properties(T) 121 Gamma = PropagationConstant(f,self.a1,air,1) 122 sinhL = numpy.sinh(Gamma*L) 123 coshL = numpy.cosh(Gamma*L) 124 TM = numpy.array([[coshL,sinhL],[sinhL,coshL]]) 125 return TM

126

127 -class Cone(Waveguides):

128 """ Truncated cone transmission matrix calculation class 129 130 This is an implementation of the formulation presented in [1] 131 132 References: 133 134 1. Transfer matrix of conical waveguides with any geometric 135 parameters for increased precision in computer modeling, 136 The Journal of the Acoustical Society of America, 137 vol. 122, pp. EL179-EL184, November 2007. 138 """ 139

140 - def __init__(self,d0,de,L):

141 """ Initialize the conical waveguide data structure 142 143 @param d0: diameter at the beginning [m] 144 @type d0: float 145 @param de: diameter at the end [m] 146 @type de: float 147 @param L: length of the cone [m] 148 @type L: float 149 150 151 """ 152 153 self.a1 = d0/2. 154 self.ae = de/2. 155 self.L = L 156 157 self.x_in = self.a1 * L / ( self.ae - self.a1 ) 158 # distance from apex to beginning of cone 159 160 self.x_out = self.x_in + L

161 # distance from apex to end of cone 162

163 - def __str__(self):

164 return "Cone(d0=%.4f,de=%.4f,L=%.4f)"%(2*self.a1,2*self.ae,self.L)

165

166 - def halfangle(self):

167 " Return the half angle of the cone" 168 return numpy.arcsin(self.a1/self.x_in)

169

170 - def radius(self,x):

171 """ Radius of the cone at distance x from its beginning """ 172 self.checkx(x) 173 return (self.x_in + x) * self.a1 / self.x_in

174

175 - def _k(self,f,z):

176 """ Local propagation constant in the conical waveguide 177 at distance z from the apex 178 """ 179 self.checkx(z - self.x_in) 180 a = self.radius(z - self.x_in) 181 air = air_properties(self.T) 182 beta = 1. 183 return -1j*PropagationConstant(f,a,air,beta)

184

185 - def _kbarL(self,freqs,xi,xf):

186 """ equivalent propagation constant multiply by length of cone 187 see Kulik eq.8 188 """ 189 kbL = numpy.empty(numpy.size(freqs),dtype=complex) 190 for i,f in enumerate(freqs): 191 # we need integrate separately real and imaginary part because " 192 # quad" do not support complex numbers 193 kbarLr = quad(lambda x: numpy.real(self._k(f,x)), xi, xf, 194 epsabs=1e-10,epsrel=1e-8) 195 kbarLi = quad(lambda x: numpy.imag(self._k(f,x)), xi, xf, 196 epsabs=1e-10,epsrel=1e-8) 197 kbL[i] = kbarLr[0] +1j*kbarLi[0] 198 return kbL

199

200 - def _theta(self,f,z):

201 " z is the distance from the cone apex " 202 return numpy.arctan( self._k(f, z) * z )

203

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

205 self.T = T 206 TM = self.subTM(f,T,0.,self.x_out-self.x_in) 207 return TM

208

209 - def subTM(self,f,T,xi,xf):

210 self.T = T 211 212 x_in = self.x_in + xi 213 x_out = self.x_in + xf 214 215 theta_in = self._theta(f,x_in) 216 t_in = 1 / numpy.sin(theta_in); 217 218 theta_out = self._theta(f,x_out) 219 t_out = 1 / numpy.sin(theta_out); 220 221 r = x_out / x_in 222 kbL = self._kbarL(f,x_in,x_out) 223 224 TM = numpy.array([[-r*t_out*numpy.sin(kbL-theta_out), 225 1j*r*numpy.sin(kbL)], 226 [1j*r*t_in*t_out*numpy.sin(kbL+theta_in-theta_out), 227 r*t_in*numpy.sin(kbL+theta_in)]]) 228 return TM

229

230 -class CurvedCone(Cone):

231 """ Class representing a curved cone 232 233 The idea is to construct an equivalent straight cone. 234 The theory behind this have to be verify. 235 236 References: 237 238 1. Nederveen, C. J. (1998). 239 Influence of a toroidal bend on wind instrument tuning. 240 JASA, 104(3), 1616-1626. 241 242 """

243 - def __init__(self,d0,de,L,R):

244 """ Initialize the CurvedCone data structure 245 246 @param d0: diameter at the beginning [m] 247 @type d0: float 248 @param de: diameter at the end [m] 249 @type de: float 250 @param L: length of the cone [m] 251 @type L: float 252 @param R: the radius of curvature of the section 253 @type R: float 254 """ 255 B0 = d0/(2*R) 256 Be = de/(2*R) 257 B = (B0+Be)/2. # mean curvature 258 259 # Nederveen, 1998, equation 9 260 rhoB_rho = B**2/2 / (1 - sqrt(1 - B**2)) 261 LB = L*sqrt(rhoB_rho) 262 d0B = 2*sqrt( (d0/2)**2 * sqrt(1/rhoB_rho) ) 263 deB = 2*sqrt( (de/2)**2 * sqrt(1/rhoB_rho) ) 264 265 Cone.__init__(self,d0B,deB,LB)

266

267 -def CreateCylindersAndConesFromList(list):

268 """ Build a list of Cylinder or Cone from a list of tuples containing 269 position and diameters. 270 271 Example:: 272 >>> import wiat 273 >>> FluteBody = wiat.TM.CreateCylindersAndConesFromList( 274 ... [(0.,0.0188),(0.132,0.0188),(0.340,0.0152), 275 ... (0.452,0.0135),(0.583,0.0112)]) 276 >>> for s in FluteBody: print s 277 Cylinder(d=0.0188,L=0.1320) 278 Cone(d0=0.0188,de=0.0152,L=0.2080) 279 Cone(d0=0.0152,de=0.0135,L=0.1120) 280 Cone(d0=0.0135,de=0.0112,L=0.1310) 281 282 @param list: a list of tuples (position,diameter) where the 283 position is the distance along the centerline of the instrument 284 from the mouthpiece. 285 286 """ 287 sections = [] 288 p0,d0 = list[0] 289 for (p,d) in list[1:]: 290 if d==d0: 291 sections.append( Cylinder(d, p-p0) ) 292 else: 293 sections.append( Cone(d0,d,p-p0) ) 294 p0,d0 = p,d 295 return sections

296

297 -def ConeApproximatedWithCylinders(d0,de,L,N=200):

298 c = Cone(d0,de,L) 299 lc = L / float(N) 300 cylinders = [] 301 dp = d0 302 for i in range(N): 303 d = 2*c.radius(lc*(2*i+1)/2) 304 # diameter of the cone at the center of the section 305 306 cylinders.append( Jump(dp,d) ) 307 cylinders.append( Cylinder(d,lc) ) 308 dp = d 309 cylinders.append( Jump(d, de) ) 310 return cylinders

311

312 -def EquationApproximatedWithCones(equation,L,N=10):

313 """ Discretization in conical segments of an equation giving the radius of 314 the waveguide as a function of its argument. 315 316 @param equation: an equation of one argument. It will be evaluated from 317 0 to L. 318 @param L: the length of the waveguide. 319 """ 320 lc = L / float(N) 321 cones = [] 322 d1 = 2*equation(0) 323 for i in range(N): 324 x = (i+1)*lc 325 d2 = 2*equation(x) 326 cones.append( Cone(d1,d2,lc) ) 327 d1 = d2 328 return cones

329

330 -class Jump(Waveguides):

331 " Class representing a diameter jump "

332 - def __init__(self,d1,d2):

333 """ Transfer Matrix of a diameter jump 334 d1: diameter 1 335 d2: diameter 2 336 """ 337 self.a1 = d1/2. 338 self.a2 = d2/2. 339 self.L = 0.

340 341 # we should add two more terms 342 # 343 # inertance (La) and inductance (Ra) 344 # see [1] Morse and Ingard, Theorical Acoustics, McGraw-Hill, 1968. 345 # p. 488 346

347 - def subTM(self,f,T,xi,xf):

348 return self.TM(f,T)

349

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

351 S1 = numpy.pi * self.a1**2 352 S2 = numpy.pi * self.a2**2 353 n = numpy.size(f) 354 TM = numpy.array([[numpy.ones(n),numpy.zeros(n)], 355 [numpy.zeros(n),numpy.ones(n)*S2/S1]]) 356 return TM

357

358 -class ToneHole:

359 """ This class implements the transmission matrix of tone holes 360 References: 361 1. J. Dalmont et al., 362 "Experimental Determination of the Equivalent Circuit of an 363 Open Side Hole: Linear and Non Linear Behaviour," 364 Acta Acustica united with Acustica, vol. 88, pp. 567-575, 2002. 365 366 2. D.H. Keefe, 367 "Theory of the single woodwind tonehole," 368 J. Acoust. Soc. Am., vol. 72, pp. 676-687, 1982. 369 370 3. D.H. Keefe, 371 "Experiments on the single woodwind tonehole," 372 J. Acoust. Soc. Am., vol. 72, pp. 688-699, 1982. 373 374 4. van Walstijn, M., & Campbell, M. (2003). 375 "Discrete-time modeling of woodwind instrument bores using wave 376 variables," 377 J. Acoust. Soc. Am., vol. 113(1), 575-585. 378 379 5. Nederveen, C. J. (1998). 380 "Acoustical Aspects of Woodwind Instruments.", revised edition. 381 147 pages. Northern Illinois University Press. 382 383 """

384 - def __init__(self,d,t,p,d_pad=0.,h_pad=0.):

385 """ Initialization of a tonehole data structure 386 @param d: diameter of the tonehole 387 @param t: height of the tonehole (minimum) 388 @param p: position of the center of the tonehole with respect to 389 the beginning of the instrument. 390 """ 391 self.b = d/2. 392 self.t = t 393 self.p = p 394 395 self.r_pad = d_pad/2. 396 self.h_pad = h_pad 397 398 # the radiation condition for holes should be investigated 399 # as an approximation, use either Unflanged or Flanged open end 400 401 #self.ToneholeOpenEnd = Waveguides.UnflangedOpenEnd(self.b) 402 self.ToneholeOpenEnd = Radiation.FlangedOpenEnd(self.b)

403

404 - def __str__(self):

405 return "ToneHole(d=%.4f,t=%.4f,p=%.4f,d_pad=%.4f,h_pad=%.4f)"%\ 406 (2*self.b,self.t,self.p,2*self.r_pad,self.h_pad)

407

408 - def relocate(self,new_p):

409 """ Change the position of the tonehole. 410 @param new_p: new position of the center of the tonehole with respect 411 to the beginning of the instrument. 412 """ 413 self.p = new_p

414

415 - def resize(self,new_d):

416 """ Change the diameter of the tonehole. 417 @param new_d: new diameter of the tonehole 418 """ 419 self.b = new_d/2. 420 self.ToneholeOpenEnd = Radiation.FlangedOpenEnd(self.b)

421

422 - def get_position(self):

423 return self.p

424

425 - def delta(self,a):

426 """ Ratio of tone hole radius to body radius 427 a: body radius 428 429 The tone hole radius was set at initialization 430 """ 431 return self.b/a

432

433 - def area(self):

434 """ Cross sectional area of the tonehole """ 435 return numpy.pi * self.b**2

436

437 - def InnerLengthCorrection(self,delta):

438 " equation 4 from [1] " 439 return self.b * ( 0.82 - 1.4*delta**2 + 0.75*pow(delta, 2.7) )

440

441 - def MatchingVolumeLengthCorrection(self,delta):

442 """ correction for matching volume, valid for cylindical bores 443 ([1] eq. 6) 444 """ 445 return (self.b * delta / 8) * (1 + 0.207*delta**3)

446

447 - def RadiationLengthCorrection(self):

448 """ from [5] 449 """ 450 #Rc : outer radius of the hole (flange) 451 #rd : extra length correction due to the key hanging above the hole 452 453 if self.r_pad < self.b: 454 td = 0. 455 else: 456 td = self.PadLengthCorrection() 457 458 Rc = self.b # zero thickness hole... 459 return self.b*(0.831 - 0.135*(self.b/Rc) -0.073*(self.b/Rc)**4) + td

460 461 # we might ass eq. 9 Dalmont [1], for infinite flange 462

463 - def PadLengthCorrection(self):

464 " Nederveen (1998) eq. 38.6 " 465 466 # when there is no pad, radius is set to zero 467 if self.r_pad == 0.: 468 return 0. 469 470 # the range of validity of this equation is 471 if not (1. <= (self.r_pad/self.b) <= 4.): 472 print "Warning, the r_pad/b ratio is out of range: %.2f"%\ 473 (self.r_pad/self.b) 474 475 tp = 0.613*self.b*((self.r_pad/self.b)**0.18 * \ 476 (self.r_pad/self.h_pad)**0.39 - 1) 477 478 if tp < 0.: tp = 0. 479 480 return tp

481 482 # def ApproximateOpenShuntImpedance(self,f,a,T): 483 # air = Air(T) 484 # te_o = self.t + self.MatchingVolumeLengthCorrection(self.delta(a)) + \ 485 #self.RadiationLengthCorrection() 486 # M = air.rho*te_o/self.area() 487 # Zs_o = 1j*2*pi*f*M 488 # return Zs_o 489 # def ApproximateClosedShuntImpedance(self,f,a,T): 490 # air = Air(T) 491 # te_c = self.t + self.MatchingVolumeLengthCorrection(self.delta(a)) 492 # C = self.area()*te_c / (air.rho*air.c**2) 493 # Zs_c = 1./(1j*2*pi*f*C) 494 # return Zs_c 495

496 - def OpenShuntImpedance(self,f,a,T):

497 """ 498 formulation with losses, etc. 499 reference [1] 500 """ 501 air = air_properties(T) 502 k = 2*numpy.pi*f/air['c'] 503 Gamma = PropagationConstant(f,self.b,air,1) 504 505 L = self.t + self.MatchingVolumeLengthCorrection(self.delta(a)) +\ 506 self.PadLengthCorrection() 507 508 #print "t = %.4f, tm = %.4f, tp = %.4f, ti = %.4f"%\ 509 #(self.t,self.MatchingVolumeLengthCorrection(self.delta(a)), 510 #self.PadLengthCorrection(), 511 #self.InnerLengthCorrection(self.delta(a))) 512 513 Zl = self.ToneholeOpenEnd(f,T) 514 Zo = (Zl * cosh(Gamma*L) + sinh(Gamma*L)) / \ 515 (Zl * sinh(Gamma*L) + cosh(Gamma*L)) 516 Zi = 1j * k * self.InnerLengthCorrection(self.delta(a)) 517 Zc_hole = air['rho'] * air['c'] / self.area() 518 return Zc_hole*(Zi+Zo)

519

520 - def ClosedShuntImpedance(self,f,a,T):

521 """ impedance of a closed tone hole """ 522 air = air_properties(T) 523 Gamma = PropagationConstant(f,self.b,air,1) 524 L = self.t + self.MatchingVolumeLengthCorrection(a) 525 # should we take into account the impedance of the pad, 526 # which is made of lether ? 527 # we are currently assuming it has infinite impedance 528 529 Zc_hole = air['rho'] * air['c'] / self.area() 530 Zclosed = Zc_hole / tanh(Gamma*L) 531 return Zclosed

532

533 - def SeriesLengthCorrection(self,delta):

534 """ inertance 535 536 see ref [1] eq. 2, 3 and 4 537 """ 538 #t_a = -0.28 * self.b * delta**2 # eq. 2 539 #t_a = -(0.37 - 0.087*delta) * self.b * delta**2 # eq.3 540 541 t_a = -self.b * delta**2 * (2.72+0.540*delta+0.285*delta**2)**-1 542 # eq. 36 ref [4] 543 544 return t_a

545 546

547 - def SeriesImpedance(self,f,a,T):

548 """ Series impedance of the tone hole 549 """ 550 air = air_properties(T) 551 k = 2*numpy.pi*f/air['c'] 552 Zc_hole = air['rho'] * air['c'] / self.area() 553 return -1j * Zc_hole * k * self.SeriesLengthCorrection(self.delta(a))

554

555 - def TM(self,f,a,T,state=0):

556 """ Transmission matrix of the tone hole 557 558 f: frequency array 559 a: radius of the bore at tonehole location p 560 T: temperature 561 state: open or closed 562 """ 563 564 assert a > self.b, "Bore smaller than tonehole (a = %.4f, b = %.4f)"%\ 565 (a,self.b) 566 567 Za = self.SeriesImpedance(f,a,T) # same for open and closed hole 568 if state == 1: # open 569 Zs = self.OpenShuntImpedance(f,a,T) 570 else: # closed 571 Zs = self.ClosedShuntImpedance(f,a,T) 572 573 574 # note: in our notation, using P and ZU as the two variables, 575 # the matrices are made 576 # of adimensional impedances (normalized with the characteristic 577 # impedance Zc) 578 air = air_properties(T) 579 Zc_tube = air['rho'] * air['c'] / (numpy.pi * a**2) 580 TM = numpy.array([ [ ones(size(f)), (Za/Zc_tube)*ones(size(f))], 581 [Zc_tube/Zs, ones(size(f))] ]) 582 return TM

583

Source Code for Module wiat.TM.Waveguides