Difference between revisions of "Waves on a Variable Beam"

From WikiWaves
Jump to navigationJump to search
 
(33 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 
== Introduction ==
 
== Introduction ==
We consider here the equations for a non-uniform free beam. We begin by presenting the theory for a uniform beam.
+
We consider here the equations for a non-uniform free beam. We begin by presenting the theory for a uniform beam.<br\>
 +
 
 +
An example of the motion for a non-uniform beam is demonstrated below:
 +
[[Image:variable_beam_ii.gif|right|Vibration of variable beam]]
  
 
== Equations for a beam ==
 
== Equations for a beam ==
 
{{equations for a beam}}
 
{{equations for a beam}}
  
== Solution for a uniform beam in [[Eigenfunctions for a Uniform Free Beam|eigenfunctions]] ==
+
== Solution for a uniform beam in [[Eigenfunctions for a Uniform Free Beam|eigenfunctions]] with no pressure ==
 
 
If the beam is uniform the equations can be written
 
<center>
 
<math> D \frac{\partial^{4}\zeta}{\partial x^{4}}+m\frac{\partial^{2}\zeta}{\partial t^{2}}=0</math> </center>
 
We can express the deflection as the series
 
<center><math>  \zeta(x,t)=\sum_{n=0}^{\infty} \left( A_n w_n(x) \cos(\lambda_n t) + 
 
A_n w_n(x) \frac{\cos(\lambda_n t)}{\lambda_n} \right) </math></center>
 
where <math>w_n</math> are the [[Eigenfunctions for a Uniform Free Beam]] and <math>\lambda_n</math>
 
  
If we introduce the initial conditions
+
{{solution for a uniform beam in eigenfunctions}}
<center>
 
:<math>  \zeta(x,0)=f(x) \,\! </math>
 
:<math>  \frac{\partial \zeta(x,0)}{\partial t}=g(x)  </math></center>
 
Then <math>  A_n \,\!</math> and <math> B_n \,\!</math> can be found using orthogonality properties:
 
<center>
 
:<math>  A_n=\frac{\int_{-L}^{L}f(x)w_n(x)\mathrm{d}x}{\int_{-L}^{L}X_n(x)X_n(x)\mathrm{d}x} \,\! </math>
 
</center>
 
<center>
 
:<math>  B_n=\frac{\int_{-L}^{L}g(x)w_n(x)\mathrm{d}x}{\int_{-L}^{L}X_n(x)X_n(x)\mathrm{d}x}  </math></center>
 
Note that <br />
 
- these modes drop off very quickly (ie <math> v_4 \,\!</math> oscillates about zero with negligible amplitude), so the higher order vibration modes can be ignored.<br />
 
- As time progresses (<math> t \rightarrow \infty \,\!</math>), each mode will vibrate around the zero displacement line with frequency <math> \overline{\omega}_{n}\,\!</math>.<br />
 
- Having obtained eigenvalues <math> k_n \,\!</math>, the natural frequencies can easily be obtained <math>\overline{\omega}_{n}=k_{n}^{2}\sqrt{\frac{D}{m}}\,\!</math>.
 
<br />
 
  
 
== Variational Techniques ==
 
== Variational Techniques ==
As an alternative to solving the eigenvalue problem (\ref{first_eig_prob}), we can equivalently minimise the Rayleigh Quotient \cite{linton}:
+
As an alternative to solving the eigenvalue problem, we can equivalently minimise the [http://en.wikipedia.org/wiki/Rayleigh_quotient Rayleigh Quotient] ([[Linton and McIver 2001]]):
<center><math>  R[v]=\frac{\varepsilon[v]}{H[v]}  </math></center>
+
<center><math>  R[\zeta]=\frac{\varepsilon[\zeta]}{H[\zeta]}  </math></center>
where <math> \varepsilon[v] \,\!</math> is known as the energy functional, and <math> H[v] \,\!</math> is some constraint.
+
where <math> \varepsilon[\zeta] \,\!</math> is known as the energy functional, and <math> H[\zeta] \,\!</math> is some constraint.
 
<br /><br />
 
<br /><br />
In the case of a uniform beam with zero  transverse load, our energy functional is \cite{lanczos}:
+
In the case of a uniform beam with zero  transverse load, our energy functional is ([[Lanczos 1949]]):
<center><math>  \varepsilon[v]=\frac{D}{2}\int_{-L}^{L}\left(\frac{\partial^2 v}{\partial x^2} \right)^2 \mathrm{d}x    </math></center>
+
<center><math>  \varepsilon[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 v}{\partial x^2} \right)^2 \mathrm{d}x    </math></center>
 
If we choose
 
If we choose
<center><math>  H[v]=\int_{-L}^{L}m v^2 \mathrm{d}x=1  </math></center>
+
<center><math>  H[\zeta]=\int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1  </math></center>
 
Then minimizing the Rayleigh Quotient can be expressed as a variational problem <br />
 
Then minimizing the Rayleigh Quotient can be expressed as a variational problem <br />
 
<center>
 
<center>
min <math>  R[v]=\frac{D}{2}\int_{-L}^{L}\left(\frac{\partial^2 v}{\partial x^2} \right)^2 \mathrm{d}x  </math>
+
min <math>  R[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 \zeta}{\partial x^2} \right)^2 \mathrm{d}x  </math>
subject to <math> \int_{-L}^{L}m v^2 \mathrm{d}x=1  </math></center>
+
subject to <math> \int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1  </math></center>
</center>
+
which will in turn, also solve our eigenvalue problem.
which will in turn, also solve our eigenvalue problem (\ref{first_eig_prob}).
 
  
 
