声辐射阻抗

下载 C# 源代码（第 3 版）

Github 项目

引言

我想我知道你在想什么：这有什么关系，它对我感兴趣吗？嗯，可能比你想象的要多，因为它们涉及有效计算方法的许多方面。如果你是一名声学家，你得出的答案可以让你深入了解声波撞击开口时的行为。

正如截图所示，这是一个归一化阻抗，它会给你一个电抗损耗，或者实际上有时会因为撞击表面的能量相对于垂直入射波而获得增益。阻抗还将伴随着由振荡质量引起的相位位移，在这里，对于低频通常会应用端部修正。这是因为振荡质量不仅在指定区域内，而且在其周围区域内。

对于这个项目，我将依赖 NuGet 包中 Math.NET 的一些数学函数。该包缺少一些特殊函数，我将用我自己的函数来补充。

背景：圆形挡板

这个数学发展始于 100 多年前，得到了一些早期数学物理学研究者和创新者的帮助。赫尔曼·冯·亥姆霍兹（Hermann von Helmholtz）表明，振动表面的声压由两个积分给出

通过应用各种守恒定律和逻辑进行严谨的数学推理的重要性表明了这种知识的精确性。请记住，在 1900 年，还没有电气机械用于测量，因此这些逻辑大多是理论推导出来的，并辅以最基本的实验。他们进行了一些计算并从中得出了逻辑结论。

给定速度，我们可以找到对第一个积分有贡献的压力。第二个积分要求已知整个表面的压力，而这通常是未知的。但是，当放置在无限挡板中时，一个向两侧无限延伸并与振动表面相邻的重墙，这个积分就会消失。现在，让我们来看看那个时期的另一位伟大的声学家，瑞利勋爵（Lord Rayleigh）。他推导出了以下积分

对于圆形，这将给出

这可以直接用于计算辐射阻抗

代码很简单

        /// <summary>
        /// Calculates radiation impedance for a Circular piston placed in an infinite baffel.
        /// </summary>
        /// <param name="k_0">Wavenumber 2*pi*f/c_0 of surrounding propagation medium</param>
        /// <param name="a">radius of circle</param>
        /// <param name="Z_0">Specific impedance of surrounding propagation medium. If not specified it returns normalized impedance</param>
        /// <returns></returns>
        public static Complex CircularBaffle(double k_0, double a, double Z_0 = 1)
        { 
            if (k_0 == 0)    
                return Complex.Zero; 
            
            double k0a = k_0 * a;
            Complex i = new Complex(0, 1);
            return Z_0*(1 - MathNet.Numerics.SpecialFunctions.BesselJ(1,2*k0a)/k0a + i * SpecialFunctions.StruveH(1,2*k0a)/k0a);
        }

接下来，我们需要一个好的 Struve 函数实现，称为 StruveH。

Struve 函数的高效实现

我意识到这一部分本身就可以成为一篇完整的文章。它触及了我可能应该写一篇文章的主题：渐近展开。我决定用 9 位有效数字精度在 z = 0 到 z = 20,000 的范围内实现一个足够好的实现。它在这个范围之上也有效，但对于声学实现来说并不那么有趣。对精度的追求也促使我实现了测试来检查函数是否按预期工作。我直接使用 Wolfram Alpha 生成了正确的值，并用它们进行了检查。

网络上分享的大多数 Struve 函数实现似乎都过于复杂。鉴于您已经拥有了一个好的 Neumann 函数实现，该实现确实非常简单。但不幸的是，它并不十分准确，对于 Struve 函数的低阶，情况会更糟。

首先，对于小的输入参数，使用泰勒级数实现已经绰绰有余。公式如下

我实现的这个

