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 ...

How to Convert Graphics to Solid Body in CATIA

This tutorial demonstrates how to convert a Stereolithography Mesh (.stl) graphics object to solid editable body in CATIA V5. This is useful when we require to modify an STL object file in any CAD processing tool, such as, AutoCAD, Autodesk Inventor, SolidWorks, CATIA, and so forth. It is also important when we plan to do a Finite Element (FE) analysis of an object that is available only in the STL format. For the demonstration purpose, this tutorial uses an STL image of an anonymous patient having Scoliosis, which is available online. You can download the file from this site.

Processing an Image in CATIA

At first, open CATIA V5, then go to 'Start' and select ‘Shape’, then choose ‘Digitized Shape Editor’. See the screenshot below for the process.

Opening CATIA V5 Digitized Shape Editor




















After that, in the following window, just click ‘ok’ to create a new part file.

Creating Part in CATIA V5 Digitized Shape Editor




















We will now import the (.stl) by clicking the ‘Import’ option as highlighted below.

Importing a STL file in CATIA V5


















Next, click ‘Apply’ and then select ‘ok’ to import the (.stl) file, which looks like below:

Importing a STL file in CATIA V5




















Since, we have imported an STL file from third-party source, we need to check the quality of the image, and clean it up. There is an option ‘Mesh Cleaner’, click it, and the following window will appear.

Mesh cleaning in CATIA V5


















A pop-up small window will appear. Select the object, and then click ‘Analyze’. Click on the ‘Isolated Triangles’ option, and change the color, so that you can visualize it clearly. In this case, we see it is yellow.

The process of removing corrupt meshes from the STL file in CATIA V5



















As you notice that we are in the ‘Deletion’ mode of the ‘Mesh Cleaner’ window. So, we need to delete those unwanted or corrupted triangles to make it a valid solid body. Otherwise, the conversion either might fail or proceed, but in that case with lots of corrupt meshes that will hamper your analysis.

As we see, there are 60999 corrupt triangles, which need to be discarded. This may also remove some parts from the image, but if that part is not significant, then it may be alright to proceed depending on our requirement. We just need to make sure that we have our concerned region unaffected for the analysis. Next, click on the ‘Long Edges’ and continue the same process for removing the corrupt regions. The last option is the ‘Small Angles’, which we need to select and click ‘Apply’. 

The process of removing corrupt meshes from the STL file in CATIA V5



















After the process is done, the following window will appear. We see that there are three ‘Non-manifold Vertices’ that need to be removed as well in a similar way.
We need to continue the mesh cleaning process at least twice to make sure that we have removed all the corrupt segments. After doing so, we will have the following image without the unwanted meshes.

The process of removing corrupt meshes from the STL file in CATIA V5



















Next, we need to click on the ‘Structure’ tab and subsequently, click on the ‘Orientation’, and ‘Split in Connected Zones’ options to check the followings. If the color for the ‘Orientation’ stays same, as in this case, ‘Yellow’, then we are good to go. For the ‘Split in Connected Zones’, there should be nothing specified in the window as highlighted by the following two screenshots below.

The process of removing corrupt meshes from the STL file in CATIA V5











The process of removing corrupt meshes from the STL file in CATIA V5


















After the above process, we need to click Start > Shape > Quick Surface Reconstruction. The steps are depicted below.

Surface reconstruction process in CATIA V5


















Next, click Insert > Surface Creation > Automatic Surface. Just follow the following instructions shown by the image.

Surface reconstruction process in CATIA V5


















The following window will appear. Try to keep the parameters same as shown in the tiny popped-up window. Then, hit ‘ok’.

Surface reconstruction process in CATIA V5 and selection of key parameters


















We are almost done in creating the Solid Body. After the automatic surface creation, we need to close the surface, which is depicted below. Click Start > Mechanical Engineering > Part

Closing the surface after reconstruction in CATIA V5


















Then follow this sequence: Insert > Surface-Based Features > Close Surface.

