Numerical Solutions to Vibrations of Bars

• Download source - 14.1 MB

A Word of Warning

Introduction

The solution to the differential wave equation leads to a combination of two main functions with two different amplitudes, i.e., A*Sin and B*Cos. The boundary conditions imposed usually result in solving the equation for the (usually scaler) variables A and B. Typically, A and B are visualized as a standing incoming wave and leaving wave after D'Alembert's principle solution to the transverse wave equation. The problem is that the boundary conditions are often complex in nature, and there is no simple way of solving these equations, so generally, a numerical method is required. The question is just which one to use.

I will, for the most part, skip the derivations of how the equations are set up and deal with the solution straight away. All the solutions are solved using the continued fraction expansion and truncating it. This will lead to a Pade approximation to the function we are trying to solve. The denominator is removed by some simple algebraic manipulation such that a single complex polynomial root-finding algorithm can be used. In this case, that is the Jenkins-Traub algorithm. The implementation in pseudo-code:

1. Replace function like:

   \(\sin(x)\)

   with a fraction:

   \( \sin(x) = \frac{a_0 + a_1 \cdot x + a_2 \cdot x^2 + \cdots + a_n \cdot x^n}{b_0 + b_1 \cdot x + b_2 \cdot x^2 + \cdots + b_n \cdot x^n }\)

2. Rearrange it into a polynomial, i.e., remove all fractions and set the equation equal to zero

3. Solve the truncated polynomial with a root-finding technique.

The reason that the Pade approximation is so useful in acoustics is due to the solutions to the wave equation. They are often truncated alternating series, and it is in these cases that the Pade approximation or Shanks transform is most useful. Due to the alternating power series, the solutions tend either to infinity or negative infinity very fast due to the limitations of computational numbers.

For the transverse wave equations, a different approach is taken. The starting point (as always, by the way. Including the Pade approximation) is the Taylor series. Since the function here is two know Taylor series, they are just symbolically multiplied together, and the root algorithm is set to find the roots from this.  

The physical difference between the longitudinal wave and the transverse wave can be illustrated quite simply. Assume that someone jumps at either one of the two locations:

The transfer of energy into the brown beam (bar) from jumps at the longitudinal section will result in compressions in the vertical direction. While jumping at the transverse location will result in a transverse wave along the grey beam (bar). The astute might realize that the words actually describe the difference between the direction of the initial force compared to the wave direction. The longitudinal wave moves in the same direction as the force is applied. The transverse wave travels at a 90-degree angle of the applied force. While jumping up and down, the wave travels from side to side in the gray beam.

Yes I am simplifying a bit, since jumping in the longitudinal section will also result in a transverse wave, but the explanation is hopefully simpler to understand if I just mention this afterward. In solid structures, they almost always coincide, but not in fluids. Here longitudinal waves are the only possible ones (well, mostly. You can have special situations where they do play a primary role).

Longitudinal Wave Equation

Free Mass Loaded Bar

The equation to solve for is this:

The solution using Linq and the Symbolic package:

private void xCotxSolver(int N = 11)
{
    // Continued fraction for x*cot(x)
    List<string> ContFraction_xCotx = new List<string>();
    ContFraction_xCotx.Add("1");
    for (int i = 1; i <= N; i++)
    {
        ContFraction_xCotx.Add("-x*x");
        ContFraction_xCotx.Add((2 * (i + 1) - 1).ToString());
    }

    // Solves
    // tan(k*L) = - (m/m_b)*k*L =>
    // kl*cot(kl) = -(m_b/m);
    // kl*cot(kl)  + (m_b/m) == 0;

    double m_b = 1;
    double m = 1;

    ContinuedFraction xCotx = new ContinuedFraction(ContFraction_xCotx);
    Expr equation = xCotx.Numerator + (m_b / m) * xCotx.Denomiator;
    List<Complex> roots = RealPolynomialRoots.Roots
                          (equation.ReverseRealCoefficients("x"));

    List<double> result =
        roots.
        // Only interested in positive values
        FindAll(x => (x.Real >= 0)).
        // And only real solutions
        Select(x => x.Real).
        // Make it to a List
        ToList<double>();

    // Sort low to high k*L (2*pi*f/c_0)*L solutions
    result.Sort();
}

The solutions should be \( kL = 2.03, 4.91, 7.98, \cdots \).