/// <summary>
/// A straight-forward implementation of the Struve function valid for Abs(z) < 16
/// </summary>
/// <param name="n">Order</param>
/// <param name="z">Argument</param>
/// <returns>Struve function: H_n(z) by Taylor series expansion</returns>
private static Complex StruveHTaylor(double n, Complex z)
{
    // Termwise result
    Complex TermResult = Complex.Zero;
    // For epsilon algorthm
    List<Complex> SummationTerms = new List<Complex>();
    
    // If zero just return that
    if (z == Complex.Zero) { return Complex.Zero; }
    
    // Cap the number of iterations for the loop
    double MaxIteration = 25d;
    
    // Precalculate a value that does not change
    Complex Log_z2 = Complex.Log(z * 0.5d);
    Complex iPi = new Complex(0, Math.PI);
    
    // Estimated error
    Complex error;
    
    // Accepted tolerance               
    double tol = 1e-10;
    
    // Standard Taylor series implementation except for taking the logaritmic values instead
    for (double m = 0; m < MaxIteration; m += 1d)
    {
        Complex LogarithmicTermEvaluation =
            // This is just i*pi*m since Log(-1) = i*pi
            m * iPi
            // Use precalcualted value
            + (2 * m + n + 1) * Log_z2
            // Natural logarithm of the gamma function
            - MathNet.Numerics.SpecialFunctions.GammaLn(m + 1.5d)
            - MathNet.Numerics.SpecialFunctions.GammaLn(m + n + 1.5d);
            
        // The exponential will remove the logarithm
        TermResult += Complex.Exp(LogarithmicTermEvaluation);
        
        // Termwise results
        SummationTerms.Add(TermResult);
        
        // Using the Epilon algorithm generally seems to shave
        // off one or two iteration for the same precition
        // Should be a good algorithm to use
        // since its an alternating Taylor series
        Complex Summation = EpsilonAlgorithm(SummationTerms.ToArray());
        
        // Must have at least two values to use Wynns epsilon algorithm
        if (m > 0)
        {
            // Assume that the Epsilon algortim improves the summation by at least order of 1.
            // So the error is estimated by simple substraction
            error = SummationTerms.Last() - Summation;
            
            // Compare magnitude of error with accepted tolerance
            if (tol > Complex.Abs(error))
            {
                //Debug.WriteLine("Number of iterations: " + m.ToString() +
                //" and error " + Complex.Abs(error).ToString("N10") + " with Wynn: ");
                return Summation;
            }
        }
    }
    
    // Tolerance not reached within maximum iteration
    return EpsilonAlgorithm(SummationTerms.ToArray());
}

我将在此使用 R. B. Paris 的论文《"The asymptotics of the Struve function H_ν(z) for large complex order and argument "》中的公式。我从 R. B. Dingle 的工作中开始

此公式仅适用于 $|arg \, z| < \frac{\pi}{2}$。但是， 通过使用关系

这个积分可以用高斯-拉盖尔求积法计算，因为积分已经是这种形式了。我在 Matlab 中做到了这一点，但只在非常高的值（z = 10,000 以上）得到了渐近收敛。因此，取而代之的是，这个积分通过反幂级数展开。所有这些都基于 R. B. Dingle 的工作和文章，你甚至可以看看他的一些著作（见 Paris 文章中的参考文献）。这本书已经绝版，但可以从迈克尔·贝瑞教授 (F.R.S.) 的主页下载。Dingle 得出的公式如下

现在，我们将积分替换为反幂级数，其中 10 个系数在 Paris 文章中给出。

        /// <summary>
        /// Steepest decent evaluation values 
        /// </summary>
        /// <param name="q">The value n/Z is used for the asymptotic expansion</param>
        /// <returns>All first 10 expansions</returns>
        private static Complex[] C_k(Complex q)
        {
            // I find the formulas easier to read by pre-calculating the powers
            Complex q2 = q * q;
            Complex q3 = q2 * q;
            Complex q4 = q3 * q;
            Complex q5 = q4 * q;
            Complex q6 = q5 * q;
            Complex q7 = q6 * q;
            Complex q8 = q7 * q;
            Complex q9 = q8 * q;
            Complex q10 = q9 * q;

            return new Complex[]{
                1d,
                2d * q,
                6d * q2 - 0.5d,
                20d * q3 - 4d * q,
                70d * q4 - 45d / 2d * q2 + 3d / 8d,
                252d * q5 - 112d * q3 + 23d / 4d * q,
                924d * q6 - 525d * q4 + 301d / 6d * q2 - 5d / 16d,
                3432d * q7 - 2376d * q5 + 345d * q3 - 22d / 3d * q,
                12870d * q8 - 21021d / 2d * q6 + 16665d / 8d * q4 - 1425d / 16d * q2 + 35d / 128d,
                48620d * q9 - 45760 * q7 + 139139d / 12d * q5 - 1595d / 2d * q3 + 563d / 64d * q,
                184756d * q10 - 196911d * q8 +61061d * q6 - 287289d / 48d * q4 + 133529d / 960d * q2 - 63d / 256d
            };

        }

