Sections24-05.mws

Unit 5: Boundary Value Problems for Partial Differential Equations

Chapter 24: Wave Equation

Sections24.5: longitudinal vibrations in an elastic rod

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Intializations

> restart;

> with(plots):
interface(showassumed=0);

Warning, the name changecoords has been redefined

>

Introduction

A rod of length L and density rho (mass per unit volume) is deformed by pulling or pushing in opposite directions at its ends. The rod is said to stretch of compress linearly along its axis if the restorative forces induced in the rod obey Hooke's law. In effect, the rod behaves like a spring which is both extensible and compressible. Such a rod is called elastic if after a deformation, it returns to its original state.

Under compression, the distance between contiguous cross-Sections in the rod decreases, while under expansion (rerefaction), it increases. Compression or expansion experienced by neighboring plane cross-Sections can propagate through the rod as a wave which obeys the wave equation

u[tt] = c^2*u[xx]

Since the disturbance travels along the axis of the rod, the wave is said to be longitudinal . Thus, a steel rod, hammered on its end, will support longitudinal waves of compression and rarefaction traveling the length of the rod. Similarly, the rubber engine mounts supporting the engine in an automobile, when subjected to shock, experience the same passage of distortions progressing as plane waves through the medium.

Although the longitudinal vibrations in a rod are governed by the same wave equation as the vibrating string, the meaning of the displacement u(x,t) is different, and subtly so. In the rod, u(x,t) measures the displacement, at time t , of the plane in the rod that was at x initially. Since these displacements are relative to where the plane Sectionswas initially, there is an additivity which makes the complete motion much harder to fathom than the vibrations in the finite string where the displacement was directly observable.

In addition, the free endpoint condition whereby an end is free to move is formalized in terms of the strain , a new concept for the reader not familier with the theory of elasticity. Hence, our approach will be to state an appropriate boundary value problem, nurture intuition by illustrating the solution, and finally, develop the theory and derive the governing partial differential equation.

>

Example 24.6

Problem Statement

Formulate and solve an appropriate boundary value problem for the longitudinal elastic vibrations in a rod of length Pi whose left and right ends are both free. Assume c = 1 as the wave speed, and uniformly stretch the rod by one unit, releasing it without imparting any initial velocity.

>

Problem Formulation

If the rod is placed between x = 0 and x = Pi on an x -axis, and if a copy of the x -axis is etched on the rod at equilibrium, then the desired boundary value problem consists of the wave equation

u[tt] = u[xx]

the boundary conditions

u[x](0,t) = 0

u[x](Pi,t) = 0

and the initial conditions

u(x,0) = f(x) = x/Pi

u[t](x,0) = 0

The function u(x,t) measures the displacement, at time t , of the plane Sectionsin the rod that was at location x initially. The endpoint conditions are derivatives because, as we will see later whenwe derive the governing equations, "free ends" mean no forces act on the ends. In fact, the elastic force on a face at location x is given by u[x](x,t) .

To determine the initial distribution of displacements in the rod, we must first obtain the constant value of the strain , defined as the change in length, per unit length. Thus, we compute the constant strain as

Delta*L/L = ([Pi+1]-[Pi])/Pi = 1/Pi

The initial distribution of displacements is then found by accumulating the strain along the rod, obtaining

u(x,0) = Int(1/Pi,s = 0 .. x) = x/Pi

Finally, the initial velocity, namely, u[t](x,0) , is zero.

This boundary value problem is new because we have not yet solved the wave equation with homogeneous Neumann conditions (as derivative conditions at the endpoints are called.)

>

Solution by Separation of Variables

Assume u(x,t) = X(x)*T(t) , that is,

> U := X(x)*T(t);

U := X(x)*T(t)

>

as the separated form for the solution. As in Sections 24.1 and 24.2, the wave equation then yields the two ODEs

Apply this separation assumption to the endpoint conditions, obtaining

`X''`(X)-lambda*X(x) = 0

and

`T''`(t)-lambda*T(t) = 0

However, applying the separation assumption to the endpoint conditions now gives