The Freely Vibrating Bar - General Boundary Condition

The equation that one seeks the solution to:

The code that does this:

private void FreeVibrationBar_GeneralBC()
      {
          // Chapter 3.6 page 76
          double rho_L = 1;
          double C = 1;
          Complex Zm0 = new Complex(10, 1);
          Complex Zm = new Complex(10, 1);

          Complex u = new Complex(0, 1);
          u *= (Zm0 / (rho_L * C) + (Zm / (rho_L * C))) /
               (1 + (Zm0 / (rho_L * C)) * (Zm / (rho_L * C)));

          // Continued fraction for Tan(x)
          List<string> ContFraction_Tan_x = new List<string>();
          ContFraction_Tan_x.Add("0");
          ContFraction_Tan_x.Add("x");
          ContFraction_Tan_x.Add("1");

          // Make sure k is odd, 25 good enough for 15 roots approximately depending on u
          for (int i = 2; i < 25; i++)
          {
              ContFraction_Tan_x.Add("-x*x");
              ContFraction_Tan_x.Add((2 * (i) - 1).ToString());
          }
          ContinuedFraction Tanx = new ContinuedFraction(ContFraction_Tan_x);
          var TanXnumerator = Tanx.Numerator.Coefficients(Expr.Variable("x"));
          var TanXDenominator = Tanx.Denomiator.Coefficients(Expr.Variable("x"));

          var PolynomialAtanExpansion = (Tanx.Numerator - Tanx.Denomiator *
              new MathNet.Numerics.Complex32((float)u.Real, (float)u.Imaginary)).
              Coefficients(Expr.Variable("x")).
              ToList().
              Select(x => x.ComplexNumberValue);

          List<Complex> roots_Atanx = ComplexPolynomialRoots.Roots
                                      (PolynomialAtanExpansion.Reverse<Complex>().ToArray());
          List<Complex> Positive_roots_Atanx =
          roots_Atanx.
          // Only interested in positive values
          FindAll(xxr => (xxr.Real >= 0)).
          // Make it to a List
          Distinct<Complex>().ToList();
          // To ignore erroneous roots
          //.Find(x => Complex.Tan(x) - u).ToList();

          List<Complex> Atan = new List<Complex>();
          int mm = 0;
          foreach (var item in Positive_roots_Atanx)
          {
              Atan.Add(item);
              Atan.Add(Math.PI * mm + Complex.Atan(u));
              //tantan.Add(Complex.Tan(tantan[tantan.Count-1]) - u);
              mm += 1;
          }
      }

Forced Vibrations on a Bar

The equation we seek the solution to is:

By making the substitution, some algebra can be performed:

private void ForcedVibrationOnABar()
       {
           // Chapter 3.7 page 76-77 in Kinser & Fray
           // Continued fraction for Tan(x)
           List<string> ContFraction_Tan_x = new List<string>();
           ContFraction_Tan_x.Add("0");
           ContFraction_Tan_x.Add("x");
           ContFraction_Tan_x.Add("1");
           // Make sure k is odd, good enough for 15 roots approximately
           for (int i = 2; i < 25; i++)
           {
               ContFraction_Tan_x.Add("-x*x");
               ContFraction_Tan_x.Add((2 * (i) - 1).ToString());
           }
           ContinuedFraction Tanx = new ContinuedFraction(ContFraction_Tan_x);

           // Z_m/rho_L*c = r + j*x

           double x1 = 10;
           double r1 = 0;

           //Trying to solve the equation
           // x*tan^2(k*L) + (r^2+x^2-1)*tan(k*L) - x = 0

           var er1 = Expr.Parse(x1.ToString()).Multiply
                     (Tanx.Numerator.Multiply(Tanx.Numerator));

           double val = r1 * r1 + x1 * x1 - 1;

           var er2 = Expr.Parse((val).ToString()).Multiply
                     (Tanx.Numerator).Multiply(Tanx.Denomiator);
           var er3 = Expr.Parse(x1.ToString()).Multiply(Tanx.Denomiator.Pow(2));

           var tyty = er1
               .Add(er2)
               .Subtract(er3)
               // Group together equal powers
               .RationalSimplify(Expr.Variable("x"));

           var TanX = tyty
               // Get the coefficients
               .Coefficients(Expr.Variable("x"))
               // Get number in front of coefficient
               .Select(x => x.RealNumberValue)
               .ToList<double>()
               // Reverse array from ascending to descending powers
               .Reverse<double>()
               .ToArray<double>();
           List<Complex> roots_tanxx = RealPolynomialRoots.Roots(TanX);
           var tw = roots_tanxx.
               Select(x => x.Real).
               ToList<double>().
               FindAll(x => x > 0).
               Distinct<double>().
               ToList<double>();
           tw.Sort();

           // Solution check
           List<double> tw2 = new List<double>();
           for (int uuu = 0; uuu < tw.Count(); uuu++)
           {
               tw2.Add(Math.PI * uuu + 0.099669);
               tw2.Add(Math.PI * uuu - 1.4711);
           }

           tw2 = tw2.FindAll(x => x > 0).Distinct().ToList();
           tw2.Sort();

           List<double> test = new List<double>();
           for (int hg = 0; hg < tw.Count; hg++)
           {
               test.Add(tw[hg]);
               test.Add(tw2[hg]);
           }
       }