假设您有一个足够好（专业）的 Neumann 函数实现。那么泰勒级数，连同 $x \to \infty$ 时的渐近展开，将为您提供一个非常直接的实现，在 0-20,000 的范围内（甚至更高，但我没有测试，因为我不需要）具有 9 位数字的精度。

        /// <summary>
        /// Struve Neumann function that returns H_n(z) - Y_n(z) by asymptotic series
        /// </summary>
        /// <param name="n">Order</param>
        /// <param name="z">Argument</param>
        /// <returns> H_n(z) - Y_n(z). Uses steepest decent evaluation </returns>
        public static Complex StruveHnNeumannYn(double n, Complex z)
        {
            Complex Result = Complex.Zero;

            // This term is actually multiplied with the summation term. I just want to take the logaritm and 
            // try to cancel out all high order terms 
            Complex ConstantTerm = (n - 1) * Complex.Log(0.5 * z)
                            - 0.5 * Complex.Log(Math.PI)
                            - MathNet.Numerics.SpecialFunctions.GammaLn(n + 0.5);

            // The coefficients from H_n(z) => n/Z
            Complex[] Ck = C_k(n / z);

            for (int k = 0; k < Ck.Length; k++)
            {
                // Trying to avoid massive devisions
                Complex item = ConstantTerm + Complex.Log(Ck[k])
                    + MathNet.Numerics.SpecialFunctions.GammaLn(k + 1)
                    - k * Complex.Log(z);

                // Remove logarithm and add the asymptotic series
                Result += Complex.Exp(item);
            }

            return Result;
        }

整合所有内容

/// <summary>
/// Struve function
/// </summary>
/// <param name="n">Order</param>
/// <param name="z">Argument</param>
/// <returns>H_n(z)</returns>
public static Complex StruveH(double n, Complex z)
{
    if (Complex.Abs(z) < 20)
        // Taylor series
        return StruveHTaylor(n, z);
    else
        // Asymtotic expansion
        return StruveHnNeumannYn(n, z) + MathNet.Numerics.SpecialFunctions.BesselY(n, z);
}

这个实现对于您在实践中可能需要的几乎任何范围都是有效的。MathNet 确实有一个修改后的 Struve 函数实现，这有点可惜，因为您可以用修改后的 Struve 函数以及带有相位偏移和负虚数的复数参数来表示它。用数学公式写出来是

椭圆活塞

Mechel 也推导出了椭圆活塞的辐射阻抗公式，如下所示

积分取决于“椭圆函数”B

解决这个问题的通用方法是将积分展开成幂级数并进行迭代。对于 BesselJ1 函数的级数

对于包含 StruveH1 函数的虚部

这些方法存在一个问题，那就是交错级数快速增长的幂，因此如果 x 的值足够大，我最终会陷入困境。一如既往，我们需要在 $x \to \infty$ 时的渐近展开。

