Transient behavior analysis of a Simple Pendulum using MATLAB

AIM: To simulate the transient behavior of a simple pendulum using MATLAB.

GOVERNING EQUATION: The ODE (Ordinary Differential Equation) represents the equation of motion of a simple pendulum with damping as follows:

`frac{d^2 theta}{dt^2} + frac{b}{m}. frac{d theta}{dt} + frac{g}{L} .sin theta = 0` ... (1)

where;

`g` = Gravity (m/s²)

`L` = Length of the pendulum in (m)

`m` = Mass of the ball (kg)

`b` = Damping coefficient

`theta` = Angle of pendulum (radian)

`t` = time (sec)

OBJECTIVE:

1. To write a program for second-order ODE that describes the position of the ball with respect to time.

L = 1 m, m = 1 kg, b = 0.15, g = 9.81 m/s²

2. To plot the animation of simple pendulum motion for the conditions:

Between 0-20 sec, for angular displacement = 0, angular velocity = 3 rad/sec at time t=0.

THEORY: A simple pendulum may be described ideally as a point mass suspended by a massless string from some point about which it is allowed to swing back and forth in a place. A simple pendulum can be approximated by a small metal sphere that has a small radius and a large mass when compared relatively to the length and mass of the light string from which it is suspended. If a pendulum is set in motion so that is swings back and forth, its motion will be periodic. The time that it takes to make one complete oscillation is defined as the period T. The frequency f of the oscillations is the number of oscillations that occur per unit time and is the inverse of the period, f = 1/T. Similarly, the period is the inverse of the frequency, T = l/f. The maximum distance that the mass is displaced from its equilibrium position is defined as the amplitude of the oscillation.

When a simple pendulum is displaced from its equilibrium position, there will be a restoring force that moves the pendulum back towards its equilibrium position. As the motion of the pendulum carries it past the equilibrium position, the restoring force changes its direction so that it is still directed towards the equilibrium position.

SOLUTION:

Programming Language: MATLAB

Procedure/Steps:

1. Assign suitable variables to store the given input values.

2. ode45 is a function in Matlab that solves ordinary or partial differential equations using a numerical method.

[t, y] = ode45(odefun, tspan, y0)

where tspan = [t0 tf], integrates the system of differential equations y′ = f(t, y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

3. Figure 1 represents the plot for the graphical behavior of Angular Displacement and Angular Velocity.

4. Create a for loop for generation of frames and declare a variable named counter (ct) to start with.

5. The first loop (i) will be ranging from 1 to the maximum value of t_span.

6. For one count of loop i and compute the values (x1,y1,θ,Time) as shown in the code.

7. The images are plotted with the help of plot function.

8. Finally, the counter will increment by 1 (ct+1) until loop i is finished and step #4 and #5 are repeated.

9. getframe() captures the individual frame and stores it in an array M.

10. To create an animation of all the frames generated,

- assign the frame movie in array M

- assign appropriate file name and extension in videofile()

- open() begins to extract the files

- writeVideo writes data from an array to the video file

- close() terminates the operation

Code: The MATLAB program to simulate a simple pendulum is shown below.

% Program to plot transient behavior of a simple pendulum by solving second-order ODE.

% Inputs:

b = 0.15;  % Damping coefficient
g = 9.81;  % Acceleration due to gravity (m/s^2)
L = 1;     % Length of the pendulum (m)
m = 1;     % Mass of the ball (kg)

% Initial condition
theta_x = [0;3];   % Angular Displacement
t_span = linspace(0,20,400);   % Time span


% Solution:

% ODE45 Function
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,L,m), t_span, theta_x);


% Plot 1:

figure(1)
plot(t,results(:,1),'linewidth',2)
hold on
plot(t,results(:,2),'linewidth',2,'linestyle','--')
ylabel('Angular Displacement (rad) / Angular Velocity (rad/sec)')
xlabel('Time (sec)')
title('Transient Behavior of a Simple Pendulum')
legend('Angular Displacement','Angular Velocity')
grid on

ct=1;
for i=1:length(t_span)
    j = results(i,1);
    x0 = 0;
    y0 = 0;
    x1 = L*sin(j);
    y1 = -L*cos(j);
    
    % Assigning the parameters
    txt = ['   θ = ' , num2str(results(i,1)) , ' deg'];
    txt1 = ['Frictionless Pivot'];
    txt2 = ['   Time = ' , num2str(t(i,1)) , ' sec'];
    txt3 = ['  Massless String'];
    txt4 = ['     x1 = ' , num2str(x1) , ' , ' , 'y1 = ' , num2str(y1)];
    txt5 = ['Bob'];
    

    % Plot 2:

    figure(2)
    plot([-1 1], [0 0],'linewidth',6,'color','g')
    hold on
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b')
    plot([0 0],[0 -1],'linestyle','--','linewidth',1,'color','k')
    plot(x1,y1,'marker','o','markersize',15,'markerfacecolor','r')
    text(x0,y0,txt,'VerticalAlignment','top')
    text(x0-1,y0,txt1)
    text(-0.4,0.3,txt2)
    text(0.5*(x0+x1),0.5*(y0+y1),txt3)
    text(x1,y1,txt4)
    text(x1-0.4,y1,txt5)
    hold off
    xlabel('X-axis (m)')
    ylabel('Y-axis (m)')
    title('Transient Behavior of a Simple Pendulum')

    grid on
    axis([-2 2 -1.5 0.5])
    pause(0.01)
    
    M(ct) = getframe(gcf);
    ct = ct+1;
end


% Creating the animation:

movie(M);
videofile = VideoWriter('Simple Pendulum Motion.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

Function: The MATLAB program as a matrix function [dheta_dt] is shown below. We call the function ode_func to find the solution as a function of t,theta,b,g,L,m.

Converting second-order ODE into first-order ODE to solve in MATLAB.

`frac{d^2 theta}{dt^2} + frac{b}{m}. frac{d theta}{dt} + frac{g}{L} .sin theta = 0` ... (1)

Let `theta = theta_1`

`:. frac{d theta} dt = frac{d theta_1}{dt} = theta_2`

`:. frac{d^2 theta} {dt^2} = frac{d^2 theta_1}{dt^2} = frac{d theta_2}{dt}`

Put in equation (1);

`:. frac{d theta_2}{dt} + frac{b}{m}.theta_2 + frac{g}{L}. sin theta_1 = 0`

`:. frac{d theta_2}{dt} = - frac{b}{m}.theta_2 - frac{g}{L}. sin theta_1` ... (2)

Also `frac{d theta_1}{dt} = theta_2` ... (3)

In matrix form;

`frac{d}{dt} [(theta_1),(theta_2)] = [(theta_2),(-frac{b}{m}.theta_2 - frac{g}{L}.sin theta_1)]`

function [dtheta_dt] = ode_func(t,theta,b,g,L,m)

% Inputs:

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


% Solution:

dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

Output: The program is executed as per given parameters as shown below.

Plot: Angular Displacement (━) and Angular Velocity (- -) vary with respect to the time as shown below.

Animation:

CONCLUSION: From the above results, we can infer that:

- Angular Displacement and Angular Velocity converges in a periodic fashion as a function of time due to the presence of damper in the system.

Thus, the program for a simple pendulum was coded and the position of the hanging ball (bob) was simulated in MATLAB plot as a video.

REFERENCES:

Ode45 Function