Acoustic Radiation Impedance

Download C# source code (version 3)

Github project

Introduction

I think I know what you are thinking: what is this about, and is it interesting to me? Well, probably more than you think, since they touch on many aspects of effective computational methods. If you are an acoustician, the answers you can bring out could give you a vital insight into the behavior of an acoustic wave as it hits an opening.

As the screenshot shows you, this is a normalized impedance that will give you a reactance loss, or, in fact, sometimes, a gain due to the energy that hits the surface compared to a normal incident wave. The impedance will also come with a phase displacement due to the oscillating mass, and here an end correction is often applied at low frequency. This is due to the fact that the oscillating mass is not only over the designated area but also in its surroundings.

For this project, I will rely on some mathematics functions that are found in Math.NET in the NuGet package. There are special functions that are missing from that package that I will supplement with my own.

Background: Circular baffle

The mathematical development of this started off well over 100 years ago, with the help of some of the earliest adopters and innovators in mathematical physics. Hermann von Helmholtz showed that the sound pressure from a vibrating surface was given by the two integrals:

The importance of sound mathematical reasoning by applying various conservation laws and logic shows how precise this kind of knowledge can be. Remember that in 1900, no electrical machinery excisted to do measuring, so most of these logics were derived theoretically with the help of the most basic experiments. They did some calculations and made logical conclusions from them.

Given the velocity, we can find the pressure that contributes to the first integral. The second integral requires that the pressure over the whole surface be known, and this cannot be known in general. But when placed in an infinite baffle, a heavy wall that extends indefinitely to each side and is adjacent to the vibrating surface, this integral vanishes. Along now comes one of the other great acousticians from this time period, Lord Rayleigh. He derives the following integral:

For a circle, this will give

This can be directly used to calculate the radiation impedance:

The code is simply:

/// <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);
}

Next, we need a good implementation of the Struve function, called StruveH.

Efficient Implementation of the Struve Function

I have come to the realization that this section could be a whole article by itself. It touches on a subject that I probably should be writing an article about: the asymptotic expansion. I decided to just do a good enough implementation with 9 significant digit accuracy from the range z = 0 to z = 20 000. It was also valid above this range, but that was not so interesting for the acoustical implementation. This search for accuracy also led me to implement tests to check that the functions did what they were supposed to do. I simply used Wolfram Alpha to generate the correct values that I used to check against.

Most implementations of the Struve function shared over the web seem to be way too complicated. Given that you already have access to a good implementation of the Neumann function, the implementation is indeed very simple. But unfortunately, it is not very accurate, and for lower orders of the Struve function, the worse it gets.

First, for the small input argument, the implementation with a Taylor series is more than good enough. The formula looks like the one below:

The implementation, which I did:

/// <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());
}

I will here use the formula from "The asymptotics of the Struve function H_ν(z) for large complex order and argument " by R. B. Paris. I start off with the work of R. B. Dingle:

This formula is only valid for $|arg \, z| < \frac{\pi}{2}$. However, it could easily be extended by the use of the relation:

This integral could be calculated with the Gaussian-Lagerre quadrature, since the integral is already of this form. I did this in Matlab but only got asymptotic convergence at very high values (from z = 10 000). So, instead, this integral is expanded with a reverse power series. All this is based on the work and articles by R. B. Dingle, and you can even take a look at some of his work in the book (see the references in the Paris article). The book is out of print but can be downloaded from Professor Sir Michael Victor Berry, FRS home page. The formula Dingle comes up with looks like this:

Now, we exchange the integral with a reverse power series, where 10 of the coefficients are given in Paris article.

/// <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
    };

}

Suppose you have a good enough (professional) implementation of the Neumann function. Then the Taylor series, together with the asymptotic expansion at $x \to \infty$ will give you a very straight-forward implementation that is accurate by 9 digits in the range 0–20 000 (and way above, but I did not test it as I did not need it).

/// <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;
}

Putting it all together:

/// <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);
}

This implementation will be valid for nearly any range that you could possibly want in practice. MathNet does have an implementation of the modified Struve function, which is a bit unfortunate since you can express the modified Struve function with the complex argument of the Struve function with some phase shift and negative imaginary number. Written out in a mathematical formula:

Elliptic piston

A formula for the radiation impedance of an elliptic piston, also derived by Mechel, looks like this:

The integral is dependent en the "elliptic function" B:

The general method for solving this is to expand the integrand into a power series and iterate. For the Series from the BesselJ1 funtion:

