Simulation of Damped Pendulum using MATLAB

Abstract :

A differential equation (D.E.) is an equation with a function and one or more of its derivatives or differential coefficients or differentials. An ordinary differential equations (ODE) is that in which all the derivatives w.r.t. a single independent variable. Given equations are examples of ODE-

`frac{d^2x}{dt^2} + a^2x = 0` (Second Order ODE)`

`x^2dx + y^2dy = 0` (First Order OD)`

In this project we aim to solve a second order ODE governing a damped pendulum and simulate the motion for 20 sec. The mass of bob, length of string, damping coefficient and initial conditions which are angular velocity(3 rad/sec) and angular displacement(0 radians) at time t=0sec are given in our problem statement. We will be using 'ode45' in-built function in MATLAB to solve our linear differential equation.

The difference between order and degree of a D.E.

The order of D.E. is the order of the highest derivative occurring in it whereas the degree of a D.E. is the degree of the highest order derivative occuring in it, after the equation has been expressed in a form free from radicals and fractions as far as derivative is concerned. It will be better explained with some example -

1. `yfrac{dy}{dx} = x(frac{dy}{dx})^2 + x^2`

The above D.E. is of first order (`frac{dy}{dx}`) but second degree (the power of highest order derivative).

2. `[1+(frac{dy}{dx})^2]^-5=a^3(frac{d^2y}{dx^2})^3`

The above D.E. is of second order (`frac{d^2y}{dx^2}`) and third degree (the power of highest order derivative).

Governing Equation for a damped pendulum:

The figure above is a FBD diagram of a simple pendulam. The weight of the pendulum acts as a restoring torque about the fixed(hinged) end. The most of the real life systems has some kind of friction or resisting force which reduces the system energy or we can say dampens the motion. The system slowly comes to a halt after dissipating all its energy. To depict such system we use a damper and thus called as damped pendulum.

The equation is as followed - `frac{d^2theta}{dt}+frac{b}{m}frac{d theta}{dt} + frac{g}{L}sintheta = 0`

`theta` = angular displacement from equilibrium position

b = damping coefficient (N-s/m) **Damping force is directly proportional to velocity of the mass

m = mass of the bob

g = acceleration due to gravity

L = length of string

Solving ODE in MATLAB

MATLAB solves only first order D.E. and therefore in order to solve higher order D.Es. we have to split them into different first order D.E. Taking governing equation into consideration following steps are followed-

1. Let `theta = theta_1 => frac{d theta}{dt}= frac{d theta_1}{dt} = theta_2` - (i) This equation is our first D.E. for 'Displacement`
2. Let `frac{d theta_1}{dt}=theta_2=> frac{d^2theta}{dt} = frac{d^2theta_1}{dt}=frac{d theta_2}{dt}`
3. Now substitute equation from 2. in our governing equation. We will get-`frac{d theta_2}{dt}= -frac{b}{m}theta_2-frac{g}{L}sin theta_1` - (ii) This equation is our second DE for 'Velocity'

From (i) and (ii) we will generate our matlab function named "myodefunction". We will convert these two equation in matrix form-

`frac{d}{dt}[[theta_1],[theta_2]]=[[theta_2],[-frac{b}{m}theta_2-frac{g}{L}sintheta_1]]`

Left hand side can be clubed and termed as `frac{d [theta]}{dt}` - (iii)

1. Code for generating a function named 'myodefunction'

function [dtheta_dt] = myodefunction(t,theta,b,l,m,g)
    
  theta1 = theta(1);   %theta1 is the ang. position of pendulam from mean position(in rad)
  theta2 = theta(2);   %theta2 is derivative of theta1 and becomes ang.velocity(in rad/sec)
    
  dtheta1_dt = theta2;                               %1st diff. equation
  dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);    %2nd diff. equation
    
  dtheta_dt = [dtheta1_dt ; dtheta2_dt];
    
end

The above function has 6 input arguments that generates a vector named 'dtheta_dt'. 'theta' i.e. `theta`in our equation (iii) is a vector with 2 elements.

NOTE. Syntax for using ode45

`[tsol,ysol]=ode45((At sign)(t,y) fname(t,y),[t_o t_(end)],y_o`

tsol, ysol : solution vector; Matlab returns ysol for each tsol

fname : function that returns dy = F(t,y)

`t_o,y_o` : initial conditions representing `y(t_o) = y_o`

`t _(end)` : final value until which solution is desired. If we put a vector of a 1xn dimension instead of `[t_o t_(end)]` , then we will receive result for each value of that vector.

2. Code for main body for solving and ploting the simulation

%Defining given geometric parameters
b = 0.05;             %damping coefficient (N-s/m)
l = 1;                %length of pendulam
m = 1;                %mass of bob
g = 9.81;             %acceleration due to gravity

%Initial Conditions
theta_0 = [0;3];      %1st element is for the displacement that is zero at time= 0
                      %2nd element is for angluar velocity at zero displacement

%time span for which solution is calculated
t_span = linspace(0,20,500);

%solve ode using ode45
[time, solution] = ode45(@(t,theta)myodefunction(t,theta,b,l,m,g), t_span,theta_0);

%plotting the curves
figure(1)
plot(time, solution(:,1))
hold on
plot(time, solution(:,2))

%PROGRAM FOR PENDULAM SIMULATION

%fixed position of the pendulam
x0 = 0; y0 = 0;
%Using for loop to generate plot at different angular position
ct = 1;        
for i = 1:length(solution)
    
    ang_disp = solution(i,1); %angular displacement of bob
    
    %bob position for the angular dispalce
    x1 = l*sin(ang_disp);
    y1 = -l*cos(ang_disp);
    
    %ploting the simulation
    figure(2)
    plot([x0 x1],[y0 y1],'LineWidth',3,"Color",'r')
    hold on
    plot(x1,y1,'-o','markersize',20,'markeredgecolor','b','markerfacecolor',"y");
    grid on
    axis([-2.5 2.5 -2 2])
    line([-2 2],[0 0],'Linewidth',5,"Color",'g')
    pause(0.12)
    hold off
end

RESULT :

1. Plot of ang.displacement and ang. velocity w.r.t time[0-20]

We can observe in the plot that as time increases the velocity and angular position are decreasing because of dampening effect.

2. Link of animation of the pendulam

https://youtu.be/t8XXpjJe-2k

Image of pendulum at one of the value of `theta`

Explaination of above code-

1. Under comment %Defining given geometric parameters

We assaign values for different variales given in our problem statement.

2. Under comment %Initial Conditions

We generate a vector for our inital given conditions.

3. Under comment %time span for which solution is calculated

We divide our 20 sec time frame into 500 equispaced parts using 'linspace'.

4. Under comment %solve ode using ode45

We will use ode45 to solve our differential equation. This returns two different arrays i.e. time(500x1) and solution(500x2). The time array contains all values of time for which solution is found. The solution array 1st column contains values of displacement(`theta`) and 2nd column contains values of velocity(`frac{d theta}{dt}`) for each value of time.

5. Under comment %plotting the curves

The plot is given above.

6. Under comment %PROGRAM FOR PENDULAM SIMULATION

Firs we will fix our hinge joint at (0,0). Next, We generate a for loop to plot the position of the string along with bob at each value of angular displacement at the given time. While plotting the order should be kept in mind. Otherwise the images will overlap over each other.

CONCLUSION :

From plot given above it can be easily concluded that our damped system will slowly loose its velocity and displacement as time passes. The more we increase damping coefficient faster will be the decrease.

REFERENCES :

https://skill-lync.freshdesk.com/support/solutions/articles/43000524980-hold-on-and-hold-off-command