15,030,908 members
Articles / Programming Languages / C#
Article
Posted 19 Mar 2021

3.8K views
4 bookmarked

Numerical Solutions to Vibrations of Bars

Rate me:
Numerically solves the problem of vibrations on bars with various boundary conditions using Jenkins-Traub algorithm
The acoustic problem of vibrations of bars is usually solved with graphic methods and most books that have these equations in them just mention that you could use the Newton-Raphson method to solve it. The approach in this article is different. This article uses the NuGet symbolic mathematics package MathNet.Symbolics together with my C# implementation of the Jenkin-Traub root-finding algorithm of polynomials to solve these problems.

A Word of Warning

I might as well warn you from the start. If you don't know the acoustic equations that I'm solving here, this article is not going to be very useful for you. Or you might be a mathematician that has a similar math problem that requires solving an equation using a polynomial root finding method. For now, it's just a short article, and I might upgrade it later with some nice solution plots.

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

The equation to solve for is this:

$$$\tan(k L) = - \frac{m}{m_b} k L$$$

The solution using Linq and the Symbolic package:

C#
private void xCotxSolver(int N = 11)
{
// Continued fraction for x*cot(x)
List<string> ContFraction_xCotx = new List<string>();
for (int i = 1; i <= N; i++)
{
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:

$$$\tan(k L) = j\frac{\left( Z_{m0} / \left(\rho_L c\right) \right) + \left( Z_m / \left( \rho_L c \right) \right)}{1+\left( Z_{m0} / \left(\rho_L c\right) \right)\left( Z_m / \left( \rho_L c \right) \right)}$$$

The code that does this:

C#
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>();

// Make sure k is odd, 25 good enough for 15 roots approximately depending on u
for (int i = 2; i < 25; i++)
{
}
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)
{
mm += 1;
}
}


Forced Vibrations on a Bar

The equation we seek the solution to is:

$$$Z_{m0} = \rho_L c \frac{\left( Z_m / \rho_L c \right) + j \tan(k L)}{1 +\left( Z_m / \rho_L c \right) j \tan(k L) }$$$

By making the substitution, some algebra can be performed:

C#
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>();
// Make sure k is odd, good enough for 15 roots approximately
for (int i = 2; i < 25; i++)
{
}
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
.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 = tw2.FindAll(x => x > 0).Distinct().ToList();
tw2.Sort();

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


Transverse Wave Equation

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

$$$y(x,t) = \left[ A \cdot \cosh(g \cdot x) + B \cdot \sinh(g \cdot x) + C \cdot \cos(g \cdot x ) + D \cdot \sin(g \cdot x )\right] \cos (\omega t + \phi )$$$

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

$$$\frac{d \, y(x,t)}{dx} = g \cdot \left[ A \cdot \sinh(g \cdot x) + B \cdot \cosh(g \cdot x) - C \cdot \sin(g \cdot x) + D \cdot \cos(g \cdot x) \right] \cos(\omega t + \phi )$$$

Second derivative:

$$$\frac{d^2 \, y(x,t)}{{dx}^2} = g^2 \cdot \left[A \cdot \cosh(g \cdot x) + B \cdot \sinh(g \cdot x) - C \cdot \cos(g \cdot x) - D \cdot \sin(g \cdot x) \right] \cos (\omega t + \phi )$$$

Third derivative:

$$$\frac{d^3 \, y(x,t)}{{dx}^3} =g^3 \cdot \left[A \cdot \sinh(g \cdot x) + B \cdot \cosh(g \cdot x) + C \cdot \sin(g \cdot x) - D \cdot \cos(g \cdot x) \right] \cos (\omega t + \phi )$$$

These equations are used together with the appropriate boundary conditions:

Fixed end:

$$$y = 0 \qquad \text{and} \qquad \frac{\partial y}{\partial x} = 0$$$

Free end:

$$$\frac{\partial^2 y}{{\partial x}^2} = 0 \qquad \text{and} \qquad \frac{\partial^3 y}{{\partial x}^3} = 0$$$

Simply supported end:

$$$y = 0 \qquad \text{and} \qquad \frac{\partial^2 y}{{\partial x}^2} = 0$$$

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:

$$$\cos(g \cdot L) \cosh(g \cdot L) = -1$$$

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.

$$$\cos(g \cdot L) \cosh(g \cdot L) = 1$$$

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

C#
 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
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:

$$$\frac{\tan(g L)} {\tanh(gL)} = 1$$$

Reference

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

History

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

Share

 Engineer Norway
No Biography provided

 First Prev Next
 Is there an easier way? Member 1172068122-Mar-21 7:38 Member 11720681 22-Mar-21 7:38
 Re: Is there an easier way? Kenneth Haugland22-Mar-21 14:37 Kenneth Haugland 22-Mar-21 14:37
 Re: Is there an easier way? Member 1172068123-Mar-21 10:23 Member 11720681 23-Mar-21 10:23
 Re: Is there an easier way? Kenneth Haugland23-Mar-21 11:00 Kenneth Haugland 23-Mar-21 11:00
 Re: Is there an easier way? Member 1172068123-Mar-21 11:47 Member 11720681 23-Mar-21 11:47
 What is the application of your algorithm? Mirzakhmet Syzdykov21-Mar-21 23:38 Mirzakhmet Syzdykov 21-Mar-21 23:38
 Re: What is the application of your algorithm? Kenneth Haugland22-Mar-21 14:34 Kenneth Haugland 22-Mar-21 14:34
 Last Visit: 31-Dec-99 18:00     Last Update: 18-Sep-21 9:26 Refresh 1