== Rayleigh-Ritz method ==
 
== Rayleigh-Ritz method ==
We can essentially replace the variational problem of finding a <math> y(x) \,\!</math> that extremises <math> J[y] \,\!</math> to finding a set of constants <math> a_1, a_2, ..., a_N \,\!</math> which extremise <math> J[a_1, a_2, ..., a_N] \,\!</math>, through using the Rayleigh-Ritz method.<br /><br />
+
We can essentially replace the variational problem of finding a <math> y(x) \,\!</math> that extremises <math> J[y] \,\!</math> to finding a set of constants <math> a_1, a_2, ..., a_N \,\!</math> which extremise <math> J[a_1, a_2, ..., a_N] \,\!</math>, through using the [http://en.wikipedia.org/wiki/Rayleigh-Ritz_method Rayleigh-Ritz method].<br /><br />
 
We solve
 
We solve
 
<center><math> \frac{\partial}{\partial a_k}J[a_1, a_2, ..., a_N]=0 \,\!</math></center>
 
<center><math> \frac{\partial}{\partial a_k}J[a_1, a_2, ..., a_N]=0 \,\!</math></center>
Line 55: Line 35:
 
<center><math> \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\!</math></center>
 
<center><math> \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\!</math></center>
  
== Non-uniform free beam ==
+
== Solution for a Non-uniform free beam in eigenfunctions ==
 
In the case of a non-uniform free beam, the Euler-Bernoulli beam equation becomes:
 
In the case of a non-uniform free beam, the Euler-Bernoulli beam equation becomes:
<center><math> \frac{\partial^2}{\partial x^2}\left( D(x) \frac{\partial^{2}v}{\partial x^{2}}\right) +m(x)\frac{\partial^{2}v}{\partial t^{2}}=0 \,\!</math></center>
+
<center><math> \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2}\zeta}{\partial x^{2}}\right) +\gamma(x)\frac{\partial^{2}\zeta}{\partial t^{2}}=0 \,\!</math></center>
Using separation of variables, where <math> v(x,t)=\widehat{X}(x)\widehat{T}(t) \,\!</math>, we obtain a corresponding eigenfunction problem (with eigenvalues denoted by <math> \mu_n=\widehat{k}_n^4 \,\!</math>):
+
Using separation of variables, where <math> \zeta(x,t)=\widehat{X}(x)\widehat{T}(t) \,\!</math>, we obtain a corresponding eigenfunction problem (with eigenvalues denoted by <math> \mu_n=\widehat{k}_n^4 \,\!</math>):
<center><math> \frac{\partial^2}{\partial x^2}\left( D(x) \frac{\partial^{2} \widehat{X}_n}{\partial x^{2}}\right)=m(x)\mu_n \widehat{X}_n \,\!</math></center>
+
<center><math> \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2} \widehat{X}_n}{\partial x^{2}}\right)=\gamma(x)\mu_n \widehat{X}_n \,\!</math></center>
 
which leads us to the following variational expression:
 
which leads us to the following variational expression:
 
<center>
 
<center>
min <math>  J[\widehat{X}_n]=\frac{1}{2}\int_{-L}^{L} \bigg\{D(x)(\widehat{X}_n^{''})^2 -\mu_n m(x) \widehat{X}_n^2 \bigg\}dx   </math>
+
min <math>  J[\widehat{X}_n]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)(\widehat{X}_n^{''})^2 -\mu_n \gamma(x) \widehat{X}_n^2 \bigg\}\mathrm{d}x   </math>
 
</center>
 
</center>
 
We can approximate this using Rayleigh-Ritz and obtain:
 
