• Abstract
• Introduction
• Animation
  • Maple code for a mass-spring system
  • Maple code for a rectangular pulse
• Free Harmonic Oscillations
  • Undamped Motion
  • Overdamped Motion
  • Critically Damped Motion
  • Underdamped Motion
• Forced Harmonic Oscillations
  • The response of a mass-spring system to an impulse function
  • The displacement resonance in forced oscillations
• Summary

Visualizing Free and Forced Harmonic Oscillations using Maple V R5

Harald Pleym

Telemark College, Department of Technology

Kjolnes Ring 56, 3914 Porsgrunn, Norway

harald.pleym@hit.no

Abstract

This paper focuses on how the symbolic, numerical and graphical power of the computer algebra system Maple V Release 5 can be used to explore and visualize with animation free harmonic motion with and without damping and vibration of a mass-spring system forced by an external sharp blow. Dirac's delta function is used to model the external force that acts for a very short period of time. It is also demonstrated how displacement resonance in a forced oscillations readily can be obtained from the transfer function of the system represented by a second order differential equation. And the information contained in the system frequency response may also be conveniently displayed in animated graphical form.

Introduction

At Telemark College, Department of Technology in Porsgrunn, Norway we have incorporated the use of the computer algebra system Maple V in the engineering mathematics curriculum from August 1995. With the powerful software program including graphical, symbolic and numerical techniques, Maple V Release 5 is the ideal program to make mathematics more relevant and motivating for engineering students and an excellent demonstration tool for the classroom. The use of Maple has changed the way I teach and hopefully the way the students learn. I now have the opportunity to emphasize the learning of concepts, the visualizatione of concepts and the solution of more realistic engineering applications.

Many of the phenomena we observe in everyday life have a periodic motion. And in many cases the periodic motion is maintained by a periodic driving force. For these forced oscillations the amplitude depends on the frequency of the driving force. There is usually some frequency of the driving force at which the oscillations have their maximum amplitude. This occur when the natural frequency of the vibrating system and the frequency of the driving force are approximately equal. This phenomenon of resonance plays an important role in almost every branch of physics. And the avoidance of destructive resonance vibrations is an important factor in the design of mechanical systems of all types.

Application of differential equations plays an important role in science and engineering. And the most important step in determining the natural frequency of vibration of a system is the formulation of its differential equation. The purpose of this article is to demonstrate how Maple can be used to investigate and visualize with animations:

- simple harmonic motion of a mass-spring system
- the response of a mechanical system (mass-spring system) to an impulse function
- the displacement resonance in forced oscillations

Animation

Coding in Maple does not require expert programmimg skills, so writing a Maple program deals with putting a proc ( ) and an end around a sequence of commands. And it is pretty easy to put prewitten routines together from Maples powerful bulding blocks and plotting facilities. The plottools package provides many useful commands for producing plotting objects, which can be scaled, translated and rotated. The following Maple codes will be used in the coming sections. The first prosedure animate a spring-mass system with and without a dashpot, and the second animate the response to a rectangular pulse by a linear differential second order equation.


Because it is impossible to show animation on paper, the animation figures shows 4 consecutive frames. The animations in this paper can be seen live on the Journal's Web site.

Maple code for a mass-spring system

m: mass, r: damping constant, k: spring constant, x1,x2: positions of the spring
x0 = x(0), xp0=D(x)(0): initial conditions, f: driving force, n: number of frames

