Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Finite Difference Approach by MATLAB for the First and Second Derivatives

The following MATLAB program determines the first and second derivatives of the data given in the problem applying the finite difference schemes and developing a custom user defined function firstsecondderivatives(x,y).

The finite difference schemes used are the following:

For first derivative:
  • At first point (three point forward)
  • At last point (three point backward)
  • At all other points (two point central)

For second derivative:
  • At first point (four point forward)
  • At last point (four point backward)
  • At all other points (three point central)

DATA SET
 x -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
 f(x) -3.632 -0.3935 1 0.6487 -1.282 -4.518 -8.611 -12.82 -15.91 -15.88 -9.402 9.017


function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.

% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n = length (x);
dy = zeros;
ddy = zeros;

% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.

% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.

dy(1) = (-3*y(1) + 4*y(2) - y(3)) / 2*(x(2) - x(1)); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

for i = 2:n-1
   
dy(i) = (y(i+1) - y(i-1)) / 2*(x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;

end

dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / 2*(x(n) - x(n-1));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;

figure(1)
subplot (3,1,1)
plot (x, y, 'linewidth', 2, 'color', 'k'), grid
title('x VS f(x)')
ylabel ('f(x)')

subplot (3,1,2)
plot (x, dy, 'linewidth', 2, 'color', 'b'), grid
title('x VS dy')
ylabel ('dy: First Derivative')

subplot (3,1,3)
plot (x, ddy, 'linewidth', 2, 'color', 'r'), grid
title('x VS dy')
ylabel ('ddy: Second Derivative'), xlabel('x')

figure(2)
plot (x, y, 'linewidth', 2, 'color', 'k');
hold on;
plot (x, dy, 'linewidth', 2, 'color', 'b');
hold on;
plot (x, ddy, 'linewidth', 2, 'color', 'r'), legend
ylabel ('y, dy, ddy'), xlabel('x')
end


The following operations are done in MATLAB command window to run the above function.

Program Outputs:

>> x = -1:0.5:4.5;
>> y=[-3.632 -0.3935 1 0.6487 -1.282 -4.518 -8.611 -12.82 -15.91 -15.88 -9.402 9.017];
>> firstsecondderivatives(x,y)

ans =

  Columns 1 through 11

    2.0805    2.3160    0.5211   -1.1410   -2.5833   -3.6645   -4.1510   -3.6495   -1.5300    3.2540   12.4485

  Column 12

   12.1947

We can also plot the derivatives, below is the figure for the function, first and second order derivatives respectively.

Plotting the first and second order numerical derivatives
Plotting the first and second order numerical derivatives

Learning Mathematica, Lesson 2: Solving Euler-Bernoulli Beam Equation

In this lesson, I would like to show the advantages of the Mathematica built-in solver to evaluate the analytical solution of a differential equation. For example, if we want to solve the well-known fourth order Euler-Bernoulli equation to solve a problem of a cantilever beam, the Mathematica code has very user-friendly features to do so. In following figure, we see that the beam has one fixed end and the other end has a concentrated load acting downward.

The following Mathematica program determines the analytical displacement function and plots the result. Most of the lines are commented as explanations for the line of codes to understand the respective operation.

Cantilever beam with a fixed end and concentrated load at free end
Cantilever beam with fixed end


(*Determination of the Exact Solution of an Euler-Bernoulli Beam*)
Clear[y, P, x, EI, L]
y[x]; (*The Final Solution, y*)
(*Defining the constant parameters, just assuming unit values*)
h = 1; (*Beam height*)
P = 1; (*Concentrated load at the end*)
L = 10; (*Beam length*)
b = 1; (*Beam width*)
Ie = 1/ 12; (*Moment of inertia*)
Ee = 10 000; (*Young's Modulus*)
EI = Ee* Ie; (*Flexural Modulus*)
th = D[y[x], x]; (*Slope*)
V = EI* D[y[x], {x, 3}]; (*Shear force*)
M = EI* D[y[x], {x, 2}]; (*Bending moment*)
(*Defining the boundaries of the problem*)
y1 = y[x] /. x → 0;
y2 = y[x] /. x → L;
th1 = th /. x → 0;
th2 = th /. x → L;
M1 = M /. x → 0;
M2 = M /. x → L;
V1 = V /. x → 0;
V2 = V /. x → L;
(*Using DSolve a built-in solver to solve beam equation*)
s = DSolve[{EI* y''''[x] ⩵ 0, y1 ⩵ 0, th1 ⩵ 0, V2 ⩵ P, M2 ⩵ 0}, y, x];
Print["The Exact Displacement Function"]
y = Simplify[y[x] /. s[[1]]]
Plot[y, {x, 0, 5}, PlotStyle → {Green, Thick}]
Print["Displacement at the tip of the beam in Y direction"]


Mathematica generates the analytical solution function which is seen below. This feature is very useful as you may solve any differential equations analytically and get the results symbolically as well like this problem.


analytical solution function

Linear Least Squares Regression Analysis by a MATLAB program

A MATLAB program is developed to determine the coefficients by linear least squares regression where the function is, y = mx + b. Here,


x
-6
-4
-1
0
3
5
8
y
18
13
6
4
-1
-8
-15


% A function defined for least square regression
function [c,R2] = linearregression(x,y)

% Least-squares fit of data to y = c(1)*x + c(2)
% Here, c(1) = m; c(2) = b;
% Inputs:
% x,y = Vectors of independent and dependent variables
% Outputs:
% c = Coefficients of the given equation

if length(y)~= length(x)
    error("x and y are not compatible");
end

% Least squares algorithm implementation
x = x(:); y = y(:); % x and y are column vectors
A = [x ones(size(x))]; % m-by-n matrix of the system
c = (A'*A)\(A'*y); % Solving normal equations

if nargout>1 % Checking number of function outputs
r = y - A*c;
R2 = 1 - (norm(r)/norm(y-mean(y)))^2;
end


Program Output:

>> x = [-6 -4 -1 0 3 5 8];
>> y = [18 13 6 4 -1 -8 -15];
>> linearregression(x,y)

ans =

   -2.3140

    4.0814

A MATLAB Program to Implement the Jacobi Iteration

A MATLAB Program to Implement Jacobi Iteration to Solve System of Linear Equations:

The following MATLAB codes uses Jacobi iteration formula to solve any system of linear equations where the coefficient matrix is diagonally dominant to achieve desired convergence.

%-----------------------------------------------------------------%
function [x, rel_error] = jacobimethod(A, b, x0, tol, iteration)

% Inputs: A - Coefficient matrix
%         b - Input matrix
%       tol - Defining tolerance for solution
% iteration - Number of iterations

% Outputs: x - Solutions
%  rel_error - Relative error

D = diag(diag(A));   % Making coefficient matrix diagonal
R = A - D;           % Construction of another matrix "R"
N = 1;               % iteration counter
x = x0;
rel_error = tol * 2; % norm(x - x0)/norm(x);

% Implementation of Jacobi method to solve Ax = b
while (rel_error>tol && N <= iteration)
   
    xprev = x;   
    x = inv(D)*(b - R*xprev);
    rel_error = norm(x - xprev)/norm(x);

    Z1=x(1);
    disp(Z1);
    Z2=x(10);
    disp(Z2);
    Z3=x(16);
    disp(Z3);
   
    fprintf('\n Iteration %i: Relative error =%d ',N, rel_error)   
    N = N + 1;       
end
%----------------------------------------------------------------%

Now, define your input matrices in the MATLAB command window and use the above developed function to compute Jacobi iteration.

MATLAB Command Prompt:

>> A = [-2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2];

>> b = [-1; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; 1];

>> x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];

>> tol=10^-5;

>> iteration=350;

>> jacobimethod(A, b, x0, tol, iteration)


ans =

   14.8499
   28.7009
   40.5541
   50.4104
   58.2708
   64.1360
   68.0066
   69.8830
   69.7653
   67.6537
   63.5477
   57.4473
   49.3516
   39.2600
   27.1715
   13.0852

A MATLAB Program to Determine the Roots of Equation by Secant method

Problem: Solve the following problem by Secant Method.

f(x) = x - 2e^-x

Determine the roots of the above function where x1 = 0 and x2 = 1.


Solution: The following MATLAB program solves this problem.

% Definition of a function "secantmethod" to solve equation by Secant Method

function [x,ea] = secantmethod(X,X0,etol)

format long;

% Input: X=0, etol=Tolerance definition for error
% Output x=Solution of equation, ea=Calculated error in loop

% Program Initialization

% Iterative Calculations

while (1)
   
    % Implementation of Secant Method
         
    solution = X-((X-(2*exp(-X)))*(X0-X))/((X0-(2*exp(-X0)))-(X-(2*exp(-X))));    
   
    solutionprevious=X;
    X0=X;
    X=solution;
                 
if (solution-solutionprevious)~=0          % Approximate percent relative error
                                
   ea=abs((solution - solutionprevious)/solution)*100;
      
end

if ea<=etol                 
  
    break,               
end

%x0=solutionprevious;
x = solution;

%disp(x0);
disp(x);

end

% Display of output parameters

disp(x);                     

disp(ea);

%disp(solution-(0.1*abs(ea)*solution));

end


Program Outputs:

>> secantmethod(1,0,1e-10)

   0.883298154248460

   0.851584208519236

   0.852612692886131

   0.852605503703827

   0.852605502013723

   0.852605502013723

     3.125167797776517e-13


ans =

   0.852605502013723


#SecantMethod

My First Code in Python: Adding all Prime Numbers

As Python programming is highly hyped now-a-days and a frequent buzzword in machine learning and AI communities, I have decided to get myself familiar in this area. As I learn and decipher the Python coding day-by-day, I will be more than happy to share my coding in this blog. Here is my first program in Python - to develop a logic that adds all prime numbers between 0 and a given number.

We know that a prime number is an integer, which is greater than one and which is only divisible by one and itself. The following shows an algorithm that adds all prime numbers between 0 and a given number.

Development of an algorithm to add all prime numbers between 0 to n

I. First, we define the lower and upper number to initialize the program. For example, we set the lower number 0 and upper number 10; and we are interested to find the summation of prime numbers in this range.
II. Prime number is greater than one, so our logic begins that it must be greater than one. Therefore, we begin with number 2 and list all the numbers in our range.
III. Since, prime numbers are divisible only by itself and one; we create a loop, which confirms the prime numbers, and display them. We implement this by checking if the remainder of the number is zero or not. If the remainder is zero, then it is not a prime number.
IV. If the remainder of the number is not zero, then it is a prime number.
V. Display the prime numbers in the defined range and group the numbers in a list or array.
VI. Finally, we add all the prime numbers in the list.

The following Python code is developed based on the above algorithm, which displays all the prime numbers in a given range and add those numbers:


lower = 0

upper = 10     #Defining the range

sum1 = 0       #Setting Summation of prime numbers to zero

print("Prime numbers between", lower, "and", upper, "are:")

for num in range(lower, upper + 1):

   if num > 1:                 #Prime number is greater than 1

       for i in range(2, num):

           if (num % i) == 0:  #Calculates & checks remainder of division

               break
       else:
         
           print(num)  # Displays prime numbers
         
           sum1 += num
         
print(sum1)


Program Outputs:

Prime numbers between 0 and 10 are:
2
3
5
7
17




#Blog #Blogger #Python #PrimeNumber #Algorithm

Looking for an Industrial Partner to Implement a Lemon Shaped Guide to Minimize Rubbing

A NOVEL LEMON SHAPED GUIDE TO MINIMIZE EXCESSIVE RUBBING IN ROTATING MACHINES

In industries, rotating machines are widely used, because rotation offers a great way to transfer power from one point to another and convert motion to different planes by gears, belts, shaft etc. A rotating machine typically includes a rotor, bearings and a support structure. There are critical relations among these components where each component of a system influences the overall dynamic behavior of a machine. For example, rotor-to-guide rub degrades a mechanical system over the years and may even cause fatal accidents earlier. It is, therefore, paramount in industries, to run rotating machines, operating at high speeds smoothly and reliably. The primary reasons of rubbing between a rotor and a guide are due to a manufacturing error, excessive imbalance, misalignment, bearing wear, smaller radial clearance between the rotating shaft and casing, bad assembly, etc. Rub occurs in rotor casing, seals, unlubricated journal bearings, loose rotor guide attached to restrict a large deflection. The problem is prevalent in the industry and demonstrated in several literatures. To deal with this problem, the industry has already been using circular shaped guide to minimize excessive vibration between a rotor and a stator. Although, the circular shaped guide may reduce the vibration, but the rub, between rotor and guide, may still be present, which may lead to the permanent damage of a mechanical system. The lemon shaped guide is not only effective in minimizing rubbing between rotor and stator, but also it suppresses the excessive vibration similar to the circular shaped guide.

The following video shows an experiment on how to minimize rubbing between a rotor and a guide that typically happens in a high speed rotating machine with a newly developed lemon shaped bearing or guide. We see here as the speed goes high, the system enters into its resonance frequency where it vibrates excessively. The lemon/elliptical shaped guide helps prevent the excessive rubbing while the rotor-bearing assembly is in the natural frequency zone. It has been found from our research that the lemon shaped guide minimizes rubbing better than the circular shaped guide where there the rotational speed is very high (https://doi.org/10.1115/1.4043817). 



The next video shows an experiment on how to minimize rubbing between a rotor and a guide that typically happens in a high speed rotating machine with a circular bearing or guide. This is a traditional approach, which has been in operation in industries for long. We see in the video, as the speed goes high, the system enters into its resonance frequency where it vibrates excessively. The circular guide helps prevent the excessive vibration while the rotor-bearing assembly is in the natural frequency zone. However, if we notice carefully, the rubbing is still present between the rotor and guide, which deteriorates the system gradually.


A US provisional utility patent (US Patent App. 62/956,833) has already been filed for this design. I am looking for industrial partners or investors if they are interested in investing this product.

Learning Mathematica, Lesson 1: Plotting a Function

This is the very first lesson or tutorial on Mathematica. As I am in the process of learning this tool, I will gradually post more articles on this, ranging from basic to advanced level problems. So, the first question is, what is Mathematica? It is simply a tool for computing, but it has an advantage that the symbolic expression is much user-friendly and more interactive than MATLAB.

Although, MATLAB is a much bigger platform than Mathematica, because MATLAB has numerous toolboxes and libraries that are designed for specific fields. Nevertheless, Mathematica, a product from Wolfram Research, is great for symbolic and interactive computing with a very neat interface. You may try it for free here just to see its environment. Let's start, how Mathematica looks like. After you install it (which isn't complex, it's pretty straightforward, if you just follow the instructions), click on your desktop shortcut and it will look like the following:

how Mathematica looks like

Then, you need to click on the "Documentations" to proceed which will bring the following:

how Mathematica looks like

Now, if you like to start writing you very first code on Mathematica, then click File, and select Notebook.

how Mathematica looks like

Now, write you first code here, and execute it. Let's say, we like to plot a function which looks as follows,

plotting a function in Mathematica

To execute the program above, you need to click on Evaluation and then, select Evaluate Notebook. The variable X varies from 0 to 5.



#Mathematica #Matlab #NoteBook #MatlabvsMathematica #Blog #Blogger

A SIMULINK Model to Solve a Simple Shaft-Disk Dynamics Problem

The following figure 1 represents pretty simple model where two circular disks of inertia J1 and J2 are mentioned. Torque (T) is applied to the disk 1. The shaft has its own stiffness which is K. 

two disks are connected by a shaft
Fig.1: Showing a very simple model where two disks are connected by a shaft and torque is applied to one of them.











At first, let's find the state variables for this problem and then write the state equations.

State Variables and State Equations:

State variables



State equations


Next, we assume the following key parameters,

Initial Conditions:

Initial conditions




Parameter Values:

J1 = 100, J2 = 100, K = 100, and T = 10000. (All values are considered unitless for simplicity)

The SIMULINK model is formed by implementing the state equations. And, the above parameter values are considered. The following simulation results are for the acceleration, velocity and position of the disks which are essentially the outputs from the SIMULINK block diagram.

SIMULINK model of shaft-disk system

Fig.2: SIMULINK model of the shaft-disk system.


Showings the results of angular acceleration, velocity and displacement respectively of disk 1
Fig.3: Showings the results of angular acceleration, velocity and displacement respectively of disk 1.



angular acceleration, velocity and displacement respectively of disk 2
Fig.4: Showings the results of angular acceleration, velocity and displacement respectively of disk 2.




#Simulink #Matlab #ShaftDisk #BlockDiagram #Blog #Blogger

Theory of Energy Conversion in Wind Turbine

In wind turbine, the wind energy is converted to first mechanical energy, and then this energy is converted to electrical energy. Damping is an essential part for a generator. However, it is not the key factor for energy conversion. If there were no damping or loss, we would have 100% efficient conversion. This is not possible in real life, as we would have certain losses during the energy conversion process. These losses are represented by non-conservative forces, such as friction, viscous damping etc. which are the essential parts to be considered in the energy conversion process. I want to show the fundamentals behind energy conversion in DC motor.

DC motor converts electrical energy (input voltage) to mechanical energy (shaft rotation). This electromechanical conversion involves Faraday’s law of induction and Ampere’s law for force generated on the conductor moving in a magnetic field. In ideal situation, the torque (T) developed on the motor shaft is proportional to the input current (I) and the induced electromotive force (EMF) (V) or back EMF is proportional to the speed (W) of the motor. This can be expressed as [1];

T = K1 I ..........................................................(1)

V = K2 W .......................................................(2)

Where, K1 and K2 are the proportionality constant.
The electrical power (Pe) input to the motor is the product of the induced EMF and current.

Pe = VI = K2 W T / K1  ..................................(3)

And, the mechanical power output (Pm) is the product of the speed of the motor and torque.

Pm = T W .......................................................(4)

Now, by comparing equation (3) and (4), the following relation is obtained.

Pe = (K2 /K1) Pm  ...........................................(5)

From Ohm’s law, it is known that,

E - V = I R  ..................................................(6)

Where, E is the input voltage to the motor, and R is the resistance of the motor armature.
Moreover, we also know that torque produced at the motor shaft is equal to the product of the inertia of the load (J) and rate of change of angular velocity or angular acceleration.

T = J (dW/dt)  .............................................(7)

Now, from equations (1), (6) and (7), it is found that

J (dW/dt) = K1 I = K1 / R (E - V)  ...........................................(8)

Using equation (2) further, the following expression can be established.

dW/dt = (K1 K2 / J R) W + (K1 / J R) E  ..................................(9)

The above equation refers to the first order linear differential equation model where ‘W’ represents the state of the system and ‘E’ is the external control input. This first order equation is good enough to predict the output speed of the motor. However, in terms of measuring the position’ it is necessary to add the following equation.

W = dθ / dt  ...........................................................................(10)

Where θ is the output position of the DC motor and refers to another state of the system. Therefore, the model has one control input and two state variables (position and velocity).

By the above mathematical analysis, I want to specify, that the electrical energy is converted to mechanical energy by a gyrator which has a constant ratio K. Here, resistance of the armature is the energy loss during the conversion of electrical to mechanical energy which is also mechanical equivalence of a damper. A damper not only signifies the energy dissipation/loss from a system, but also helps to make a system stable by removing oscillations.


Reference
[1] Control system design : an introduction to state-space methods / Bernard Friedland.—Dover ed.