We can approximate this using Rayleigh-Ritz and obtain:
<center><math> J[\vec{a}]=\frac{1}{2}\int_{-L}^{L} \bigg\{D(x)[\sum_{i=1}^{N}a_{i}{X}_{i}^{''}]^2 -\mu_i m(x) [\sum_{i=1}^{N}a_{i}{X}_{i}]^2 \bigg\}dx \,\!</math></center>
+
<center><math> J[\vec{a}]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)[\sum_{i=1}^{N}a_{i}{X}_{i}^{''}]^2 -\mu_i \gamma(x) [\sum_{i=1}^{N}a_{i}{X}_{i}]^2 \bigg\}\mathrm{d}x \,\!</math></center>
 
Then extremising <math> J[\vec{a}] \,\!</math> as follows
 
Then extremising <math> J[\vec{a}] \,\!</math> as follows
 
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=0 \,\!</math></center>
 
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=0 \,\!</math></center>
 
for all <math> k=1,2,...,N \,\!</math>, allows us to obtain:
 
for all <math> k=1,2,...,N \,\!</math>, allows us to obtain:
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=\sum_{i=1}^{N} \left\{ \int_{-L}^{L}D(x) {X}_{i}^{''} {X}_{k}^{''}dx -\mu_i \int_{-L}^{L} m(x){X}_{i}{X}_{k}dx\right\}a_i=0 \,\!</math></center>
+
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=\sum_{i=1}^{N} \left\{ \int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x -\mu_i \int_{-L}^{L} \gamma(x){X}_{i}{X}_{k}\mathrm{d}x\right\}a_i=0 \,\!</math></center>
for all <math> k=1,2,...N \,\!</math> (here <math> \mu_i \,\!</math> denotes both the Lagrange multiplier and eigenvalues of a non-uniform beam)
+
for all <math> k=1,2,...N \,\!</math> (here <math> \mu_i \,\!</math> denotes both the [http://en.wikipedia.org/wiki/Lagrange_multipliers Lagrange multiplier] and eigenvalues of a non-uniform beam)
  
 
== Converting to Matrix Form ==
 
== Converting to Matrix Form ==
 
If we define
 
If we define
 
<center>
 
<center>
:<math> K_{ik}=\int_{-L}^{L}D(x) {X}_{i}^{''} {X}_{k}^{''}dx \,\!</math>
+
:<math> K_{ik}=\int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\!</math>
:<math> M_{ik}=\int_{-L}^{L} m(x){X}_{i}{X}_{k}dx \,\!</math>
+
:<math> M_{ik}=\int_{-L}^{L} \gamma(x){X}_{i}{X}_{k}\mathrm{d}x \,\!</math>
 
</center>
 
</center>
 
then for all <math> k=1,2,...N \,\!</math>:
 
then for all <math> k=1,2,...N \,\!</math>:
<center><math> \frac{\partial}{\partial a_k}J[\vec{a}]=\sum_{i=1}^{N}(K_{ik}-\mu_i M_{ik})a_i=0 \,\!</math></center>
+
<center><math> \frac{\partial}{\partial a_k}J[\mathbf{a}]=\sum_{i=1}^{N}(K_{ik}-\mu_i M_{ik})a_i=0 \,\!</math></center>
Let us create the matrix <math> K \,\!</math>, with elements <math> K_{ik} \,\!</math>, the matrix <math> M \,\!</math>, with elements <math> M_{ik} \,\!</math>, the sparse matrix <math> \Lambda \,\!</math> with <math> \mu_i \,\!</math> terms on the diagonal, and the vector <math> \vec{a} \,\!</math>.  We can consequently express the above sum in the following way:
+
Let us create the matrix <math> K \,\!</math>, with elements <math> K_{ik} \,\!</math>, the matrix <math> M \,\!</math>, with elements <math> M_{ik} \,\!</math>, the sparse matrix <math> \Lambda \,\!</math> with <math> \mu_i \,\!</math> terms on the diagonal, and the vector <math> \mathbf{a} \,\!</math>.  We can consequently express the above sum in the following way:
<center><math> (K-\Lambda M)\vec{a}=\vec{0} \,\!</math></center>
+
<center><math> (K-\Lambda M)\mathbf{a} = 0 \,\!</math></center>
So we can solve for <math> \vec{a} \,\!</math> in any of the problems below to obtain coefficients <math> a_i \,\!</math>:
+
So we can solve for <math> \mathbf{a} \,\!</math> in any of the problems below to obtain coefficients <math> a_i \,\!</math>:
 
<center>
 
<center>
:<math> (K-\Lambda M)\vec{a}=\vec{0} \,\!</math>
+
:<math> (K-\Lambda M)\mathbf{a} = 0  \,\!</math>
:<math> K\vec{a}=\Lambda M\vec{a}  \,\!</math>
+
:<math> K\mathbf{a}=\Lambda M\mathbf{a}  \,\!</math>
:<math> (M^{-1}K)\vec{a}=\Lambda \vec{a} \,\!</math>
+
:<math> (M^{-1}K)\mathbf{a}=\Lambda \mathbf{a} \,\!</math>
 
</center>
 
</center>
  
 
== Quadrature ==
 
== Quadrature ==
Rather than tediously evaluate an integral for each element of the matrices <math> K \,\!</math> and <math> M \,\!</math>, we can break up the integral into subintervals with weights <math> w_h \,\!</math> using Composite Simpson's Rule:
+
Rather than tediously evaluating an integral for each element of the matrices <math> K \,\!</math> and <math> M \,\!</math>, we can break up the integral into subintervals with weights <math> w_h \,\!</math> using [http://en.wikipedia.org/wiki/Composite_simpson%27s_rule#Composite_Simpson.27s_rule Composite Simpson's Rule]:
 
<center>
 
<center>
 
:<math> \int_a^b f(x) dx \approx \frac{h}{3} \left[ f(x_0)+2 \sum_{j=1}^{n/2-1}f(x_{2j})+4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n)\right] \,\!</math>
 
:<math> \int_a^b f(x) dx \approx \frac{h}{3} \left[ f(x_0)+2 \sum_{j=1}^{n/2-1}f(x_{2j})+4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n)\right] \,\!</math>
Line 97: Line 77:
 
So we can express <math> K_{ik} \,\!</math> as follows:
 
So we can express <math> K_{ik} \,\!</math> as follows:
 
<center>
 
<center>
:<math> K_{ik}=\int_{-L}^{L}D(x) {X}_{i}^{''} {X}_{k}^{''}dx \,\!</math>
+
:<math> K_{ik}=\int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\!</math>
:<math> \approx \sum D(x_h){X}_{i}^{''}(x_h) {X}_{k}^{''}(x_h)w_h \,\!</math>
+
:<math> \approx \sum \beta(x_h){X}_{i}^{''}(x_h) {X}_{k}^{''}(x_h)w_h \,\!</math>
 
:<math> \approx \vec{X}_{i}^{''} H \vec{X}_{k}^{''T}  \,\!</math>
 
:<math> \approx \vec{X}_{i}^{''} H \vec{X}_{k}^{''T}  \,\!</math>
 
</center>
 
</center>
Line 105: Line 85:
 
<math> H =
 
<math> H =
 
\begin{bmatrix}
 
\begin{bmatrix}
   D(x_1)w_1 & 0 & ... &0 \\
+
   \beta(x_1)w_1 & 0 & ... &0 \\
   0 & D(x_2)w_2 & ... &0 \\
+
   0 & \beta(x_2)w_2 & ... &0 \\
 
   \vdots &  & \ddots&\vdots \\
 
   \vdots &  & \ddots&\vdots \\
   0 & ... & ...&D(x_h)w_h
+
   0 & ... & ...&\beta(x_h)w_h
 
\end{bmatrix}
 
\end{bmatrix}
 
\,\!</math>
 
\,\!</math>
Line 132: Line 112:
 
<math> J =
 
<math> J =
 
\begin{bmatrix}
 
\begin{bmatrix}
m(x_1)w_1 & 0 & ... &0 \\
+
\gamma(x_1)w_1 & 0 & ... &0 \\
0 & m(x_2)w_2 & ... &0 \\
+
0 & \gamma(x_2)w_2 & ... &0 \\
 
\vdots &  & \ddots&\vdots \\
 
\vdots &  & \ddots&\vdots \\
0 & ... & ...&m(x_h)w_h
+
0 & ... & ...&\gamma(x_h)w_h
 
\end{bmatrix}
 
\end{bmatrix}
 
\,\!</math>
 
\,\!</math>
Line 174: Line 154:
 
<center><math> \sum_{n=0}^{\infty}A_n\widehat{X}_n(x)=f(x) \,\!</math></center>
 
<center><math> \sum_{n=0}^{\infty}A_n\widehat{X}_n(x)=f(x) \,\!</math></center>
 
multiplying by <math> m(x)\widehat{X}_n(x) \,\!</math> and integrating allows us to obtain:
 
multiplying by <math> m(x)\widehat{X}_n(x) \,\!</math> and integrating allows us to obtain:
<center><math> A_n=\frac{\int_{-L}^{L}f(x)m(x)\widehat{X}_n(x)dx}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)dx} \,\!</math></center>
+
<center><math> A_n=\frac{\int_{-L}^{L}f(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\!</math></center>
 
Then using the second initial condition, and applying the same technique yields
 
Then using the second initial condition, and applying the same technique yields
<center><math> B_n=\frac{1}{\widehat{k}_n^2}\frac{\int_{-L}^{L}g(x)m(x)\widehat{X}_n(x)dx}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)dx} \,\!</math></center>
+
<center><math> B_n=\frac{1}{\widehat{k}_n^2}\frac{\int_{-L}^{L}g(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\!</math></center>
 
Consequently,
 
Consequently,
 
<center><math> v(x,t)=A_0\widehat{X}_0(x)+A_1\widehat{X}_1(x)+\sum_{n=2}^{\infty} A_n\widehat{X}_n(x)\cos(\widehat{k}_n^2 t)+\sum_{n=2}^{\infty} B_n\widehat{X}_n(x) \sin(\widehat{k}_n^2 t)  \,\!</math></center>
 
<center><math> v(x,t)=A_0\widehat{X}_0(x)+A_1\widehat{X}_1(x)+\sum_{n=2}^{\infty} A_n\widehat{X}_n(x)\cos(\widehat{k}_n^2 t)+\sum_{n=2}^{\infty} B_n\widehat{X}_n(x) \sin(\widehat{k}_n^2 t)  \,\!</math></center>
 
where <math> \widehat{X}_0(x) \,\!</math> and <math> \widehat{X}_2(x) \,\!</math> are the two rigid modes of the non-uniform beam.  Note that <math> B_0 \,\!</math> and <math> B_1 \,\!</math> do not exist.
 
where <math> \widehat{X}_0(x) \,\!</math> and <math> \widehat{X}_2(x) \,\!</math> are the two rigid modes of the non-uniform beam.  Note that <math> B_0 \,\!</math> and <math> B_1 \,\!</math> do not exist.
  
 +
== Matlab Code ==
 +
 +
A program to solve for a variable beam can be found here {{variable beam}}
 +
 +
=== Additional code ===
 +
 +
This program requires
 +
* {{eigenvalues beam}}
 +
* {{eigenvectors beam}}
 +
* {{double derivative eigenvectors beam}}
 +
 +
[[Category:Simple Linear Waves]]
  
== To Do List ==
+
==References==
MATLAB image of subplots <br />
+
[[Lanczos 1949]],
MATLAB code wikiwaves, em_beam_dd.m,ev_beam,em_beam needed <br />
+
[[Linton and McIver 2001]],
References for authors/sources, as well as for formulas inside a page (ie refer to formula 33)<br />
+
[[Rao 1986]]
  
[[Category:Simple Waves]]
+
[[Category:Complete Pages]]

Latest revision as of 23:45, 2 July 2009

Introduction

We consider here the equations for a non-uniform free beam. We begin by presenting the theory for a uniform beam.<br\>

An example of the motion for a non-uniform beam is demonstrated below:

Vibration of variable beam

Equations for a beam

There are various beam theories that can be used to describe the motion of the beam. The simplest theory is the Bernoulli-Euler Beam theory (other beam theories include the Timoshenko Beam theory and Reddy-Bickford Beam theory where shear deformation of higher order is considered). For a Bernoulli-Euler Beam, the equation of motion is given by the following

[math]\displaystyle{ \partial_x^2\left(\beta(x)\partial_x^2 \zeta\right) + \gamma(x) \partial_t^2 \zeta = p }[/math]

where [math]\displaystyle{ \beta(x) }[/math] is the non dimensionalised flexural rigidity, and [math]\displaystyle{ \gamma }[/math] is non-dimensionalised linear mass density function. Note that this equations simplifies if the plate has constant properties (and that [math]\displaystyle{ h }[/math] is the thickness of the plate, [math]\displaystyle{ p }[/math] is the pressure and [math]\displaystyle{ \zeta }[/math] is the plate vertical displacement) .

The edges of the plate can satisfy a range of boundary conditions. The natural boundary condition (i.e. free-edge boundary conditions).

[math]\displaystyle{ \partial_x^2 \zeta = 0, \,\,\partial_x^3 \zeta = 0 }[/math]

at the edges of the plate.

The problem is subject to the initial conditions

[math]\displaystyle{ \zeta(x,0)=f(x) \,\! }[/math]
[math]\displaystyle{ \partial_t \zeta(x,0)=g(x) }[/math]

Solution for a uniform beam in eigenfunctions with no pressure

If the beam is uniform the equations can be written as

[math]\displaystyle{ \beta \frac{\partial^{4}\zeta}{\partial x^{4}} + \gamma \frac{\partial^{2}\zeta}{\partial t^{2}}=0 }[/math]

We can express the deflection as the series

[math]\displaystyle{ \zeta(x,t)=\sum_{n=0}^{\infty} A_n X_n(x) \cos(k_n t) + \sum_{n=2}^{\infty}B_n X_n(x) \frac{\sin(k_n t)}{k_n} }[/math]

where [math]\displaystyle{ X_n }[/math] are the Eigenfunctions for a Uniform Free Beam and [math]\displaystyle{ k_m = \lambda^2_n \sqrt{\beta/\gamma} }[/math] where [math]\displaystyle{ \lambda_n }[/math] are the eigenfunctions.

Then [math]\displaystyle{ A_n \,\! }[/math] and [math]\displaystyle{ B_n \,\! }[/math] can be found using orthogonality properties:

[math]\displaystyle{ A_n=\frac{\int_{-L}^{L}f(x)X_n(x)\mathrm{d}x}{\int_{-L}^{L}X_n(x)X_n(x)\mathrm{d}x} \,\! }[/math]
[math]\displaystyle{ B_n=\frac{\int_{-L}^{L}g(x)X_n(x)\mathrm{d}x}{\int_{-L}^{L}X_n(x)X_n(x)\mathrm{d}x} }[/math]

Note that we cannot give the plate an initial velocity that contains a rigid body motions which is why the sum starts at [math]\displaystyle{ n=2 }[/math] for time derivative.

Variational Techniques

As an alternative to solving the eigenvalue problem, we can equivalently minimise the Rayleigh Quotient (Linton and McIver 2001):

[math]\displaystyle{ R[\zeta]=\frac{\varepsilon[\zeta]}{H[\zeta]} }[/math]

where [math]\displaystyle{ \varepsilon[\zeta] \,\! }[/math] is known as the energy functional, and [math]\displaystyle{ H[\zeta] \,\! }[/math] is some constraint.

In the case of a uniform beam with zero transverse load, our energy functional is (Lanczos 1949):

[math]\displaystyle{ \varepsilon[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 v}{\partial x^2} \right)^2 \mathrm{d}x }[/math]

If we choose

[math]\displaystyle{ H[\zeta]=\int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1 }[/math]

Then minimizing the Rayleigh Quotient can be expressed as a variational problem

min [math]\displaystyle{ R[\zeta]=\frac{\beta}{2}\int_{-L}^{L}\left(\frac{\partial^2 \zeta}{\partial x^2} \right)^2 \mathrm{d}x }[/math]

subject to [math]\displaystyle{ \int_{-L}^{L}\gamma \zeta^2 \mathrm{d}x=1 }[/math]

which will in turn, also solve our eigenvalue problem.

Rayleigh-Ritz method

We can essentially replace the variational problem of finding a [math]\displaystyle{ y(x) \,\! }[/math] that extremises [math]\displaystyle{ J[y] \,\! }[/math] to finding a set of constants [math]\displaystyle{ a_1, a_2, ..., a_N \,\! }[/math] which extremise [math]\displaystyle{ J[a_1, a_2, ..., a_N] \,\! }[/math], through using the Rayleigh-Ritz method.

We solve

[math]\displaystyle{ \frac{\partial}{\partial a_k}J[a_1, a_2, ..., a_N]=0 \,\! }[/math]

for all [math]\displaystyle{ k=1,2,...,N \,\! }[/math].

In our example of non-uniform beams, the assumption is made that we are able to approximate the eigenfunctions for a non-uniform beam as a linear combination of the eigenfunctions for a linear beam:

[math]\displaystyle{ \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\! }[/math]

Solution for a Non-uniform free beam in eigenfunctions

In the case of a non-uniform free beam, the Euler-Bernoulli beam equation becomes:

[math]\displaystyle{ \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2}\zeta}{\partial x^{2}}\right) +\gamma(x)\frac{\partial^{2}\zeta}{\partial t^{2}}=0 \,\! }[/math]

