Analysis of damped motion of a simple pendulum

Aim of the project:- Analyze the damped motion of a simple pendulum based on 2nd order ordinary differential equations and create a program to compute the solution for the case study based on the same phenomenon in matlab.

Theory:- A simple pendulum is an object suspended in a fixed point so it can swing freely. A simple pendulum can be seen below,

The assumptions of the case study of simple pendulum motion are detailed below:

1. The object is joined with a fix point through a non-extending wire and damping is caused by any viscous media such as air resistance.

2. The weight of the wire is negligible in comparison with the weight of the bob and the wire always remains taut.

3. The motion occurs in a 2-dimensional plane, the object does not trace an ellipse.

When this object is displaced from its resting equilibrium position, it is subjected to a restoring force due to the gravity that accelerates it back towards the equilibrium position. When released, the restoring force will make it oscillate around the equilibrium position, swinging back and forth.

All the objects in movement that are in contact with any fluid, suffer a force. This force is called drag (fluid resistance) and this force acts in the opposite direction of the movement. Unlike other resistive forces such as dry friction, which is nearly independent of velocity, drag forces depend on velocity.

Analysis of the simple damped pendulum motion:-

The general equation of the force on a moving object due to a fluid is given by,

`F_d = 1/2*C_d*rho_f*A*V^2`,

where `rho_f` is the density of the fluid, `A` is the projected area of the object in the fluid flow, `V` is the velocity of the object and `C_d` is the drag coefficient. The drag coefficient is a parameter that depends on the Reynolds number. Reynolds equation and the equation of the drag coefficient can be seen below,

`Re = (rho_f*V*L_c)/mu` and `C_d ~~ 24/(Re)+6/(1+sqrt(Re))+0.4`,

where `mu` is the dynamic viscosity of the fluid and `L_c` is the characteristic length of the object. The Reynolds number is an important number to determine the behaviour of the fluid and especially the transition between the laminar flow and the turbulent flow.

In this case study it is assumed that there is laminar flow for still air flow i.e. reynolds number is very small. Now,

`F_d = 1/2*24/(Re)*rho_f*A*V^2`

Substituting the value of `Re` in the above equation,

`F_d = 1/2*24/((rho_f*V*L_c)/mu)*rho_f*A*V^2`

or, `F_d = 12*mu*A/L_c*V`

This equation is the stokes equation if the object is spherical. This equation is appropriate for objects moving through a fluid at relatively slow speeds where there is no turbulence and for this reason this equation can be used.

In this case study and with the purpose of simplifying the equations, the expression `12*mu*A/L_c` will be called as `b` (drag co-efficient for a laminar flow as the damping factor). The final equation will be framed as,