Closing the surface after reconstruction in CATIA V5


















We are all done. Finally, we need to save it as a 'CATIA part' file and also export it to (.igs) format in case if somebody wants to perform FE analysis.

The final reconstructed 3D solid body ready for further modification


















The following image shows the closed view of the processed solid 3D Scoliosis skeleton that may be used for any analysis that needs geometry manipulation. 

Reconstructed 3D Solid Body in CATIA showing the Scoliosis Skeleton

























#Catia #ReverseEngineering #3Drendering #3Dmodel #ImageProcessing #Blog #Blogger

How to Solve a System of Partial Differential Equations in MATLAB?

In this tutorial, we are going to discuss a MATLAB solver 'pdepe' that is used to solve partial differential equations (PDEs). Let us consider the following two PDEs that may represent some physical phenomena. Sometimes, it is quite challenging to get even a numerical solution for a system of coupled nonlinear PDEs with mixed boundary conditions. This example demonstrates how we may solve a system of two PDEs simultaneously by formulating it according to the MATLAB solver format and then, plotting the results. This method is relatively easier and saves time while coding.

∂y₁/∂t = 0.375 ∂²y₁/∂x² + A(y₁ - y₂)          (1)

∂y₂/∂t = 0.299 ∂²y₂/∂x² - A(y₁ - y₂)          (2)

Where, A is a function of sinusoidal x, meaning that A = sin (x). The above two equations have two derivative terms. The time derivative, generally represents parabolic equation, and the spatial derivative defines elliptic equation. So, we have both forms of PDEs and the highest order is two in the space. There is one limitation of the 'pdepe' that you need to have the parabolic term in the PDEs in order to solve by the 'pdepe'. Now, let's assume, we have the following initial conditions:

y₁(x,0) = 1

y₂(x,0) = 0

The boundary conditions are,

∂/∂x y₁(0,t) = 0

y₂(0,t) = 0

∂/∂x y₂(1,t) = 0

y₁(1,t) = 1

The initial and boundary conditions are true if 0 ⩽ x ⩽ 1 and t ⩾ 0. Before we move forward with the coding, we need to understand first how the 'pdepe' solver accepts the PDEs in MATLAB, which form it recognizes. The general form of PDEs that the solver understands is of the following form:

The form of PDEs that the 'pdepe' solver understands





Here,

x is the independent spatial variable.

t is the independent time variable.

y is the dependent variable being differentiated with respect to x and t. It is a two-element vector where y(1) is y₁(x,t) and y(2) is y₂(x,t).

m is the symmetry constant. For cartesian coordinate, the value is 0; for cylindrical coordinate, the value is 1; and for spherical coordinate, it is 2. For our problem, it is 0 as our coordinate system is simply cartesian.

The functions c, f, and s refer to the coefficients in the above two PDE equations (1) and (2), which are required in a form that is usually expected by 'pdepe' solver.

So, in the above mentioned form, the coefficients of the PDEs may be arranged in a matrix and the equations become as,



Now, we are ready to begin the coding. In the MATLAB, we may create three functions, for example, the first function is for the equations, the second function is for the initial conditions, and the third function is for the boundary conditions. Also, we may incorporate all these functions inside another global function that is convenient as we would have everything in a single file. As the MATLAB solvers use the finite difference approach, the time integration is done with the MATLAB ‘ode15s’ solver. So, the ‘pdepe’ takes advantage of the capabilities of the stiff ‘ode15s’ solver for solving the differential-algebraic equations (DAE), which may arise when the PDEs contain elliptic equations. It is also used for handling the Jacobians with a certain sparsity pattern. Now, if there is an error during the solution process, you may try to refine the mesh size as the initial conditions are sensitive and sometimes are not consistent with the mesh size and the solver.

MATLAB Program:

function system_of_PDEs
close all
clear
clc

% Discretization of the simulation domain
x = linspace(0,10,100);
t = linspace(0,10,100);