u[x](0,t) = 0 => `X'`(0)*T(t) = 0 => `X'`(0) = 0

u[x](Pi,t) = 0 => `X'`(Pi)*T(t) = 0 => `X'`(Pi) = 0

In addition, applying the separation assumption to the initial velocity yields

u[t](x,0) = 0 => X(x)*`T'`(0) = 0 => `T'`(0) = 0

a result we obtained each time we solved the wave equation with a similarly vanishing initial velocity.

We therefore face the Sturm-Liouville system

`X''`(x)-lambda*X(x) = 0

`X'`(0) = `X'`(Pi) = 0

and the modified initial value problem

`T''`(t)-lambda*T(t) = 0

`T'`(0) = 0

From Sections16.1, the Sturm-Liouville BVP for X(x) has the eigenvalues

lambda = -n^2 , n = 0, 1, 2, ...

and the corresponding eigenfunctions

X[n] = A[n]*cos(n*x) , n = 0, 1, 2, ...

The corresponding solutions for T(t) are then

T[n] = B[n]*cos(n*t) , n = 0, 1, 2, ...

A typical eigensolution is therefore

u[n](x,t) = a[n]*cos(n*x)*cos(n*t)

and u(x,t) is a linear combination of all such eigensolutions. Evaluating at t = 0 tells us the coefficients in the sum are the Fourier coefficients for f(x) , so we have the solution

u(x,t) = a[0]/2 + Sum(a[n]*cos(n*x)*cos(n*t),n = 1 .. infinity)

a[n] = 2/Pi Int(f(x)*cos(n*x),x = 0 .. Pi)

For f(x) = x/Pi , that is, for

> f := x/Pi;

f := x/Pi

>

the Fourier coefficients are

a[0] = 2/Pi Int(x/Pi,x = 0 .. Pi) = 1

and

a[n] = 2/Pi Int(``(x/Pi)*cos(n*x),x = 0 .. Pi) = 2 ((-1)^n-1)/(Pi^2)/(n^2)

which we obtain in Maple as follows. First, obtain a[0] with the separate integral

> q0 := (2/Pi)*Int(f,x=0..Pi);

q0 := 2*Int(x/Pi,x = 0 .. Pi)/Pi

>

whose value is

> a0 := value(q0);

a0 := 1

>

The general Fourier coefficient is given by

> qa := (2/Pi)*Int(f*cos(n*x),x=0..Pi);

qa := 2*Int(x*cos(n*x)/Pi,x = 0 .. Pi)/Pi

>

Tell Maple that n is an integer via

> assume(n,integer);

>

and evaluate this defining integral, obtaining

> a := simplify(value(qa));

a := 2*((-1)^n-1)/(Pi^2*n^2)

>

The partial sum