我知道对于圆形挡板我有一个非常好的算法，所以我将直接将其用于渐近展开。我只需将两个不同的归一化阻抗与椭圆盘的两个不同半径相加，然后取它们的平均值。当自然波数与椭圆盘的平均直径重合时，此技巧才足够好，而这大约发生在椭圆活塞的归一化实部变为 1 时。理想情况下，我们应该检查迭代是否变得难以处理，然后切换，但这将是以后要做的事情，而且无论如何都不能保证能得到准确的结果。

        /// <summary>
        /// Calculates radiation impedance for a Elliptic piston placed in an infinite baffel.
        /// </summary>
        /// <param name="k">Wavenumber 2*pi*f/c_0 of surrounding propagation medium</param>
        /// <param name="a">Long radius of ellipse</param>
        /// <param name="b">Short radius of ellipse</param>
        /// <param name="Z_0">Specific impedance of surrounding propagation medium. If not specified it returns normalized impedance</param>
        /// <returns></returns>
        public static Complex EllipticBaffle(double k, double a, double b, double Z_0 = 1)
        {
            if (k == 0)
                return Complex.Zero;

            double ka = k * a;
            double ka2 = ka * ka;
            double beta = Math.Min(a, b) / Math.Max(a, b);

            double Limit = Math.Ceiling(4 * ka) + 4d;
            double MaxLimit = 80;

            if (Limit < MaxLimit)
            {
                //Taylor series
                double Z_r = beta * (EllipticSumZ_r(ka2, beta, Limit));
                double Z_i = 4 / Math.PI / Math.PI * ka * beta * EllipticSumZ_i(ka2, beta, Limit);
                return Z_0 * new Complex(Z_r, Z_i);
            }
            else
            {
                // Asymptotic expansion
                Complex LongSideBaffel = CircularBaffle(k, a);
                Complex ShortSideBaffel = CircularBaffle(k, b);
                return Z_0 * (LongSideBaffel + ShortSideBaffel) / 2;
            }
        }

场激发矩形辐射器

接下来的推导遵循 Mechel 的方法，并加入了我自己对最终实现的修改。给定一个放置在硬挡板墙上的矩形元件，并有一个斜入射的平面波朝向表面传播，场激发的术语指的是由传播到开口的波激发的振动模式。

对象是一个硬挡板墙上的矩形 A，其边长为 a 和 b。该矩形由具有极入射角 $\theta$ 和方位角 $\phi$ 的平面波激发。表面上的平均速度模式为

根据 Mechel 的计算，我们发现 $|v(x,y)| = \text{constant}$ 在表面 A 上，并且辐射阻抗源于平均场阻抗。速度分布的傅里叶变换导致一个双重积分的形式

使用 $\alpha_1 = k_1/k_0$ 和 $\alpha_2 = k_2/k_0$

这可以转化为

积分不适合数值积分，因为积分限是 $k_0 \cdot a $ 和 $k_0 \cdot b$，它们可能非常大的值。因此，最好进行变量替换并将双重积分转换为

中间积分由下式给出

使用变量

可以通过替换余弦项来对中间积分进行积分，使用恒等式

由于每个项都乘以 $ \cos(x) $，我们可以将常数 $1/4$ 从积分中提取出来，这可以稍微简化

我们现在引入一个称为 $c_n$ 的变量，它有以下 4 种形式

给定 $n=1,2,3,4$， $c_n$ 的一般项

还可以组合原始值并简化 $c_n$ 的表达式

然而，在代码中使用按位与运算符可以使实现更加简单

// Generate sequence 1, -1, 1, -1
double mod1 = (1 & n) == 0 ? 1 : -1;

// Generate sequence 1, 1, -1, -1
double mod2 = (2 & n) == 0 ? 1 : -1;

Complex C_n = mod1 * alpha(j) + mod2 * beta(j) - 1;

结果是

使用标准积分和分部积分，这些很容易求解

整理和重排得到以下结果

生成此积分的代码可以这样实现

/// <summary>
/// Intermediate integral
/// </summary>
/// <param name="R">Upper integration limit</param>
/// <param name="C_n">Permutation group of 4 elements</param>
/// <param name="U">Constant depending size</param>
/// <param name="V">Constant depending on angle, size and wave number</param>
/// <param name="W">Constant dependent on angle only</param>
/// <returns></returns>
public static Complex I_R(Complex R, Complex C_n, Complex U, Complex V, Complex W)
{
    Complex C_n2 = C_n * C_n;
    Complex C_n3 = C_n2 * C_n;
    Complex i = new Complex(0, 1);
    Complex expR = Complex.Exp(i * C_n * R);
    
    Complex result =
        (   expR * (    - U * i * C_n2 
                        + V * (C_n - i * C_n2 * R) 
                        + W * (-i * C_n2 * R * R + 2 * C_n * R + 2 * i)
                    )
            +(  U * i * C_n2 
                - V * C_n 
                - W * 2 * i ) 
        ) / (C_n3 * 4);
        
    return result;
    
}