> mass_spring:=proc(m,r,k,x1,x2,x0,xp0,tk,f,n )
local mass,deq,init,sol,xk,xu,plt1,plt2,plt,pltxk,rect,base,spring,
spring1,rod,cylinder, piston,fluid,dashpot;
with(plots):with(plottools):
spring:=proc(x1,x2,n) #procedure for the spring
local p1,p2,p3,p4,pn_1,pn,p;
p1:=[x1,0];
p2:=[x1+0.25,0];
p3:=[x2-0.25,0];
p4:=[x2,0];
p:=i->[x1+0.25+(x2-x1-0.5)/n*i,(-1)^(i+1)*0.5];
plot([p1,p2,seq(p(i),i=1..n-1),p3,p4],thickness=2,color=aquamarine);
end:
spring1:=x2->spring(x1,x2,12);
rect:=t->translate(rectangle([-0.2,-0.2],[0.2,0.2],color=red),t,xk(t)+x2+1);
#rect: display the position of the mass on the displacement curve x(t)
mass:=x2->rectangle([x2,-1],[x2+2,1],color=red):
if x0=0 and xk(0.1)=0 then
xu:=1.5
elif x0=0 and xk(0.1)<>0 then
xu:=5.0;
else
xu:=x0;
fi;
rod:=x2->plot([[x2+2,0],[x2+2.5+2*xu,0]],color=grey,thickness=2);
cylinder:=plot([[x2+2.5+xu,-1],[x2+2.5+xu,1],[x2+3.5+3*xu,1],[x2+3.5+3*xu,-1],
[x2+2.5+xu,-1]],color=tan,thickness=2);
piston:=x2->rectangle([x2+2.5+2*xu,-1],[x2+3+2*xu,1],color=grey);
fluid:=rectangle([x2+2.5+xu,-1],[x2+3.5+3*xu,1],color=blue);
dashpot:=x2->display(rod(x2),cylinder,piston(x2),fluid);
deq:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(t);
print(deq);
init:=x(0)=x0, D(x)(0)=xp0:
sol:=dsolve({deq,init},x(t));
print(combine(simplify(sol)));
xk:=unapply(rhs(sol),t);
#xk=x(t) : the displacement of mass
pltxk:=plot([xk(t)+x2+1,x2+1],t=0..max(tk,x2+4+3*xu)+0.5,
color=[blue,grey],numpoints=400);
base:=plot(-1.1,t=0..max(tk,x2+4+3*xu)+0.5,color=brown,thickness=4);
if r=0 then #undamped case
plt1:=x->translate(display(spring1(x2+x),mass(x2+x),base),0,-3);
else #with damping
plt1:=x->translate(display(spring1(x2+x),mass(x2+x),dashpot(x2+x),base),0,-3);
fi;
plt2:=i->display(pltxk,rect(tk/n*i));
plt:=i->display(plt1(xk(tk/n*i)),plt2(i));
display(seq(plt(i),i=0..n),insequence=true,args[11..nargs]);
end:

Maple code for a rectangular pulse

m: mass, r: damping constant, k: spring constant,
a: constant in the impulse function ,
eps: width of a rectangular pulse , n: number of frames

> impulse_func:=proc(m,r,k,a,eps,n)
local h,dlign1,dlign2,sol1,sol2,plt1,plt2,pltd,plt,txt,txtd,kloss,e;
with(plots):
alias(u=Heaviside):
h:=eps->(u(t-a)-u(t-a-eps))/eps;
dlign1:=eps->m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*h(eps);
dlign2:=m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*Dirac(t-a);
print(dlign1(eps), epsilon=eps);
print(dlign2);
sol1:=eps->simplify(dsolve({dlign1(eps),y(0)=0,D(y)(0)=0},y(t)));
sol2:=dsolve({dlign2,y(0)=0,D(y)(0)=0},y(t));
e:= textplot([6,1.3,convert([101],bytes)], font=[SYMBOL,12]):
txt:=eps->textplot([6.9,1.3,cat(` = `,convert(evalf(eps,2),string))]);
txtd:=textplot([4,1.4,`Response to Dirac's delta function`],align={RIGHT,ABOVE});
pltd:=eps->plot(h(eps),t=0..5*Pi,color=green,thickness=2):
plt1:=eps->plot(rhs(sol1(eps)),t=0..5*Pi,color=red,thickness=2):
plt2:=plot(rhs(sol2),t=0..5*Pi,color=blue,thickness=2):
plt:=eps->display(pltd(eps),txt(eps),e,txtd,plt1(eps),plt2);
display(seq(plt(eps-(eps-0.1)/n*i),i=0..n),insequence=true,args[7..nargs]);
end:

Free Harmonic Oscillations

Machines with rotating components commonly involve mass-spring systems or their equivalents in which the driving force is simple harmonic. The motion of a mass attached to a spring serves as a simple example of vibrations that occur in more complex mechanical systems.

From a teaching point of view it is suitable to consider a body of mass m attached to one end of a spring that resists compression as well as stretching. A rod attached to the mass carries a disk moving in an oil-filled cylinder (a dashpot). The other end of the spring could be attached to a fixed wall and vibrating horizontally. The resultant force on the body is the sum of the restoring force - and the damping force - , is the force constant, is the distance of the body of mass from its equilibrium position, is time and and is the damping constant. We take > 0 when the spring is stretched. The differential equation of motion is therefore:

> deq1:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=0;

If we set in the motion is undamped. Otherwise the solution of presents three distinct cases of damping according to whether is greater than, equal to, or less than zero.

> eq:=d=r^2-4*m*k;

> solm:=solve(eq,m);

Let us replace with a new variable in .

> mv:=unapply(solm,d,r,k);

> deq2:=collect(4*k*subs(m=mv(d,r,k),deq1),diff(x(t),t$2));

With this differential equation, , the system is overdamped if , critically damped if and underdamped if . If we request Maple to solve for the undamped case and for each of the damped cases,subject to initial condition and = 0

> init:=x(0)=x[0],D(x)(0)=0: #initial conditions

we get:

With

With

> deq1a:=subs(r=0,deq1);

> dsolve({deq1a,init},x(t));

A typical graph of with and is shown i Figure 1. Figure 2 shows the animation of the undamped motion.

> xk:=unapply(rhs(%),m,k,x[0],t): # defines x = x(m,k,x[0],t)

> interface(showassumed=0):

> assume(d>0):

> dsol:=dsolve({deq2,init},x(t));

> subs(eq,dsol);

> xo:=unapply(rhs(%),m,r,k,x[0],t): # defines x = x(m,r,k,x[0],t)

The solution consists of two exponential terms when > . In all subsequent figures we take and . Figure 3 shows some typical graphs of the position function for the overdamped case. Figure 4 shows the animation of the motion with the damping constant .

> lhs(dsol)=limit(rhs(dsol),d=0);

> subs(r=sqrt(4*m*k),%);

> xc:=unapply(rhs(%),m,k,x[0],t): #defines x = x(m,k,x[0],t)

The solution consists of one exponential term when . The graph in Figure 5 resemble those of the overdamped case in Figure 3. The animation of a critically damped motion with is shown in Figure 6.

> assume(d<0):

> dsolve({deq2,init},x(t));

> subs(eq,%);

> xu:=unapply(rhs(%),m,r,k,x[0],t): #defines x = x(m,r,k,x[0],t)

The solution represents exponentially damped oscillations of the mass-spring system about its equilibrium position as shown in both Figure 7 and Figure 8.

Undamped Motion

> with(plots): #load the plots package

Warning, the name changecoords has been redefined

> plot(xk(1,1,1,t),t=0..19,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Figure 1 Undamped motion. , ,

> mass_spring(1,0,1,0,5,2,0,14,0,60,scaling=constrained);

Figure 2 Animation of an undamped motion. , ,

The small rectangle shows the position of the mass on the displacement curve.

Overdamped Motion

If we select , and three different -values: , named becomes:

> rk:=(d,m,k)->simplify(sqrt(d+4*m*k)):

> r=[seq(rk(d,1,1),d=[12,45,96])];

> plt1:=d->plot(xo(1,rk(d,1,1),1,2,t),t=0..20):

> txt1:=d->textplot([4,xo(1,rk(d,1,1),1,2,4),convert(r = rk(d,1,1),string)]):

> pltod:=display(seq({plt1(d),txt1(d)},d=[12,45,96]),labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]):%;

