Solving Second order ODEs for damping effect of a pendulum

AIM:

• To solve the 2^nd order ODEs for damping effect of pendulum

GOVERNING EQUATIONS:

The governing equation is the basic equation for the motion of a simple pendulum with damping

Where THETA is the angle suspended by string w.r.t axis passing through mean position of bob

b is the damping coefficient

m is the mass of bob

g is acceleration due to gravity

l is length of string to which bob is attached.

OBJECTIVES OF THE PROJECT:

• To solve the 2^nd order ODEs.
• To plot graphs of ‘angular displacement vs time’ and ‘angular velocity vs time’
• To produce an animation.

BODY OF THE CONTENT:

Damping is the process of dissipation of energy in an oscillating or vibrating body. There are different kinds of damping: -

1. Viscous damping- due to medium
2. Dry friction damping- due to friction between surfaces due to lack of lubrication
3. Hysteretic damping- due to friction between planes by slips or slide as the material deforms
4. Magnetic damping- due to magnetic effect

The pendulum's movement depends on the Newton’s second law. This law is written as a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

The code is produced to solve ODE for damping effect of pendulum. The primary quadratic equation can’t be solved directly. Hence we use an approach in which the equation is broken into two linear equations. It is done by creating a function with results d(theta1)/dt and d(theta2)/dt

The procedure consists of replacing the angle theta is given as theta1 and d(theta1)/dt is taken as theta2. Solving the equations, we get the two equations for 2^nd order ODE, as follows

1. d(theta1)/dt = theta2
2. d(theta2)/dt= -{(b/m) *theta2}-{(g/l) *sin(theta1)}

function [dtheta_dt]= dam_eqn(t,theta,g,l,b,m)

theta1=theta(1)
theta2=theta(2)

%breaking the given quadratic equation into two linear equation
dtheta1_dt = theta2;
dtheta2_dt = -b/m*theta2-g/l*sind(theta1);

%results
dtheta_dt =[dtheta1_dt; dtheta2_dt]
end

Further the solving of the equations will be done by the command “ode45”. This will give us separate value of angular displacement and angular velocity which were the results.

WORKING OF ODE45

clear all
close all
clc

%provided constants
b=0.05;
g= 9.81;
m= 1;
l= 1;

%initial conditions
int_con = [0,3];  %intial condition contains time and velocity respectively

%time span from ) to 20 sec 
time_span = linspace(0,20 ,150); 

%solving ode
[t,results]= ode45(@(t,theta)dam_eqn(t,theta,g,l,b,m),time_span, int_con);
 
figure(1)
plot(t,results(:,1))             %for angular displacement
line(t,results(:,2),'color','r') %for angular velocity
xlabel('timesteps');

%legends were added later by Open Property Inspector

ct=1                %counter
for i= 1: length(time_span)
    
    x1=0;           %initial point coordinates
    y1=0;               
    
    x2=l*sind(180+(results(i,1)));  %unfixed end of string's point coordinates
    y2=l*cosd(180+(results(i,1)));
    
    figure(2)
    plot(x2,y2,'o','markers',10,'markerfacecolor','r') %ball
    line([x1,x2],[y1,y2],'color','b')                  %the string
    line([-0.4 0.4],[0 y1])                 %the fixed base for swinging pendulum
    axis([-1 1 -2 2]);
    pause(0.03)
    M(ct)= getframe(gcf);             %to save the current frame
    ct=ct+1;                          %increment to counter
end
movie(M)
videofile=VideoWriter('pendulum.avi','Uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ode45 is a default code which uses Runge-kutta method. In this method. This method is used to solve first order odes when initial condition(f(tn,yn)) and step size(h) are provided. It uses the following 5 equations to solve the value of (n+1) th value where n is the number of time steps. The following equations are used to solve it

For the working of ode45, initial condition, time span and function are provided. This is given in following syntax.

[t,results]= ode45(@(t,theta)dam_eqn(t,theta,g,l,b,m),time_span, int_con)

Here, t and results on the left are solutions. The ode consists of t and theta as the variables. Function is called, further, time span and initial conditions are given.

Initial conditions - int_con = (0,3) [displacement, velocity]

Time span - time_span= linspace (0,20,150) [time is chosen between 0 sec to 20 sec]

Since the constants are provided, the calculation is carried out for 150 time steps. The results are denoted as angular displacement and angular velocity at different time steps.

ANIMATION

The input for i=1. The points are calculated for corresponding values provided earlier using the formulae and plot is formed. At this point frame is captured for recent counter(ct=1) and the further command provides the increment. The equations are to form the points

x1=0;

y1=0;

x2=l*sind(180+(results(i,1)))

y2=l*cosd(180+(results(i,1)))

The first point is the equation for fixed point of string. The later point (x2, y2) is the other end of string. To keep the pendulum in 3^rd and 4^th quadrants, 180deg were added to the angle for their trigonometric functions. Since the

Plot command is used to draw lines and circle for the animation.

plot(x2,y2,'o','markers',10,'markerfacecolor','r') [the bob of pendulum]

line([x1,x2],[y1,y2],'color','b') [the string]

line([-0.4 0.4],[0 y1]) [the fixed base for swinging pendulum]

Meanwhile, before every increment in counter, the frame is captured.

The captured frames are combined to form a movie with command “movie(M)”. further video is formed in uncompressed avi format.

RESULTS

The results are as follows: -

1. Graph of angular velocity and angular displacement vs the time step

The damping effect can be seen as the trough and crest position diminishes with time.

1. animation

The animation is formed in .avi format and is been attached here

https://drive.google.com/file/d/1Qw1ejKSATRkxA7GsRf3kou4F6dYCdS_Z/view?usp=sharing

CONCLUSION -

The pendulum's second order ODE was solved using ode45 command and command working was explained.

The plot for the angular momentum and angular displacement has been plotted and damping can be observed with time.

The animation gives a better insight of the situation.