Using separation of variables, where [math]\displaystyle{ \zeta(x,t)=\widehat{X}(x)\widehat{T}(t) \,\! }[/math], we obtain a corresponding eigenfunction problem (with eigenvalues denoted by [math]\displaystyle{ \mu_n=\widehat{k}_n^4 \,\! }[/math]):

[math]\displaystyle{ \frac{\partial^2}{\partial x^2}\left( \beta(x) \frac{\partial^{2} \widehat{X}_n}{\partial x^{2}}\right)=\gamma(x)\mu_n \widehat{X}_n \,\! }[/math]

which leads us to the following variational expression:

min [math]\displaystyle{ J[\widehat{X}_n]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)(\widehat{X}_n^{''})^2 -\mu_n \gamma(x) \widehat{X}_n^2 \bigg\}\mathrm{d}x }[/math]

We can approximate this using Rayleigh-Ritz and obtain:

[math]\displaystyle{ J[\vec{a}]=\frac{1}{2}\int_{-L}^{L} \bigg\{\beta(x)[\sum_{i=1}^{N}a_{i}{X}_{i}^{''}]^2 -\mu_i \gamma(x) [\sum_{i=1}^{N}a_{i}{X}_{i}]^2 \bigg\}\mathrm{d}x \,\! }[/math]

