wiat.MM.Solver

1 # Copyright 2008 Antoine Lefebvre 2 # 3 # This file is part of "The Wind Instrument Acoustic Toolkit". 4 # 5 # "The Wind Instrument Acoustic Toolkit" is free software: 6 # you can redistribute it and/or modify it under the terms of 7 # the GNU General Public License as published by the 8 # Free Software Foundation, either version 3 of the License, or 9 # (at your option) any later version. 10 # 11 # "The Wind Instrument Acoustic Toolkit" is distributed in the hope 12 # that it will be useful, but WITHOUT ANY WARRANTY; without even 13 # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 14 # See the GNU General Public License for more details. 15 # 16 # You should have received a copy of the GNU General Public License 17 # along with "The Wind Instrument Acoustic Toolkit". 18 # If not, see <http://www.gnu.org/licenses/>. 19 # 20 # All rights reserved. 21 22 """ Solver.py - Multimodal wave propagation solver. 23 24 This is an implementation of the theory presented in references [1] and [2] 25 with boundary layer loss model from reference [3]. 26 27 References: 28 29 1. V. Pagneux, N. Amir, and J. Kergomard, 30 "A study of wave propagation in varying cross-section waveguides by 31 modal decomposition. Part I. Theory and validation," 32 J. Acoust. Soc. Am., vol. 100, 1996, pp. 2034-2048. 33 34 2. Amir, N., Pagneux, V., & Kergomard, J. (1997). 35 A study of wave propagation in varying cross-section waveguides by 36 modal decomposition. Part II. Results. 37 The Journal of the Acoustical Society of America, 101(5), 2504-2517. 38 39 3. Bruneau, A. M., Bruneau, M., Herzog, P., & Kergomard, J. (1987). 40 "Boundary layer attenuation of higher order modes in waveguides." 41 Journal of Sound and Vibration, 119(1), 15-27. 42 43 44 Example:: 45 46 >>> import wiat 47 >>> din = 0.02 48 >>> dout = 0.025 49 >>> f = numpy.arange(50.,500.,10.) 50 >>> air = wiat.util.air_properties(25.) 51 >>> c = wiat.MM.Waveguides.QuarterCircleHorn(din, dout) 52 >>> rad = wiat.MM.Radiation.RadiationFlanged(dout) 53 >>> Z5 = wiat.MM.Solver.SolveImpedance(c,rad,f,air,5) 54 55 """ 56 57 __all__ = ['SolveImpedance', 'SolveAdmittance'] 58 59 import time 60 61 import numpy 62 from numpy import sin, cosh, sqrt, real, imag, pi, dot, transpose 63 import scipy.special 64 from scipy.special import j0,j1,jn_zeros,jvp 65 66 from wiat.util import air_properties 67 from Functions import flat2matrix, matrix2flat 68

69 -def calculate_gamma(N):

70 return numpy.concatenate(([0.],scipy.special.jn_zeros(1,N)))

71

72 -def _Q(N):

73 """ Q is the constant matrix defined in reference 1, equation 28 74 75 @param N: number of modes 76 """ 77 gamma = calculate_gamma(N) 78 Q = [] 79 for n in range(N): 80 Q.append([]) 81 gamman2 = gamma[n]**2 82 for m in range(N): 83 if m==n: 84 Q[n].append(0.) 85 else: 86 gammam2 = gamma[m]**2 87 Q[n].append(gammam2 / (gammam2 - gamman2)) 88 return numpy.array(Q)

89

90 -def _K(N,a,k,air):

91 """ Multimodal propagation constant 92 (with boundary layer thermo-viscous losses) 93 94 See [2], equations 5 and 6 as well [3] 95 96 """ 97 gamma = calculate_gamma(N) 98 lv = air['mu']/(air['rho']*air['c']) 99 lt = lv/air['nu']**2 100 epsilonv = complex(1.,1.) * sqrt(k/2) * sqrt(lv) 101 epsilont = complex(1.,1.) * sqrt(k/2) * (air['gamma']-1)*sqrt(lt) 102 K = [] 103 for i in range(N): 104 alpha = gamma[i]/a 105 epsiloni = (1. - alpha**2/k**2)*epsilonv + epsilont 106 K.append( k**2 - alpha**2 + (2j*k/a)* 107 (numpy.imag(epsiloni) -1j*numpy.real(epsiloni))) 108 return K

109

110 -def _Zprime(Zflat,z,geometry,k,air,N,Q):

111 """ Matrical Riccati equation for the generalized impedance matrix Z 112 113 This is an implementation of equation 32 of reference [1], which is 114 the same as equation 1 of reference [2]. 115 116 """ 117 Z = flat2matrix(Zflat,N) 118 Zc = air['rho']*air['c']/geometry.S(z) 119 K = numpy.diag(_K(N,geometry.radius(z),k,air), 0) 120 tmp = dot(K,Z) 121 tmp = dot(Z,tmp)/(k*Zc) 122 tmp2 = numpy.diag(numpy.ones((N))*k*Zc, 0) 123 Zp = 1j*(-tmp2+tmp)+geometry.Sp_S(z)*(dot(Q,Z) + dot(Z,transpose(Q))) 124 return matrix2flat(Zp)