Figure 3 Overdamped motion, , , , : damping constant

Figure 3 shows three typical graphs of the position function for the overdamped case. Figure 4 shows the animation of the motion with .

> mass_spring(1,4,1,0,5,2,0,14,0,60,scaling=constrained);

Figure 4 Animation of an overdamped motion of a mass-spring system with dashpot.
, , ,

Critically Damped Motion

> plt2:=plot(xc(1,1,2,t),t=0..20):

> txt2:=textplot([3,xc(1,1,2,3),convert(r = rk(0,1,1),string)]):

> pltcd:=display(plt2,txt2):%;

Figure 5 Critically damped motion, , , , : damping constant

> mass_spring(1,2,1,0,5,2,0,14,0,60,scaling=constrained);

Figure 6 Animation of an critically damped motion. , , ,

In this critically damped case, then resistance of the dashpot is just large enough to damp out any oscillations.

Underdamped Motion

We select , and two values .

> plt3:=plot([xu(1,1/2,1,2,t),xu(1,1/5,1,2,t)],t=0..20):

> txt3:=textplot([ [6.5,xu(1,1/2,1,2,6.5),`r=1/2`],[6.5,xu(1,1/5,1,2,6.5),`r=1/5`]]):

> pltod:=display(plt3,txt3,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]):%;

Figure 7 Underdamped motion,  , , , : damping constant

> mass_spring(1,1/5,1,0,5,2,0,20,0,60,scaling=constrained);

Figure 8 Animation of an underdamped motion of a mass-spring system with dashpot. , ,

The action of the dashpot exponentially damps the oscillations in accord with with the time-varying amplitude. The dashpot also decreases the frequency of the motion from in the undamped to case in the underdamped motion with the same mass and force constant.


Comparison of the graphs in Figure 2, 5 and 7 (Figure 9) shows that when the motion is critically damped the mass reaches its equilibrium posistion in a shorter time than when the damping is larger. And for any damping constant less than in our particular example the motion becomes oscillatory.