Then extremising [math]\displaystyle{ J[\vec{a}] \,\! }[/math] as follows

[math]\displaystyle{ \frac{\partial}{\partial a_k}J[\vec{a}]=0 \,\! }[/math]

for all [math]\displaystyle{ k=1,2,...,N \,\! }[/math], allows us to obtain:

[math]\displaystyle{ \frac{\partial}{\partial a_k}J[\vec{a}]=\sum_{i=1}^{N} \left\{ \int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x -\mu_i \int_{-L}^{L} \gamma(x){X}_{i}{X}_{k}\mathrm{d}x\right\}a_i=0 \,\! }[/math]

for all [math]\displaystyle{ k=1,2,...N \,\! }[/math] (here [math]\displaystyle{ \mu_i \,\! }[/math] denotes both the Lagrange multiplier and eigenvalues of a non-uniform beam)

Converting to Matrix Form

If we define

[math]\displaystyle{ K_{ik}=\int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\! }[/math]
[math]\displaystyle{ M_{ik}=\int_{-L}^{L} \gamma(x){X}_{i}{X}_{k}\mathrm{d}x \,\! }[/math]

then for all [math]\displaystyle{ k=1,2,...N \,\! }[/math]:

[math]\displaystyle{ \frac{\partial}{\partial a_k}J[\mathbf{a}]=\sum_{i=1}^{N}(K_{ik}-\mu_i M_{ik})a_i=0 \,\! }[/math]