% Application of the Matlab partial differential equation solver 'pdepe'
m = 0;
sol = pdepe(m,@pde_func,@pde_ics,@pde_bcs,x,t);

% The solution matrices
y1 = sol(:,:,1);
y2 = sol(:,:,2);

figure(1)
surf(x,t,y1)
title('y_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure(2)
surf(x,t,y2)
title('y_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

% Equations to solve
function [c,f,s] = pde_func(x,t,y,dydx) 
% Equations arranged for the 'pdepe' solver
c = [1; 1];
f = [0.375; 0.299].*dydx;
A = sin(x);
s = [A; -A];
end

% Initial Conditions
function u0 = pde_ics(x) 
u0 = [1; 0];
end

% Boundary Conditions
function [pl,ql,pr,qr] = pde_bcs(xl,ul,xr,ur,t) 
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end

end

Results:

The following results show the variations of the two dependent variables, y₁, and y₂, with respect to both space (x) and time (t). We may use the MATLAB 'Surface Plot' feature to do that.

Showing the variation of the variable y1 with respect to space and time

















Showing the variation of the variable y2 with respect to space and time


















#PDE #Matlab #FiniteDifference #HyperbolicEquation #PDEPE #Blog #Blogger

Application of the Taylor Series to Determine the Unknown Coefficients of a Polynomial

Let's say, we have a polynomial equation of the following form. 

y = ax^b + 5

Where, a, and b are the unknown coefficients of the polynomials. In this tutorial, we will discuss a method to determine the unknown coefficients using the Taylor series expansion assuming the data points are provided.

From the above polynomial equation, y is the dependent variable, and x is the independent variable. Let's assume that we know ten data points for the polynomial so that the following equations may be written based on the data points: 

y₁ = ax₁^b + 5

y2 = ax2^b + 5

y3 = ax3^b + 5

.
.
.

y10 = ax10^b + 5

This model seems to have a nonlinear dependence on the parameters a and b. Nonlinear regression can be used here to estimate these parameters. Nonlinear regression is based on determining the values of the parameters that minimize the sum of the squares of the residuals and the solution is preceded in an iterative manner. The Gauss-Newton method is one of the algorithms which minimize the sum of the squares of the residuals between nonlinear equation and data. A Taylor series expansion is used to express the nonlinear equation is an approximate linear form. After that, least square method is applied to estimate the parameters, which move towards the minimization of the residuals.

yi = f(xi; a,b,c) + ei;          i = 1, 2, 3, ... ... ... 10

Where, yi is a measured value of the dependent variable,  f(xi; a,b,c) is the nonlinear function of the independent variable, xi, with two unknown parameters a, b, and one given parameter c, which is 5 according to the given correlation. The last term ei is the random error. The model may be represented in a short form without the parameters for convenience as,

yi = f(xi) + ei                               (1)

For the two parameters case, a Taylor series expansion is written for the first two terms as,

Taylor series expansion for first two terms





Here, j is the initial guess, j+1 is the prediction, ∆a = aj+1 - aj and ∆b = bj+1 - bjEquation (2) is substituted to equation (1) which yields,

Taylor series expansion for first two terms




In matrix form, equation (3) may be written as,

Taylor series in matrix form




Where, Zj is the matrix of the partial derivatives of the function at the initial guess j.

Matrix of the partial derivatives


















The vector {D} contains the difference between the measurements and function values at 10 points.

Contains the difference between the measurements and function values











And, the vector {A} contains the changes in the parameters’ values as,

Contains the changes in the parameters






Applying least square method to equation (4) results the following normal equation,

Applying the least square method




Equation (5) is implemented to solve for {∆A} to estimate the improved values of the parameters a and b,

aj+1 = aj + ∆a

bj+1 = bj + ∆b

Thus, this procedure is continued until the solution converges.



#TaylorSeries #Polynomial #NonlinearRegression #LeastSquare #Matlab #Blog #Blogger

Discussion of the Secant Method to Solve an Equation

Let's consider the following equation:

eË£ = 2x² + 1

We may write the above expression as,

f(x) = eË£ - 2x² - 1

To find the value of x or roots of the equation, we may apply the Secant method. In this method, the derivative may be approximated by a backward finite divided difference approach. So, the function may be expressed as,
Secant method expression
This method is convenient when evaluating derivative for some function is difficult in Newton-Raphson method. Now, the above equation may be substituted in Newton-Raphson’s formula to get the following algorithm.
Secant formula
In contrast to Newton-Raphson method, two initial assumptions are made for the two unknowns’ xi-1 and xi.
For the proper initial guess for x, graphical method is applied to visualize the characteristics of the function. The MATLAB commands generate the plot of the function f(x). This function is plotted with respect to the values of x from -10 to +10. The following figure shows the function plot. From the plot, the function changes sign when x is between 2 to 3. So, the root lies in this interval. And, in the program developed for Secant method, two initial guesses for xi-1 and xi are 2 and 3 respectively.

Plot of the function f(x) to identify the interval where root of the function lies














However, more close observation reveals two roots, which are shown in the following figure. In this time, the range of x is decreased significantly for precise investigation between -2 to +2 as this interval seems unclear in the previous figure for the presence of any roots. 

Shows two more roots which lie in between 0.1 to 1, and at origin















So, the above figure shows the existence of other two roots, which lie in between 0.1 to 1 and at the origin. The origin is obvious for this case, and there is no need for any further estimation by numerical methods. But, for the interval 0.1 to 1, the accurate and precise approximation is necessary. For this, the two initial guesses are set to 0.1 and 1 for xi-1 and xi accordingly.
The following Secant formula is implemented to approximate the two roots lie in the intervals 2 to 3 and 0.1 to 1. The results are shown after the program.


Secant Formula Implementation by a MATLAB Program

% Definition of a function "secant" to solve equation by Secant Method
 
function [x,ea] = secant(X,X0,etol)
 
format long;
 
% Input: X,X0, etol=Tolerance definition for error
% Output: x, ea=Calculated error in loop
% Iterative Calculations
 
while (1)
    
    % Implementation of Secant Algorithm 
          
solution = X-((exp(X)-(2*X^2)-1)*(X0-X))/((exp(X0)-(2*X0^2)-1)-(exp(X)-(2*X^2)-1));     
    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
 
x = solution;
 
disp(x);
 
end
 
% Display of output parameters
 
disp(x);                      
 
disp(ea);
 
%disp(solution-(0.1*abs(ea)*solution));
end


Outputs:

>> secant (3, 2, 1e-5)
% Iterations

   2.597424571529533
   2.796706966294471
   2.858723109017894
   2.841817367110900
   2.842657964680533
   2.842673428005184
   2.842673428005184

>> secant(1, 0.1, 1e-5)
% Iterations

   0.308929151717719
   0.570046187158243
   1.157486960694727
   0.682985699401547
   0.723807551583747
   0.742065590882227
   0.740827136083267
   0.740850398734311
   0.740850398734311


So, the solutions for the equation are 2.842673428005184, 0.740850398734311 and 0.



#SecantMethod #RootsofEquations #Matlab #NumericalMethod #Blog #Blogger

Modelling and Simulation of an Electric Vehicle in SIMULINK

In this tutorial, I am going to present you a very simple way to design and control your very first electric vehicle in SIMULINK. You need to have MATLAB version 2019/2020 with full features. Most of the features that I will be using are from Simulink and Simscape libraries. Simscape is a very powerful tool that Matlab has integrated lately. This tutorial will be straightforward - showing you the diagrams and explaining afterwards. At first, let's begin with a simple question, how an overall electric vehicle system looks like? Below is the image of the electric vehicle system. 

Complete Simulink-Simscape block diagram of an electric vehicle system












As it appears, this is a Simulink block diagram, which contains many blocks from the library browser. Let's see what are the fundamental components of an electric vehicle system. From the diagram, there are five main blocks:

1. Vehicle Body 
2. DC Motor
3. DC Motor Driver
4. Battery
5. Controller

1. Vehicle Body

In this modelling, we can consider some basic properties of a vehicle. Below, we see the model of the vehicle body, which has four tires, a simple gear, and vehicle property block. When you click on each block, you will find the respective properties that you need to choose or enter independently. Just remember, these blocks represent the governing dynamics that Matlab has already prepared for us in the Simulink blocks. You do not need to spend hours for deriving system equations. Rather, here, we ought to be paying attention while connecting the blocks and providing appropriate parameters.

Representing vehicle body by Simscape blocks
















When you click on the vehicle body block, it will pop-up a new window as below. For this tutorial, I just chose the default parameters, but you can certainly change these for your specific problem.

Simscape vehicle body block parameters
















The following is the tire property block.

Simscape vehicle tire model block


















This is how a simple gear block appears when you double click the block.

Simscape simple gear box model block
















2. DC Motor

The purpose of the DC motor in electric vehicle is to drive the tires/wheels. It should be noted that from the Simscape library, you can choose different types of motors. The one side of the DC motor block is connected to the vehicle body blocks, and the other side is connected to the DC motor driver, which controls and rotates the motor. The following image shows the DC motor block in Simulink.

DC motor block properties


3. DC Motor Driver

To precisely control a DC motor, we are required to design a driver for that. There are numerous ways you may develop the motor controller. I used an H-Bridge block from Simscape library along with the PWM voltage controller block. In the following image, you see that there are some other blocks, such as, current and voltage sensors, electrical reference etc. 

DC motor driver system in Simulink - Simscape















4. Battery

For any electric vehicle, battery is a very important part. Still, in industries, people are struggling to find a better solution for the battery that can run an electric vehicle with more mileage before it drains fully. Here, I have chosen a simple battery model from Simscape library where you can obviously tune the parameters for your design requirement. To implement the battery and connect it to the rest of the system, a controlled current source block is added. As before, we need to have an electrical reference for the current source block, which is depicted in the subsequent figure.

Simscape battery system blocks

5. Controller

In newer versions Matlab with full features, you will find a block named 'Longitudinal Driver' that may act as a controller for the electric vehicle system. Below, we see the block parameters. I am using all the default properties for this block. You can decide the type of the controllers at first, then shift type, units, output gear signal, and so forth. There is a whole lot option to tune this block parameter for your specific goals.

Showing longitudinal driver block properties





















We begin this tutorial with the final version of the Simulink system model of an electric vehicle. There were some blocks that were masked by some images, below you can see the original version of the complete system diagram for your first virtual electric vehicle model. As you can see, in this model, most of the blocks are from the Simscape library, and you cannot connect them directly to any other Simulink blocks. But, with a converter, you can easily make the connections. In Simscape, there are two blocks - PS-Simulink and Simulink-PS converter - which are used to convert physical signals from Simscape block into Simulink compatible signal and vice versa.

Complete Simulink-Simscape block diagram of an electric vehicle system













Some Results

Everything is done. Now, if you would like to see the responses of different parts of the electric vehicle, you can do so. You need to add the oscilloscope blocks from the Simulink library to those parts. 

Here is the source signal that is fed into the system to run the vehicle.

Signal generated from the source
















This is the signal that is an output from the longitudinal driver block.

Signal from the longitudinal driver block






















This oscilloscope shows the current from a current source block.

Shows the current from a current source block.



















That's all. We have designed a virtual model for an electric car. Certainly, this is not a Tesla as the picture of the vehicle body block shows in the first image. Nevertheless, we can at least have a little sense of accomplishment that a simple model may even be designed in our laptops. 



#ElectricVehicle #Simulink #Simscape #Modelling #Simulation #DCmotor #Blog #Blogger