两个积分现在可以通过柯西求和（Cauchy sum）轻松实现。这与黎曼和（Riemann sum）相同，只是柯西求和使用有限个元素，而黎曼设定元素数量为无穷大。

// Integral one
for (double i = 0; i < arctan; i += deltaPhi)
{
    for (int n = 0; n < 4; n++)
    {
        // Generate sequence 1, -1, 1, -1
        double mod1 = (1 & n) == 0 ? 1 : -1;
        
        // Generate sequence 1, 1, -1, -1
        double mod2 = (2 & n) == 0 ? 1 : -1;
        
        Complex C_n = mod1 * alpha(i) + mod2 * beta(i) - 1;
        result += I_R(k0a / Math.Cos(i), C_n, U, V(i), W(i));
    }
}

// Integration two
for (double j = arctan; j < Math.PI / 2; j += deltaPhi)
{
    for (int n = 0; n < 4; n++)
    {
        // Generate sequence 1, -1, 1, -1
        double mod1 = (1 & n) == 0 ? 1 : -1;
        
        // Generate sequence 1, 1, -1, -1
        double mod2 = (2 & n) == 0 ? 1 : -1;
        
        Complex C_n = mod1 * alpha(j) + mod2 * beta(j) - 1;
        result += I_R(k0b / Math.Sin(j), C_n, U, V(j), W(j));
    }
}

如果您倾向于使用类似高斯求积法的牛顿-科茨式求和，您应该记住，使用它的要求是函数行为现实上良好。这里没有这样的保证，至少不是事先的。这就是为什么我选择了更安全、更直接的求和方法。即使对于像 100x100 米那样规模巨大的开口，这段代码似乎也非常稳定。

矩形活塞

矩形活塞也有直接的实现，并且包含在内作为参考。但不要使用它，而是使用场激发的矩形活塞，并使用垂直入射的入射波（theta = 0 和 phi = 0）。这会更稳定。为了清晰起见，代码如下

            List<Complex> RetangularPistonResult = new List<Complex>();
            for (double k = 0; k <= UpperMode; k += deltaWaveNumber)
            {
                RetangularPistonResult.Add(Acoustics.Radiation.FieldExcited.Rectangular(k, 0, 0, a, b));
                //   RetangularPistonResult.Add(Acoustics.Radiation.Pistons.RetangularBaffel(k, a, b));
            }

宽窄条带

通过更新的 Struve 函数和使用宽度是高度 10 倍的矩形活塞，这些现在对于非常大的范围来说也是准确的。

端部修正

辐射阻抗的虚部实际上告诉我们开口中振荡的质量。在较低频率下，质量实际上大于仅包含在开口中的质量，因此我们有称为端部修正的东西。它们取决于开口的形状和可以振荡的周围可用区域，并且必须添加它们以获得更真实的虚阻抗。

进一步的工作

还有许多其他形状活塞的实现，如球体、圆柱体、甜甜圈等。也有几个小型穿孔开口协同工作的案例。但对于这些，我将留给您查阅参考文献以作进一步研究。

参考文献

• F. P. Mechel 等人的《声学公式》第二版。
• Tor Erik Vigran 的《建筑声学》。

我还使用了一些这些书籍引用的文章材料，因此此处未包含它们。

我还要感谢已故的 F. P. Mechel 先生的慷慨帮助，并感谢他提供他在《声学公式》中使用过的代码。

历史

• 2024 年 1 月 14 日：初始版本
• 2024 年 1 月 17 日：更新了解释和代码
• 2024 年 1 月 25 日：Struve 更新至 9 位精度，并为数学函数添加了测试。