Let us create the matrix [math]\displaystyle{ K \,\! }[/math], with elements [math]\displaystyle{ K_{ik} \,\! }[/math], the matrix [math]\displaystyle{ M \,\! }[/math], with elements [math]\displaystyle{ M_{ik} \,\! }[/math], the sparse matrix [math]\displaystyle{ \Lambda \,\! }[/math] with [math]\displaystyle{ \mu_i \,\! }[/math] terms on the diagonal, and the vector [math]\displaystyle{ \mathbf{a} \,\! }[/math]. We can consequently express the above sum in the following way:

[math]\displaystyle{ (K-\Lambda M)\mathbf{a} = 0 \,\! }[/math]

So we can solve for [math]\displaystyle{ \mathbf{a} \,\! }[/math] in any of the problems below to obtain coefficients [math]\displaystyle{ a_i \,\! }[/math]:

[math]\displaystyle{ (K-\Lambda M)\mathbf{a} = 0 \,\! }[/math]
[math]\displaystyle{ K\mathbf{a}=\Lambda M\mathbf{a} \,\! }[/math]
[math]\displaystyle{ (M^{-1}K)\mathbf{a}=\Lambda \mathbf{a} \,\! }[/math]

Quadrature

Rather than tediously evaluating an integral for each element of the matrices [math]\displaystyle{ K \,\! }[/math] and [math]\displaystyle{ M \,\! }[/math], we can break up the integral into subintervals with weights [math]\displaystyle{ w_h \,\! }[/math] using Composite Simpson's Rule:

[math]\displaystyle{ \int_a^b f(x) dx \approx \frac{h}{3} \left[ f(x_0)+2 \sum_{j=1}^{n/2-1}f(x_{2j})+4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n)\right] \,\! }[/math]
[math]\displaystyle{ \approx \frac{h}{3}\left[ f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+...+4f(x_{n-1})+f(x_n)\right] \,\! }[/math]

So we can express [math]\displaystyle{ K_{ik} \,\! }[/math] as follows:

[math]\displaystyle{ K_{ik}=\int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\! }[/math]
[math]\displaystyle{ \approx \sum \beta(x_h){X}_{i}^{''}(x_h) {X}_{k}^{''}(x_h)w_h \,\! }[/math]
[math]\displaystyle{ \approx \vec{X}_{i}^{''} H \vec{X}_{k}^{''T} \,\! }[/math]

where

[math]\displaystyle{ \vec{X}_i^{''}=[X_i^{''}(x_1),X_i^{''}(x_2),...,X_i^{''}(x_h)] \,\! }[/math]

[math]\displaystyle{ H = \begin{bmatrix} \beta(x_1)w_1 & 0 & ... &0 \\ 0 & \beta(x_2)w_2 & ... &0 \\ \vdots & & \ddots&\vdots \\ 0 & ... & ...&\beta(x_h)w_h \end{bmatrix} \,\! }[/math]

and weights [math]\displaystyle{ w_h \,\! }[/math] are defined as: [math]\displaystyle{ w_1=h/3, w_2=4h/3, ... , w_h=h/3 \,\! }[/math].

We can extend this concept to the the full matrix [math]\displaystyle{ K \,\! }[/math] if we form:

[math]\displaystyle{ K=X^{''}_{mat}H X^{''T}_{mat} \,\! }[/math]

where

[math]\displaystyle{ X^{''}_{mat} = \begin{bmatrix} {X}_{1}^{''}(x_1) & {X}_{1}^{''}(x_2) & ... & {X}^{''}_{1}(x_h) \\ {X}_{2}^{''}(x_1) & {X}_{2}^{''}(x_2) & ... & {X}^{''}_{2}(x_h) \\ \vdots & \vdots & & \vdots \\ {X}_{N}^{''}(x_1) & {X}_{N}^{''}(x_2) & ... & {X}^{''}_{N}(x_h) \end{bmatrix} \,\! }[/math]