`F_d = b*V` , the drag force acts in the opposite direction of the movement of the simple pendulum`

Variation of the velocity, acceleration and tension in the wire along the motion of the simple pendulum can be seen as below,

The new force diagram taking the drag force into account can be seen as below,

Using the second law of Newton yields,

`m*a_n = W - mg*costheta` (along normal direction)

`m*a_t = - mg*sintheta - b*V` (tangential direction)

The equation for obtaining the motion of the pendulum is the equation of the tangential direction. Now substituting the linear velocity and the tangential acceleration in the above equation for its angular expression and writing that in a differential form,

`m*alpha*L = -mg*sintheta - b*omega*L`

or `m*L*(d^2(theta))/(dt^2) = -mg*sintheta - b*L*(d(theta))/dt`

or `(d^2(theta))/(dt^2) + (mg*sintheta)/(m*L) + (b*L)/(m*L)*(d(theta))/dt = 0`

And the final equation for the damped motion of a simple pendulum to be used for matlab programming can be written as,

`(d^2(theta))/(dt^2) + (g)/(L)*sintheta + (b)/(m)*(d(theta))/dt = 0`

This is a case of 2nd order Oridinary differential equation (ODE) as the power of the derivative in the equation is 2. An ODE is an equation that contains one independent variable (e.g. time) and one or more derivatives with respect to that independent variable. In the time domain, ODEs are initial-value problems, so all the conditions are specified at the initial time t = 0.

Matlab has several different functions (built-ins) for the numerical solution of ODEs. These solvers can be used with the following general syntax:

Matlab has different solvers available for different application problem statements,

In general defining an ODE function in an matlab file (m-file) is as,

[t,state] = ode45(@dstate,tspan,ICs,options)

1. Define tspan, ICs and options in one file (e.g. call_dstate.m), which sets up ode45

2. Define your constants and derivatives in another file (e.g. dstate.m) or as a function dstate within the call file

3. Run call_dstate.m

4. Analyze the results as a plot of state vs. t

For problems involving higher order ODE such as 2nd order ODE in this case study of simple damped pendulum motion, the higher order ODE should be divided in number of first order ODEs.

The 2nd order ODE for this case study can be re-written in two first order ODEs as,

`(d^2(theta))/(dt^2) + (g)/(L)*sintheta + (b)/(m)*(d(theta))/dt = 0` (2nd order ODE equation)

Let `theta = theta_1`

or `(d(theta))/dt = (d(theta_1))/dt = theta_2 (say)`

or `(d^2theta)/dt = (d^2theta_1)/dt = d/dt((d(theta_1))/dt) = (d(theta_2))/dt`

Now,

`(d^2(theta))/(dt^2) + (b)/(m)*(d(theta))/dt + (g)/(L)*sintheta = 0` (this equation is a non-linear equation because of the presence of `sintheta` term)

or `(d^2(theta_1))/(dt^2) + (b)/(m)*(d(theta_1))/dt + (g)/(L)*sintheta_1 = 0`

or `(d(theta_2))/(dt) + (b)/(m)*theta_2 + (g)/(L)*sintheta_1 = 0`

The 2nd oder ODE can be written in two first order ODE as,

`(d(theta_1))/dt = theta_2` and `(d(theta_2))/dt = -b/m*theta_2 - g/L*sintheta_1`

Writing the above first order ODEs in matrix form as general form of every entity in matlab is in matrix format,

`d/dt [(theta_1), (theta_2)] = [(theta_2), (-b/m*theta_2 - g/L*sintheta_1)]` , where `d/dt[(theta_1), (theta_2)]` is the solution vector for the ODE.

Objectives of the project:-

1. To analyse the case study problem statement and to create the program that solves the 2nd order ODE equation.

2. Plotting the parameters of the damping pendulum motion w.r.t time.

3. Create a animation movie of the simple damped pendulum motion.

4. Create a animation movie of the pendulum motion along with the damping parameters of the pendulum motion in real time.

Tools Used:- Matlab

Case Study:- The simple damped pendulum motion is given by `(d^2(theta))/(dt^2) + (b)/(m)*(d(theta))/dt + (g)/(L)*sintheta = 0` where `g` is gravity in `m/sec^2` i.e. 9.81 `m/sec^2`, `L` is length of string of the pendulum in metre i.e. 1 metre, `m` is mass of the ball of the pendulum in Kg i.e. 1 Kg, `b` is damping co-efficient i.e. 0.05. Conditions for the case study are:- stimulation time for the pendulum motion is 0 - 20 secs, initial angular displacement of the pendulum when at rest is zero, initial angular velocity is 3 rad/sec at time t=0.

The matlab program for the same can be written as,

In the program the parameters are defined and the ode45 solver is used for solving the ODE in the given case study and in the ode45 solver, a function call is present that intimates the function definition with the initial set parameters to use as inputs to solve the ODE. The function definition of the following matlab program can be written as below,

After the numerical computation is done to solve the ODE, the plot command for the effect of damping on angular velocity and angle of displacement can be written as,

The output of the code for plotting can be seen as below,

The animation of the damping effect on a simple pendulum motion can be created by writing the code in matlab as shown below,

The movie of the animation of the simple damped pendulum motion can be created by writing the code in matlab as below,

The output is the movie of the animation,

The snapshot of the workspace for reference of variables can be seen as below,

The code snippet for the entire matlab program can be seen as below,

%Aim of the program is to solve 2nd order differential equation of a simple pendulum problem

close all 
clear all
clc

b =0.05;     %Damping Co-efficient
g = 9.81;    %acceleration due to gravity, m/sec^2
l = 1;       %length of the string of the pendulum is 1 metre
m = 0.1;     %mass of the bob of the simple pendulum is 1 kg

%Initial condition
%Initial values for the 2 first order ode i.e initial condition for
%the first order odes are initial displacement as 0 and initial angular velocity 
%as 3 rad/sec
theta_0 = [0;3];

%time points 
t_span = linspace(0,20,200);    %for 20 secs of time span

%solve ODE
[t, results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

omega_n = sqrt(g/l);
T = 2*pi/omega_n;

figure(1)   %Plot of the damping oscillation of the pendulum 
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Magitude')
xlabel('time (s)')
legend('angle of displacement (rad)','angular velocity (rad/sec)')
grid on

ct=1;
figure(2)
for i = 1:length(results(:,1)) %for loop is excuted till number of entries in results variable
    
    %Origin at the fixed end of string i.e the first point of start of length of
    %string
    s_x0 = 0;
    s_y0 = 0;
    
    %Every time the loop executes the incremented value of the angle of
    %displacemnt is stored in a_d to locate the other plot points
    a_d = results(i,1);
    
    %End point of length of string for every amgle of displacement 
    s_x1 = l*sin(a_d);
    s_y1 = -l*cos(a_d);
    
    %Plot for the fixed end
    %Starting point of fixed end 
    f_x0 = -1;
    f_y0 = 0;
    %Ending point of the fixed end 
    f_x1 = 1;
    f_y1 = 0;
    
    plot([f_x0, f_x1],[f_y0, f_y1],'color','black','linewidth',20)
    axis([-1, 1, -1.5, 1]) %axis limits of the plot
    hold on
    
    %Fixed point of pendulum
    plot(0,0,'-o','markers',8,'markerfacecolor','white')
    axis([-1, 1, -1.5, 1]) %axis limits of the plot
    hold on
    
    %Plot for the string 
    plot([s_x0, s_x1],[s_y0, s_y1],'color','b','linewidth',3)
    axis([-1, 1, -1.5, 1]) %axis limits of the plot
    hold on
    
    %Plot for the bob of the pendulum at the end 
    %creating a circle at the end of the length of the string 
    plot(s_x1,s_y1,'-o','markers',30,'markerfacecolor','r') 
    axis([-1, 1, -1.5, 1]) %axis limits of the plot
    text(-0.38,0.75,'Damped pendulum motion')
    title([num2str(t(i)),'s'])
    pause(0.01)  %time interval to display every plot
    hold off   %update every cycle individually
    M(ct) = getframe(gcf);
    ct = ct+1;

end


movie(M)    %creates a movie from all the frames captured
videofile = VideoWriter('pendulum_motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

The animation movie of the pendulum motion along with the damping parameters of the pendulum motion in real time can be seen for a better visual of the case study, the output of the matlab program for the same yields the following output,

Conclusion:- Similarly, there are hundreds of ODEs that describe key phenomena and the ability to solve these equations will make oneself able to simulate any of these systems. The objectives mentioned at the start of this project are carried out and the analysis of the damped motion of the simple pendulum is done, which yielded results same as compared to natural occuriing phenomenon.