u[20](x,t) = 1/2 + 2 Sum(``(((-1)^n-1)/(Pi^2)/(n^2))*cos(n*x)*cos(n*t),n...

that is,

> u20 := a0/2 + sum(a*cos(n*x)*cos(n*t),n=1..20):

>

is graphed as a surface in Figure 24.18, below.

> plot3d(u20, x=0..Pi, t=0..4*Pi, axes=frame, color=black, labels=[x,t,u], labelfont=[TIMES,ITALIC,12], scaling=constrained, tickmarks=[3,6,2], orientation=[-45,50]);

[Maple Plot]

>

The wave-like properties of the solution are apparent, but the detailed interpretation of the meaning of these waves requires more investigation. Remember, u(x,t) gives the displacement, at time t , of the plane Sectionsthat was at x initially. What we want to see is a representation of the physical motion of the elastic rod, and unlike the transverse waves on the finite string, u(x,t) does not directly show these physical motions.

An animation of the solution u(x,t) is a dynamic interpretation of the solution surface, but is still not a representation of the physical motions in the rod. It shows the displacements u(x,t) which each face that was at location x is now experiencing at time t . It tells us that whatever is happening is happening linearly, and it takes place from the endpoints in, and back, with a wave-front like a shock wave. Several frames from this animation appear in Figure 24.19, but the full animation is generated in Maple by executing the following command.

> p1 := animate(u20, x=0...Pi, t=0..4*Pi, frames=80, labels=[x,u], color=black, labelfont=[TIMES,ITALIC,12], xtickmarks=3, ytickmarks=[0,1], scaling=constrained):
p1;

[Maple Plot]

>

What does this animation mean? What is it telling us about the longitudinal vibrations in the rod? For openers, the animation shows the displacements u(x,t) which each face that was at location x is now experiencing at time t . Unfortunately, it is still hard from this animation to comprehend the actual motions taking place in the rod.

However, it does tell us that whatever is happening is happening linearly, and it takes place from the endpoints in, and back, with a wave-front like a shock wave.

A more revealing animation simulates the actual rod by imagining a set of eleven equispaced scratch marks on the rod, and following the motion of these scratch marks in real time. Figure 24.20 shows the bar at t = 3.7 but the full animation is generated by the following Maple commands.

> uu := unapply(u20, x, t):
F1 := (X,T) -> [[X+uu(X,T),0],[X+uu(X,T),1]]:
F2 := T -> {seq(F1(k*Pi/10,T), k = 0..10)}:
F3 := T -> plot(F2(T), color=black):
p2 := display([seq(F3(k*Pi/20), k = 0..79)], insequence=true, scaling=constrained, ytickmarks=2, labels=[x,``]):
p2;

[Maple Plot]

>

The uniformly stretched rod relaxes to its natural length, with the collapse occurring from the "outside, in." The center of the rod is the only part of the rod that does not move. Then, the collapse continues past the point of reaching natural length, and the rod contracts by an amount equal to the original stretch. This continuing contraction now takes place from "inside, out" as the center of the rod first experiences the compression, with the outside of the rod feeling the compression last. The rod then expands from this state of full compression to natural length, from "outside, in." Having achieved natural length, the rod continues expanding to its initial stretched state, this time from the "inside, out."

In addition to this simulation of the rod itself, it is also instructive to examine an animation of u[x](x,t) , which is the localized (per-face) value of the strain, Delta*L/L . This animation shows more clearly how the compression and rarefaction waves travel with a "front" through the rod.

To compute the derivative u[x](x,t) we have to differentiate the infinite sum representing the solution u(x,t) . The derivative of a Fourier series converges even more slowly than the original series. Where we obtained acceptable results with a twenty-term approximation for u(x,t) itself, it will take 100 terms for the derived series to converge. This translates into a modest increase in computation time.

Figure 24.21 shows, at times t = 1 and t = 4 , the strain in the rod. However, executing the following Maple commands will generate the complete animation.

> u100 := a0/2 + sum(a*cos(n*x)*cos(n*t),n=1..100):
ux100 := diff(u100,x):
p3 := animate(ux100, x=0..Pi, t=0..4*Pi, frames=80, color=black, labels=[x,strain], xtickmarks=3, ytickmarks=[-.3,0,.3]):
p3;

[Maple Plot]

>

The animation shows, as a function of time, the changing strain in the rod. The strain moves longitudinally through the rod, as a shock wave. It starts out uniform across the rod, and drops to zero with a shock-front. When the strain is completely gone, the rod has reached its natural length. As the rod then enters a compressive mode, the strain becomes negative, and the region of compression expands from the center out to the ends. At this juncture the strain graph is a constant negative amount.

The compression then starts to vanish from the ends inwards until the rod again reaches its natural length. This corresponds to the vanishing of the strain. As the rod enters an expansion phase, the strain appears above the x -axis, growing, with a shock-front, outwards from the center. The expansion continues until the rod returns to its original stretched state.

Finally, we place side-by-side the three animations of displacement, physical motion, and strain. Run continuously, and slowly enough, the resulting animation allows the three facets of the solution to be compared and coordinated.

> display(array([p||(1..3)]));

[Maple Plot]

>

Theory and Derivation

A rod which obeys a linear law of elasticity is said to obey Hooke's law , wherein forces are proportional to displacement. This linear relationship between force and displacement, expressed by the formula F = k*x , is a staple of elementary calculus and general physics, especially in any discussion of the behavior of a spring. In fact, the linear spring was an essential element in Chapters 5 and 6. Thus, F , the force required to stretch a spring, is proportional to x , the amount of the stretch. The constant of proportionality, k , is called the spring constant.

In the context of elasticity, Hooke's law is sometimes stated as

stress is proportional to the strain

where stress is the force per unit cross-sectional area and

strain = Delta*L/L

is defined as the elongation per unit length. In the context of elasticity, Hooke's law is then

stress = E strain

where E is Young's modulus, the constant of proportionality between stress and strain . This gives Young's modulus the units of

E ~ stress/strain ~ ``(force/area)/``(length/length) = force/area

In fact, starting with F = k*x , and dividing through by A , the cross-sectional area, we have

F/A = k/A x

The displacement x is the elongation Delta L from the definition of strain, namely, strain = Delta*L/L . Hence, the spring law becomes

stress = F/A = k/A Delta L

= k/A L*strain

= E*strain

where E , Young's modulus, is k*L/A . Hence, the two statements of Hooke's law are equivalent.

Since longitudinal waves travel in the direction of the rod's axis, we want a coordinate system along that axis. However, we must assume a linear elasticity which keeps plane Sections both plane and parallel. Thus, forces within the rod obey the linear Hooke's law and preserve the orientation of all plane Sections which are perpendicular to the axis of the rod.

Put the rod (whose equilibrium length is L ) on a ruler, aligning the left end of the rod with the left end of the ruler. At equilibrium, paint a copy of the ruler's coordinates on the rod. For clarity, take this linear coodinate system to be an x -coordinate system.

Identify a face (i.e., a plane section) within the rod by the coordinate it had at equilibrium. If the equilibrium coordinates are painted on the rod, then these coordinates travel with the rod should the rod be stretched or compressed. In addition, leave the ruler on the table so the original meaning of an x -coordinate is preserved.

Consider R , a segment of the rod bounded at equilibrium by faces at x and x+dx . This segment has length dx and is shown in the following figure, FIgure 24.22.

> p4 := plot([10+2*cos(t),3+3*sin(t),t=0..2*Pi],color=black):
p5 := plot([2*cos(t),3+3*sin(t),t=Pi/2..3*Pi/2],color=black):
p6 := plot([3+2*cos(t),3+3*sin(t),t=Pi/2..3*Pi/2],color=black):
p7 := plot([5+2*cos(t),3+3*sin(t),t=Pi/2..3*Pi/2],color=black):
p8 := plot([5+2*cos(t),3+3*sin(t),t=-Pi/2..Pi],color=black, linestyle=2):
p9 := plot(6,x=0..10,color=black):
p10 := plot(0,x=0..12,color=black):
p11 := textplot({[2,-.5,`x`], [5,-.5,`x+dx`], [12.5,0,`x`], [10,-.5,`L`]}, font=[TIMES,ITALIC,12]):
p12 := textplot([2.4,3,`R`], font=[TIMES,ITALIC,12]):
display([p||(4..12)],scaling=constrained,view=[-2..13,-1..6],axes=none);

[Maple Plot]

>

Let u(x,t) denote the displacement experienced, at time t , by the face, which at equilibrium, was located at x . It helps to remember that this face still bears the coordinate x . The displacement undergone by this face is the quantity u(x,t) , and is measured along the axis of the rod. Unlike the vibrating string where the displacements are very obvious, the lateral displacements of faces in the rod are much more difficult to visualize.

The partial differential equation governing longitudinal vibrations in the rod will be a consequence of F = m a , Newton's second law of motion. It's not hard to guess that a , the acceleration of the face with the coordinate x painted on it, is given by a = u[tt](x,t) . Determining F , the net force acting on this face, is more subtle.

At time t , the left and right faces of the segment R are now located respectively at

> X[left] := x+u(x,t);
X[right] := (x+dx)+u(x+dx,t);

X[left] := x+u(x,t)

X[right] := x+dx+u(x+dx,t)

>

The length of the deformed segment is then

> Delta := X[right]-X[left];

Delta := dx+u(x+dx,t)-u(x,t)

>

The change in length of this segment is therefore the relative displacement between the left and right faces of the segment, namely

> Delta - dx;

u(x+dx,t)-u(x,t)

>

the relative displacement between the left and right faces of the segment. The strain associated with the face which, at equilibrium, was at location x , is given by the limiting ratio

Limit((Delta-dx)/dx,dx = 0) = Limit((u(x+dx,t)-u(x,...

Notice that strain, an extensive property, is ascribed to a point by this limit. This definition parallels the definition of pressure which is also an extensive property, force per unit area, ascribed to a point.

Thus, the strain on a single face with coordinate x painted on it is

> q3 := Limit((u(x+dx,t)-u(x,t))/dx,dx=0);

q3 := Limit((u(x+dx,t)-u(x,t))/dx,dx = 0)

>

so that the value of the strain is strain = u[x](x,t) , that is,

> strain = convert(value(q3),diff);

strain = diff(u(x,t),x)

>

This is a fundamental result in the theory of elasticity. Getting this right at the beginning is well worth the time invested in getting it. Thus, at time t and location x , the strain caused by the displacement u(x,t) , is positive if the applied force tends to increase the size of the region R, and negative otherwise.

The force on a plane Sectionscan be obtained from the strain, using Hooke's law, stress = E strain . With A as the cross-sectional area of the rod, we have

force/area = stress = E strain

so that

force = E A strain = E A u[x](x,t)

Thus, the elastic force on the face whose equilibrium coordinate was x , is given by

F(x,t) = E*A*u[x](x,t)

To represent this force in Maple, define

> F := (x,t) -> E*A*D[1](u)(x,t);

F := proc (x, t) options operator, arrow; E*A*D[1](...

>

so that

> 'F(x,t)' = convert(F(x,t),diff);

F(x,t) = E*A*diff(u(x,t),x)

>

The net force on segment R is given by the sum of the forces acting on the left and right faces of R , namely,

E*A*u[x](x+dx,t)-E*A*u[x](x,t)

which, is Maple, is

> `net force` := F(x+dx,t)-F(x,t);

`net force` := E*A*D[1](u)(x+dx,t)-E*A*D[1](u)(x,t)...

>

The total mass of segment R is given by

m = [rho*A*dx]

since the volume of the segment is A , the cross-sectional area, times the length, dx .

Taking the acceleration of a face as u[tt](x,t) , we can write the i -component of Newton's second law, F = m a , as

E*A*u[x](x+dx,t)-E*A*u[x](x,t) = rho*A*dx*u[tt](x,t...

that is, as

> q4 := `net force` = rho*A*dx*diff(u(x,t),t,t);

q4 := E*A*D[1](u)(x+dx,t)-E*A*D[1](u)(x,t) = rho*A*...

>

Divide by A*rho*dx to get

E/rho (u[x](x+dx)-u[x](x,t))/dx = u[tt](x,t)

that is,

> q5 := simplify(q4/(A*rho*dx));

q5 := -E*(-D[1](u)(x+dx,t)+D[1](u)(x,t))/(rho*dx) =...

>

which, in the limit as dx*``-``*`>`*0 , becomes

E/rho u[xx](x,t) = u[tt](x,t)

that is,

> q6 := convert(map(limit,q5,dx=0),diff);

q6 := E*diff(u(x,t),`$`(x,2))/rho = diff(u(x,t),`$`...

>

Upon defining c^2 = E/rho , the partial differential equation becomes u[tt] = c^2*u[xx] , that is

> algsubs(E/rho=c^2,q6);

diff(u(x,t),`$`(x,2))*c^2 = diff(u(x,t),`$`(t,2))

>

With a slight rearrangement and the introduction of subscript notation for derivatives, this equation can be written as

u[tt] = c^2*u[xx]

the very same wave equation we derived for the vibrating string. The difference here is that u(x,t) represents the displacement of the face that was at x at time t , and this displacement is lateral, along the axis of the rod.

>