Waves on a Variable Beam

From WikiWaves
Jump to: navigation, search


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] \partial_x^2\left(\beta(x)\partial_x^2 \zeta\right) + \gamma(x) \partial_t^2 \zeta = p [/math]

where [math]\beta(x)[/math] is the non dimensionalised flexural rigidity, and [math]\gamma [/math] is non-dimensionalised linear mass density function. Note that this equations simplifies if the plate has constant properties (and that [math]h[/math] is the thickness of the plate, [math] p[/math] is the pressure and [math]\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] \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] \zeta(x,0)=f(x) \,\! [/math]
[math] \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] \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] \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]X_n[/math] are the Eigenfunctions for a Uniform Free Beam and [math]k_m = \lambda^2_n \sqrt{\beta/\gamma}[/math] where [math]\lambda_n[/math] are the eigenfunctions.

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

[math] 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] 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]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] R[\zeta]=\frac{\varepsilon[\zeta]}{H[\zeta]} [/math]

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

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

[math] \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] 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] 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}\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] 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.

We solve

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

for all [math] 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] \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] \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] \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]):

[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]

which leads us to the following variational expression:

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]

We can approximate this using Rayleigh-Ritz and obtain:

[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]

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

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

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

[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]

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)

Converting to Matrix Form

If we define

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

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

[math] \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] 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:

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

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

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


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 Composite Simpson's Rule:

[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] \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] K_{ik} \,\![/math] as follows:

[math] K_{ik}=\int_{-L}^{L}\beta(x) {X}_{i}^{''} {X}_{k}^{''}\mathrm{d}x \,\![/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] \vec{X}_i^{''}=[X_i^{''}(x_1),X_i^{''}(x_2),...,X_i^{''}(x_h)] \,\![/math]

[math] 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] w_h \,\![/math] are defined as: [math] w_1=h/3, w_2=4h/3, ... , w_h=h/3 \,\![/math].

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

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


[math] 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] M \,\![/math] we use an identical approach:

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


[math] 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] 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]


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

[math] 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] \mu_n \,\![/math]), and each column of the matrix Ai_values represents the coefficients [math] a_1, a_2,... \,\![/math] corresponding to [math] \widehat{X}_n \,\![/math]:

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

We can obtain the matrix [math] \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] \widehat{T}_n^{''}+\mu_n\widehat{T}_n=0 \,\![/math]

which has solutions of the form

[math] \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] v(x,0)=f(x) \,\![/math]
[math] \frac{\partial v(x,0)}{\partial t}=g(x)\,\![/math]

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

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

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

[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]

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

[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]


[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]

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.m

Additional code

This program requires


Lanczos 1949, Linton and McIver 2001, Rao 1986