Transverse Wave Equation

The general solution to the transverse wave equation for a bar:

In order to solve this equation, we need its derivative, second derivative and third derivative:

Second derivative:

Third derivative:

These equations are used together with the appropriate boundary conditions:

Fixed end:

Free end:

Simply supported end:

Boundary conditions (BC) for the clamped end. This is usually the fixed beam since y = 0 and the derivative is also zero. This means that the initial position does not change at all, and the beam is always at a straight angle from its fixed point regardless of the vibrations present in the rest of the beam. The involves solving the equation:

This is almost equal to the equation with a vibrating bar free at both ends. Think of this as the frequency that is generated by throwing a bar at the floor.

With these two equations, the Pade approximation is not used (although you could). Instead, the Taylor (or rather the MacLaurin) series is used directly. 

 private void ClamedFreeBarVibration()
        {
            // Solutions for Bar clamed at one end and bar free on both ends
            // Clamped in one end gives Cos(x)*Cosh(x) = -1
            // Free at both ends gives Cos(x)*Cosh(x) = 1
            Expr xr = Expr.Parse("x");
            Expr one = Expr.Parse("1");
            Expr cosh = Expr.Parse("1");
            Expr cos = Expr.Parse("1");

            // Cant really go higher than 25 due to set limit in root finding tool
            for (int t = 1; t <= 25; t++)
            {
                cosh += xr.Pow(2 * t) / Factorial(2 * t);
                cos += Math.Pow(-1, t) * xr.Pow(2 * t) / Factorial(2 * t);
            }

            // Alternative formulation: 1/4*(exp(-1i*x) + exp(i*x))*(exp/-x)+exp(x)) + 1 == 0

            // Change .Add to .Subtract in order to switch equation conditions
            Expr coshcos = cosh.Multiply(cos).Add(one).RationalSimplify(xr);
            var Coefficients = coshcos
                .Coefficients(Expr.Variable("x"))
                .Select(x => x.RealNumberValue)
                .ToList<double>()
                .Reverse<double>()
                .ToList<double>();

            List<Complex> roots_coshcos = RealPolynomialRoots.Roots(Coefficients.ToArray());

            List<Double> roots_cch =
                roots_coshcos.
                // Only interested in positive values where imaginary part is zero               
                Where(p => ((p.Imaginary <= 0.01) && (p.Real >= 0))).
                // Change from Complex to Double 
                Select(x => x.Real).
           
                // Make it to a distinct list
                Distinct<double>().
                ToList<double>().
                // Roots that dont fit due to truncation 
                // (Could do Newton-Rapson or other root polishing techniques)
                FindAll(x => Math.Abs(coshcos.Evaluate
                (new Dictionary<string, FloatingPoint> { { "x", x } }).RealValue) < 1e-4).
                // To frequencies
                Select(x => Math.Round(x * 2 / Math.PI, 3)).
                Distinct<double>().
                ToList<double>();

            // Sort low to high k*L (2*pi*f/c_0)*L solutions
            roots_cch.Sort();

        }

The simply supported-free bar’s transcendental equation is:

Reference

Fundamentals of Acoustics, 4^th edition, Kinsler, Frey, Coppers and Sanders.

History

• 20^th March, 2021: Initial release
• 24^th March, 2021: Corrected written formula (code correct) for freely vibrating bar, added referances and fixed some minor formatting issues