For [math]\displaystyle{ M \,\! }[/math] we use an identical approach:

[math]\displaystyle{ M=X_{mat}J X^{T}_{mat} \,\! }[/math]

where

[math]\displaystyle{ J = \begin{bmatrix} \gamma(x_1)w_1 & 0 & ... &0 \\ 0 & \gamma(x_2)w_2 & ... &0 \\ \vdots & & \ddots&\vdots \\ 0 & ... & ...&\gamma(x_h)w_h \end{bmatrix} \,\! }[/math]

[math]\displaystyle{ X_{mat} = \begin{bmatrix} {X}_{1}(x_1) & {X}_{1}(x_2) & ... & {X}_{1}(x_h) \\ {X}_{2}(x_1) & {X}_{2}(x_2) & ... & {X}_{2}(x_h) \\ \vdots & \vdots & & \vdots \\ {X}_{N}(x_1) & {X}_{N}(x_2) & ... & {X}_{N}(x_h) \end{bmatrix} \,\! }[/math]

Evaluation

Using the MATLAB program eig.m, we can solve the eigenvalue problem

[math]\displaystyle{ K\vec{a}=\Lambda M\vec{a} \,\! }[/math]

using the either of the commands

[Ai_values,nonuniform_eigvals] = eig(M^(-1)*K)
[Ai_values,nonuniform_eigvals] = eig(K,M)

Where the diagonal of nonuniform_eigvals denotes the eigenvalues for the nonuniform beam ([math]\displaystyle{ \mu_n \,\! }[/math]), and each column of the matrix Ai_values represents the coefficients [math]\displaystyle{ a_1, a_2,... \,\! }[/math] corresponding to [math]\displaystyle{ \widehat{X}_n \,\! }[/math]:

[math]\displaystyle{ \widehat{X}_n=\sum_{i=1}^{N}a_{i}X_{i}(x) \,\! }[/math]

We can obtain the matrix [math]\displaystyle{ \widehat{X}_{mat} \,\! }[/math] by taking Xhat_mat= Ai_valuesTX_mat.

Non-uniform beam revisited

We have already solved the eigenfunction problem. We now turn our attention to the time component arising from separation of variables:

[math]\displaystyle{ \widehat{T}_n^{''}+\mu_n\widehat{T}_n=0 \,\! }[/math]

which has solutions of the form

[math]\displaystyle{ \widehat{T}_n(t)=A_n \cos(\widehat{k}_n^2 t)+B_n \sin(\widehat{k}_n^2 t) \,\! }[/math]

Introducing the initial conditions

[math]\displaystyle{ v(x,0)=f(x) \,\! }[/math]
[math]\displaystyle{ \frac{\partial v(x,0)}{\partial t}=g(x)\,\! }[/math]

Then using the the first of these initial conditions we obtain:

[math]\displaystyle{ \sum_{n=0}^{\infty}A_n\widehat{X}_n(x)=f(x) \,\! }[/math]

multiplying by [math]\displaystyle{ m(x)\widehat{X}_n(x) \,\! }[/math] and integrating allows us to obtain:

[math]\displaystyle{ A_n=\frac{\int_{-L}^{L}f(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\! }[/math]

Then using the second initial condition, and applying the same technique yields

[math]\displaystyle{ B_n=\frac{1}{\widehat{k}_n^2}\frac{\int_{-L}^{L}g(x)m(x)\widehat{X}_n(x)\mathrm{d}x}{\int_{-L}^{L}m(x)\widehat{X}_n(x)\widehat{X}_n(x)\mathrm{d}x} \,\! }[/math]

Consequently,

[math]\displaystyle{ v(x,t)=A_0\widehat{X}_0(x)+A_1\widehat{X}_1(x)+\sum_{n=2}^{\infty} A_n\widehat{X}_n(x)\cos(\widehat{k}_n^2 t)+\sum_{n=2}^{\infty} B_n\widehat{X}_n(x) \sin(\widehat{k}_n^2 t) \,\! }[/math]

where [math]\displaystyle{ \widehat{X}_0(x) \,\! }[/math] and [math]\displaystyle{ \widehat{X}_2(x) \,\! }[/math] are the two rigid modes of the non-uniform beam. Note that [math]\displaystyle{ B_0 \,\! }[/math] and [math]\displaystyle{ B_1 \,\! }[/math] do not exist.

Matlab Code

A program to solve for a variable beam can be found here variable_beam.m

Additional code

This program requires

References

Lanczos 1949, Linton and McIver 2001, Rao 1986