> display(plt1(45),plt2,plt3,txt1(45),txt2,txt3,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Figure 9 Overdamped, critically damped and damped oscillatary motion of a mass-spring system.
, , , : damping constant

Forced Harmonic Oscillations

Mechanical systems are often acted upon by an external force of large magnitude that acts for only a short period of time. For example a spring-mass system could be given a sharp blow at some specific time . In engineering practice it is convenient to use the Dirac delta function as a mathematical model for such a blow. But it is my experience that it is often difficult for engineering students to get a real understanding of what this impulse function stands for. With the use of Maple it is easy and instructive to simulate Dirac's delta function

The response of a mass-spring system to an impulse function

As said above the Dirac delta function could serve as a mathematical model for an external force of large magntiude that acts for only a very short period of time.

> delta(t-a)=piecewise(t<a,0,t=a,infinity,t>a,0);

where

> Int(delta(t-a),t=-infinity..infinity)=int(Dirac(t-a),t=-infinity..infinity);

Obviously no real function can satisfy beeing zero except at a single point and have an integral equal to one.

In Maple = Dirac (t - a) is expressed as the derivative of the Heaviside standard unit step function

> alias(u=Heaviside):

> Diff(u(t-a),t)=diff(u(t-a),t);

It is instructive to use Maple to model such an instantaneous unit impulse by starting with the function

> delta[epsilon](t-a)=piecewise(t<a,0,t<a+epsilon,1/epsilon, 0);

> h:=unapply(convert(rhs(%),Heaviside),t,a,epsilon);

A plot of the rectangular pulse for and gives:

> plot(h(t,Pi,3),t=0..10,thickness=3);

Figure 10 Rectangular pulse with and (width and height )

The area of the rectangular pulse is equal to 1 .

> assume(epsilon>0):

> A:=Int(h(t,a,epsilon),t=a..a+epsilon):

> A=value(A);

where is Heavisides unit step function. We can compare the response of a damped mass-spring system to a rectangular pulse as with the response of the Dirac's delta function used by Maple.

The output of the Maple code impulse_func (written in section 3) in Figure 11 animate the response to the rectangular pulse by a linear second order differential equation with constant coefficients and with damping.

> impulse_func(1,3,1,Pi,3,10);

Figure 11 The behavior of the response to a rectangular pulse by a linear second order

differential equation with damping as compared to the response to Dirac's delta function.

The last animated frame in Figure 11 shows little difference between the two responses when

Figure 12 shows an animated response of the damped mass-spring system initially at rest. At time the system is suddenly given a sharp "hammerblow" modelled by

> f:=t->5*Dirac(t-Pi):

> mass_spring(1,3,1,0,5,0,0,14,f,60,scaling=constrained);

Figure 12 The motion of a mass-spring system with dashpot

under the influence of a sharp blow at provided by

The displacement resonance in forced oscillations

In every oscillating system there is dissipation of mechanical energy, which results that the motion of the mass-spring system described in the previous sections dies out. If the oscillation are to be maintained, energy must be supplied to the system. In this section we shall assume that the system is acted on by a periodic driving force. Suppose that the mass-spring system is subjected to a periodic force , where is the maximum value of the applied force and is its frequency. The equation of motion is then:

where is the mass of the system, is the distance of the body of mass from its equilibrium position, r is the damping constant and is the force constant. The frequency response of this system can be readily obtained from the system transfer function by replacing by

> f:='f':

> with(inttrans): L:=x->laplace(x,t,s): invL:=X->invlaplace(X,s,t):

> alias(X(s)=L(x(t)),F(s)=L(f(t))):

> dlign:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(t);

> simplify(L(dlign));

> X(s)=solve(%,X(s));

> subs(x(0)=0,D(x)(0)=0,%);

The transfer function of the mass-spring system is given by:

> H:=s->1/(m*s^2+r*s+k):

It is easy to show that the steady-state frequency response to an input becomes:

> x(t)=F[0]*abs('H'(I*omega))*cos(omega*t+arg('H'(I*omega)));

> value(%);

The steady-state system response is also a cosine having the same frequency as the input. And the amplitude of this response is . The variation in both the magnitude and argument as the frequency of the input cosine is varied constitute the frequency response of the system, as the following example shows.

Lets us first solve the differential equation using Maples dsolve .

> f:=(omega,t)->2*cos(omega*t);

> deq:=(m,r,k,omega)->m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(omega,t);

> sol:=dsolve({deq(1,0.25,4,2),x(0)=0,D(x)(0)=0},x(t),method=laplace);

> plot(rhs(sol),t=0..40);

Figure 13 The displacement of a mass-spring system
undergoing forced oscillations plotted against the time.
, , ,

The graph shows that the transient solution dies out as increases.

The transfer function of the system is given by :

> H:=(m,r,k,s)->1/(m*s^2+r*s+k);

With , and , we get:

> 'H'(1,r,4,I*omega)=H(1,r,4,I*omega);

> Habs:=unapply(evalc(abs(H(1,r,4,I*omega))),r,omega);

There is some frequency called the resonance frequency at which the amplitude becomes a maximum.

This resonance frequency can be recognized in many vibrating systems unless the damping force is too large.

> solve(omega^4-8*omega^2+16+r^2*omega^2=0,{r});

With zero damping force, and the resonance frequency in our example. See also Figure 15.

The phase angle is given by:

> phi:=unapply(evalc(argument(H(1,r,4,I*omega))),r,omega);

Substituting the values , and , gives the steady-state response,

> xs:=unapply(2*Habs(0.25,2)*cos(2*t+phi(0.25,2)),t):
'xs'(t)=xs(t); #steady state respons

> sol; # Maple's solution

As , .

> control:='4*sin(2*t)-4*cos(2*t-Pi/2)'=simplify(4*sin(2*t)-4*cos(2*t-Pi/2));

The system is lagging because the phase shift between input and output is = -1.57 ... .

Figure 14 shows this phase shift.

> with(plots):
plt1:=plot([f(2,t),xs(t)],t=30..39,thickness=[1,2]):
plttext1:=textplot([[39,xs(39),`steady-state response`],[39,f(2,39),`input`]],align=RIGHT):
display({plt1,plttext1},labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Figure 14 Phase shift between steady-state response and input displacement

> A:='A':

> Habsi:=i->plot(Habs(i,omega),omega=1..3,A=0..2.5) :

> plt2:=i->display(Habsi(i),text2(i)):

> text2:=i->textplot([2,Habs(i,2),convert(r=i,string)],align=ABOVE,font=[TIMES,BOLD,12]):

> display(seq(plt2(i),i=[0.1,0.25,0.5,0.75,1,1.41]),labels=[`omega`,`A`],labelfont=[TIMES,BOLD,12]);

Figure 15 Variation of the amplitude
against the angular velocity of the forced oscillations-displacement resonance for different values of the damping constant

Figure 15 shows that the smaller the value of , the sharper the resonance of the vibrating system to the applied force.

And the resonance frequency varies with the value of the damping constant . The amplitude becomes infinity when the applied frequency becomes equal to the resonance frequency and the system has zero damping.

Figure 16 animates the displacement resonance of the forced oscillations and Figure 17 shows the variation of the phase shift against the angular velocity, which is animated in Figure 18.

Figure 19 shows what happens with the mass-spring system when the damping constant and Figure 20 when the damping constant is approximately equal to zero, .

> display(seq(plt2(0.05*i),i=1..28),insequence=true,labels=[`omega`,`A`],labelfont=[TIMES,BOLD,12]);

Figure 16 Animation of Figure 15

> Harg:=i->plot(phi(i,omega),omega=1..3) :

> text3:=i->textplot([2.4,phi(i,2.4),convert(r=i,string)],align={RIGHT,ABOVE},font=[TIMES,BOLD,12]):

> plt3:=i->display(Harg(i),text3(i)):

> display(seq(plt3(i),i=[0.1,0.25,0.5,0.75,1,1.41]),labels=[`omega`,`phi`],labelfont=[TIMES,BOLD,12]);

Figure 17 Variation of the phase shift against the angular velocity of the forced oscillations-displacement resonance for different values of the damping constant

As the applied angular velocity is increased from zero to the natural angular velocity of the system, the phase angle increases from zero to . When the angle changes abruptly from to

> display(seq(plt3(0.05*i),i=1..28),labelfont=[TIMES,BOLD,12],insequence=true);

Figure 18 Animation of Figure 17

> f:=t->2*sin(2*t):

> mass_spring(1,0.25,4,0,5,0,0,20,f,60,scaling=constrained);

Figure 19 The displacement of the mass-spring system
undergoing forced oscillations with damping constant .

> mass_spring(1,0.01,4,0,5,0,0,11.5,f,60,scaling=constrained);

Figure 20 The displacement of the mass-spring system
undergoing forced oscillations with damping constant .

Summary

This paper discusses the use of the Computer Algebra System (CAS) Maple in calculating and animating the responses of free and force vibrations systems. CAS have changed the fundamental way in which many math and engineering course are taught. The user can enters a mathematical expression, which the CAS (in this paper Maple) subsequently execute and one are freed from tedious computations. This allow the student to explore how the solution of an engineering problem depends on various parameters of the problem under study. The graphics capabilities are extremely helpful for visualizing the behavior of the system under investigation. Maple's plots and plottools packages enable the user to utilize a lot of plotting objects for animating purposes. This is highlighted and demonstrated in this paper. It is hoped that this features can attract the attention of both lectures and students to visualize dynamic response, typical in engineering, more realistically.

>