For the imaginary part containing the StruveH1 function

There is a problem with these methods, and that has to do with the rapidly growing powers of the alternating series, so I will eventually struggle if the values of x gets high enough. Like always, we need an asymptotic expansion when $x \to \infty $.

I know I have a very good algorithm for the circular baffle, so I will simply use that for the asymptotic expansion. I simply add the two different normalized impedances with the two different radiuses from the elliptic disk and take the average of them. This trick only works well enough when the natural wave number coincides with the average diameter of the elliptic disk, and this roughly happens when the normalized real part of the elliptic piston becomes 1. Ideally, we should check if the iterations become unruly and switch then, but that will have to be for a later time, and that is, in any case, not guaranteed to give you an accurate result.

        /// <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;
            }
        }

Field-Excited Rectangular radiator

The following derivation follows Mechel's approach with some of my own tweaking to the final implementation. Given a rectangular element placed in a hard baffle wall with an obliquely incident plane wave moving towards the surface, the term field excited refers to a vibration pattern set in motion by a wave propagating towards the opening.

The object here is a rectangle A with side lengths a and b on a hard baffle wall. The recatangle is excited by a plan wave with a polar incident angle of $\theta$ and an azimuthal angle of $\phi$. The average velocity pattern on the surface is:

Following Mechel's calculations, we have the fact that $|v(x,y)| = \text{constant}$ on the surface $A$ and that the radiation impedance follows from the average field impedance. The Fourier transform of the velocity distribution leads to a form with a double integration:

using $\alpha_1 = k_1/k_0$ and $\alpha_2 = k_2/k_0$:

This could be translated into:

The integral is not really suited for numerical integration as the limits are $k_0 \cdot a $ and $k_0 \cdot b$, which can be very large values indeed. Therefore, it is better to perform a substitution of variables and translate the double integral into:

The intermediate integral is given by:

With the variables:

The intermediate integral could be integrated by substituting the cosine terms with the identity:

Since every term is multiplied by $\cos(x)$ we can take the constant $1/4$ out of the integral, which can be a little simplified:

We now introduce a variable called $c_n$ which has the following 4 forms:

Given that $n=1,2,3,4$ the general term for $c_n$:

There is also the possibility to combine the original values and simplify the expression for $c_n$:

However, the implementation is much simpler in code with the use of the bitwise and operator:

// 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;

Resulting in:

Using standard integration and integration by parts, these are easily solved:

Sorting out and rearranging gives the following:

The code that generates this integral can be implemented as such:

/// <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;
    
}

The two integrals can now easily be implemented by a Cauchy sum. Which is the same as the Riemann sum, only that a Cauchy sum uses finitely many elements while Riemann sets the number of elements to infinity.

// 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));
    }
}

If you are tempted to use some Newton-Coates-style summation, like Gaussian quadrature, you should remember that the requirement for using this is a realistically well-behaved function. There is no guarantee of that here, at least not in advance. That is why I opted for a more safe, straight-forward summation. This code also seems to be very stable even for massively large openings, like 100 times 100 meters.

Rectangular Piston

There is a direct implementation for a Rectangular Piston as well, and it is included as a reference. But don't use it, instead use the field excited rectangular with an incoming wave which is normal incident theta = 0 and phi = 0. This is much more stable. The code just so we are crystal clear here:

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));
}

Wide and Narrow strips

With the updated Struve function and the use of a rectangular piston with width 10 times longer than the height, these are now also accurate for a very large range.

End corrections

The imaginary part of the radiation impedance really tells you about the mass that is oscillating in the opening. At lower frequencies, the mass is actually larger than the mass only contained in the opening, so we have something called end corrections. These are dependent on the shape of the opening and the available surrounding area that can oscillate, and they must be added on for a more realistic imaginary impedance. 

Further work

There are many more implementations of different shapes of pistons, like spheres, cylinders, doughnuts, etc. There is also the case where several small perforated openings work in tandem. But for those, I will leave it to you to investigate the references for further study.

References

• Formulas of Acoustics, 2^nd ed. by F. P. Mechel et. al.
• Building Acoustics by Tor Erik Vigran

I also used some material for articles that are referenced by those books, so I did not include them here.

I would also like to thank the late F. P. Mechel for his kind help and for sending his own code used in Formulas of Acoustics. 

History

• 14^th January, 2024: Initial version
• 17^th January, 2024: Updated explanations and code
• 25^th January, 2024: Updated Struve to 9-digit accuracy and added tests to mathematical functions.