125 126

127 -def SolveImpedance(geometry,Zr,freqs,air,N):

128 freqs = numpy.asarray(freqs).flatten() 129 130 zf = geometry.get_length() 131 zs = numpy.array([0., zf]) 132 Z00 = [] 133 Zc = air['rho']*air['c']/geometry.S(zf) 134 Zcr = air['rho']*air['c']/geometry.S(0.) 135 print "Integrating system for %d modes and %d frequencies..."%(N,len(freqs)) 136 initialtime = time.clock() 137 Q = _Q(N) 138 for f in freqs: 139 140 # when frequency is 0, it is a flow problem and the method can't work 141 # we append a null impedance in that case 142 if f == 0.: 143 Z00.append(complex(0.,0.)) 144 continue 145 146 previoustime = time.clock() 147 k = 2*numpy.pi*f/air['c'] 148 Zflat,infodict = scipy.integrate.odeint( 149 _Zprime, Zcr*Zr(N,k), zs, (geometry,k,air,N,Q), 150 None, full_output=1, mxstep=15000, 151 rtol=1e-5, atol=1e-10) 152 Z = flat2matrix(Zflat[1],N) 153 Z00.append(Z[0][0]/Zc) 154 print "% 7.2f %.2f seconds"%(f, time.clock() - previoustime) 155 print "...done in %.2f seconds"%(time.clock() - initialtime) 156 return numpy.array(Z00)

157

158 -def _Yprime(Yflat,z,geometry,k,air,N,Q):

159 """ Riccati equation for admittance 160 161 Implementation of equation 4 of reference [2]. 162 """ 163 Y = flat2matrix(Yflat, N) 164 kZc = k * air['rho'] * air['c'] / geometry.S(z) 165 K = numpy.diag(_K(N,geometry.radius(z),k,air), 0) 166 Yp = -1j*K/(kZc) - geometry.Sp_S(z)*(dot(transpose(Q),Y) + dot(Y,Q)) +\ 167 1j*kZc*dot(Y,Y) 168 return matrix2flat(Yp)

169

170 -def SolveAdmittance(geometry,Yr,freqs,air,N):

171 freqs = numpy.asarray(freqs).flatten() 172 173 zf = geometry.L 174 zs = numpy.array([0., zf]) 175 Y00 = [] 176 Yc = geometry.S(zf)/(air['rho']*air['c']) 177 print "Integrating system for %d modes and %d frequencies..."%(N,len(freqs)) 178 initialtime = time.clock() 179 Q = _Q(N) 180 for f in freqs: 181 previoustime = time.clock() 182 k = 2*numpy.pi*f/air['c'] 183 Yflat,infodict = scipy.integrate.odeint( 184 _Yprime, Yr(N,k),zs,(geometry,k,air,N,Q), 185 None,full_output=1,mxstep=15000, 186 rtol=1e-5,atol=1e-10) 187 Y = flat2matrix(Yflat[1],N) 188 Y00.append(Y[0][0]/Yc) 189 print "% 7.2f %.2f seconds"%(f, time.clock() - previoustime) 190 print "...done in %.2f seconds"%(time.clock() - initialtime) 191 return numpy.array(Y00)

192 193 # to be done... 194 #def CalculateMaxima(Z,f,solver,inst,rad,number_of_maxima=None,**kwargs): 195 # derivative = numpy.diff(numpy.abs(Z)) 196 # signarray = numpy.sign(derivative) 197 # crossing = numpy.diff(signarray) 198 # indices = crossing.nonzero() 199 # maximalist = [] 200 # for i in indices[0]: 201 # # do not take into account the minimum 202 # if crossing[i] > 0: 203 # continue 204 # if len(maximalist) == number_of_maxima: 205 # break 206 # f0 = f[i-1] 207 # f1 = f[i+1] 208 # if kwargs.get('Y'): 209 # results = scipy.optimize.brent(lambda f: numpy.abs(solver.Ysolve(f,inst,rad)),brack=(f0,f1),full_output=True) 210 # else: 211 # results = scipy.optimize.brent(lambda f: 1./numpy.abs(solver.Zsolve(f,inst,rad)),brack=(f0,f1),full_output=True) 212 # fres = results[0][0] 213 # zres = 1./results[1][0] 214 # #fres = scipy.optimize.brentq(lambda f: numpy.imag(solver.Zsolve(numpy.array([f]),inst,rad)),f0,f1) 215 # if zres > 1.: # it is a maxima and not a minima 216 # maximalist.append((fres, zres)) 217 # 218 # return maximalist 219

Source Code for Module wiat.MM.Solver