0

我正在尝试重现本文https://hal.archives-ouvertes.fr/file/index/docid/683477/filename/clarinette-logique-8.pdf中描述的模型。

这是一个返回半径为 a 和长度为 L 的圆柱体的传递矩阵的方法。

# density of air
rho = 1.2250

# viscosity of air
eta = 18.5

# transfer matrix for a cylinder
def Tcyl(a, L):
    rv = a * math.sqrt(rho*omega/eta)
    Z0 = (rho * constants.c) / (constants.pi * a**2)
    Zc = Z0 * complex(1 + 0.369/rv, -(0.369/rv + 1.149/rv**2))
    gamma = k * complex(1.045/rv + 1.080/rv**2, 1 + 1.045/rv)
    return
       np.matrix([[cmath.cosh(gamma*L),Zc*cmath.sinh(gamma*L)],
       [(1.0/Zc)*cmath.sinh(gamma*L),cmath.cosh(gamma*L)]], dtype=complex)

不同频率的输入阻抗计算如下。

for f in range(1,4000):
    
    # angular frequency
    omega = 2*constants.pi*f

    # wave number
    k = omega/constants.c

    # transfer matrix for a cylinder of radius 7.5mm and length 100mm
    T = Tcyl(7.5, 100.0)

    # radiation impedance
    Zr = complex(0,0)

    # [P,U] vector (pression and velocity)
    T = np.dot(T,np.matrix([[Zr],[complex(1,0)]], dtype=complex))

    # input impedance (Z = P/U)
    Z = T.item(0)/T.item(1)

播放频率满足等式Im[Z]=0。绘制 ZI 的虚部时，得到下图：阻抗错误

这显然是错误的，因为预期的输出应该是这样的：正确的阻抗

我究竟做错了什么？谢谢你。

4

1 回答 1

0

你有

    Z0 = (rho * constants.c) / (constants.pi * a**2)

单簧管的阻抗取决于声音的速度，而不是光的速度。替换constants.c为343，您的结果将更接近。我仍然不确定它是否完全正确，但更接近。

作为一名单簧管演奏家，我试图让人们认为我的手指移动得像光一样快，但事实并非如此。

于 2021-03-21T03:10:10.980 回答