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

Basic Control System Development Approach for Electro-mechanical System

This article attempts to present a general approach in designing and developing a control system for any type of electro-mechanical system. The design of any robust controller for the electro-mechanical system (for example, an airship system having payload, positioner, and antenna) should involve the following steps.

a) Identification of the governing dynamics for the system

The very first step should identify the parts/elements of the system that play meaningful roles in the overall dynamics of the system. For our airship model, there are four rigid bodies in the system: airship, payload, positioner, and antenna. They are connected in a manner to form a total dynamical system. Before moving forward with the controller, we need to figure out the appropriate dynamical equations to represent the system that needs to be controlled by a controller. To do so, we may consider the following approach.

We may attach body fixed frame to each of the rigid body of the model. The differential equation for each rigid body is derived with respect to the body fixed reference frame. We may consider the ground as an absolute or inertial frame.

Once the differential equations for the rigid bodies are done, next we do the kinematic analyses that will be used for the transformation of the differential equations from the local or body fixed frame to the global or inertial frame for analyses.

We also need to identify the forces acting on the overall system for the governing equations. In practice, there are inertial forces, gravitational forces, buoyancy forces, aerodynamic forces, and so forth that act on the airship system. We may model these forces and incorporate in the differential equations.

Based on all the above assumptions, we may end up with an 8 degree-of-freedom (DOF) system where the airship has 6 DOF, and the positioner has 2 DOF. For the controller, we need to consider these degrees of freedom as a design requirement. 

b) Stability analysis of the system

We may need to linearize the system differential equations before heading towards the controller design. After the linearization, we need to check the stability of the system. Linearization may be necessary since the difficulty may arise in the controller design with complex nonlinear differential equations, which may prevent us to achieve the desired objectives. We may use the basic eigenvalue analysis to determine if the system is stable. At this stage, the linearization of the nonlinear airship model will be finalized for the design of the controller.

c) Controller design in MATLAB SIMULINK

We may choose two approaches to design a controller for the airship system. The first approach is open loop controller that is simple in design and easy to implement. The next approach may be a feedback control system that would be comprehensive and accurate in terms of system responses.

Open loop controller

In this open loop control system design step, we may implement the following SIMULINK model, which shows an actuator block and a plant or the airship system block. We are required to control the position and altitude of the motors that are connected to the antenna and positioner, respectively. 

Open loop control system


















The open loop system would be simple. The following schematic may represent the implementation method. Here, our task is to control the actuator so that it can guide the antenna. There are two actuators involved: one for 2 DOF positioner, and another for antenna. The required inputs are the desired position/altitude of the antenna that may be fed and controlled through the motor drivers. In this open loop condition, there is limited scope to optimize the system response based on the feedback errors.

Implementation of open loop controller
















Closed loop controller

We may also design a (Proportional, Integral , and Derivative) PID controller to control the position and altitude of the actuators. The following SIMULINK diagram shows a simple concept for this task:

Closed loop control system

















In this closed loop feedback control system, we may design and implement a PID controller that will control the actuators for the airship model. The input to the controller may be wind disturbance, position, altitude etc. that need to taken care of precise controlling. Since the goal of this project is to position the antenna, which is also connected to a 2 DOF positioner; we may fit a PID controller in between the actuators taking considerations all the input and output parameters for the motors.

Summary

I would recommend designing the overall airship system in SolidWorks as if it represents a real physical system. Then, the computer-aided-design (CAD) model may be translated into the MATLAB environment, keeping all the physical properties same from SolidWorks. After that, in MATLAB, we may have a lot of options to design a control system conveniently, but the added benefit here is the accurate physical CAD model that needs to be controlled eventually. In SolidWorks, we may run kinematic or dynamic motion analyses to ensure the movement of each rigid part even before the development of a controller.



#ControlSystemDesign   #MATLAB  #SIMULINK  #SolidWorks  #ElectroMechanicalSystem  #ModellingSimulation  #FeedbackControl  #MotionSimulation

Double Cantilever Beam Model to Predict Mode I Fracture for Pinned and Unpinned Laminates

This post shows the finite element (FE) model development with Abaqus CAE [1] for a double cantilever beam (DCB) to predict the mode I fracture of the z-pinned and unpinned laminates [2, 3]. This study first develops a DCB test model in Abaqus CAE with the geometric information from (ASTM D5528) [3]. The report is organized into three fundamental sections according to the project requirement.

1. Development of the DCB test model
2. Development of the Finite element model 
3. Results and discussions of the DCB model test 
3.1. Unpinned laminate simulation
3.2. Z-pinned laminate simulation
3.3. Comparison between standard test and model
        3.4. Parametric study


1. Development of the DCB test model
For the geometrical shape and size of the DCB test setup, we choose the length of the beam as 100 mm where there are two parts: 70 mm length for the adhesive material and rest 30 mm length. The adhesive layer has a depth of 1 mm, and the beam part has also a depth of 1 mm. The following schematic drawing shows the geometry for the DCB model.

Geometry for the DCB model








The following image depicts the implementation of the schematic geometry in the Abaqus CAE.

Geometry in Abaqus CAE








We have chosen materials for the cohesive and beam parts separately. The following table documents the major properties that have been implemented in the DCB test. The traction type is selected for the elastic material properties and its values are shown in the following table. For the damage criterion, the ‘Quads Damage’ option is chosen. In the ‘Damage Evolution’ module, we choose the type as ‘Energy’, then softening as ‘Linear’, then degradation as ‘Maximum’, then mixed mode behaviour as ‘BK’, then mixed mode ratio as ‘Energy’, and finally power set to 2.284. For the ‘Quads Damage’ option, we select the XFEM as ‘Normal’, where the tolerance is set to 0.05 and position is centroid.

Major properties in DCB test



















2. Development of the Finite element model 
After the geometry is created, next we move on to develop the FE model for this study. At first, we begin with meshing the geometry. The beam is meshed as ‘Quad’ and ‘Structured’ from the Abaqus meshing option, where a CPE4R beam element is selected for plain strain problem. For the cohesive part of the adhesive materials, for the mesh control, we select ‘Quad’, then we further select the ‘Sweep’, and consider the portion of highlighted orientation. We chose the cohesive element with COH2D4, where we set the viscosity to 0.00005. The following image shows the meshing for the whole structure in Abaqus.

Meshing for the structure








Next, we develop to create constraint in the DCB model so that the laminate is attached in between the cantilevers. We use tie constraint between the top surfaces of the cantilever beams to the bottom surfaces of the cohesive layers. In this way, we tie the cohesive material layer to both of the cantilever beams. The following image shows the tie constrain simulation in Abaqus.

Tie constrain simulation in Abaqus









3. Results and discussions of the DCB model test
3.1. Unpinned laminate simulation
We first show the results for the unpinned laminates in the DCB model. The following simulation results show the laminate separation for the unpinned configuration of the DCB test. 

Showing laminate separation for unpinned configuration of DCB test








The following result show the in-plane principal stress distribution in the separated DCB beam.

In-plane principal stress distribution in separated DCB beam










The following result shows the material orientations plots on the deformed shape of the beam.

Material orientations on deformed shape of beam







3.2. Z-pinned laminate simulation
Next, we test the condition for the Z-pins, where pins are inserted in the laminate layers of the adhesive material. The pin material is chosen copper, where the Young’s modulus is 130 GPa and Poisson’s ratio is 0.34. The following image shows the position of a pin inside the DCB laminate. The z-pin is inside the cohesive layer of the DCB setup.

Position of a pin inside DCB laminate





The following image shows the simulation result when a Z-pin is tied to the cohesive laminate.

Z-pin tied to cohesive laminate











3.3. Comparison between standard test and model
This section of the report discusses the comparison between the standard DCB test and the model test. We compare here force vs displacement plot for both the unpinned and pinned laminate configurations. From the following plots, we see that for both the unpinned (left) and pinned (right) laminates, there is a slight discrepancy between the standard and model tests. This is due to the facts the model test has slightly different geometry, different material properties, variations in the loads and boundary conditions, choice of the solvers, and so forth. However, the profiles are similar in both cases.

Force vs displacement plot for both unpinned and pinned laminate configurations











3.4. Parametric study
In this part, we do the parametric study for the Z-pin location and count within the laminate through the pin pull out test. We chose two pin materials, copper, and steel, in order to see the effects of them on the displacements of the adhesive layers with respect to the applied force between the standard test and the current model results. The following two tables represent the parametric study for the two materials. 

Parametric study for the two materials
















REFERENCES
[1] ABAQUS UNIFIED FEA (https://www.3ds.com/products-services/simulia/products/abaqus/)
[2] Blacklock, M, Joosten, MW, Pingkarawat, K, and Mouritz, AP, Prediction of mode I delamination resistance of z-pinned laminates using the embedded finite element technique, Composites: Part A 91 (2016) 283–291
[3] Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites, Designation: D5528 – 13




#DoubleCantileverBeam   #FiniteElement  #FractureAnalysis  #ABAQUS  #LaminateSimulation

A Fixed Free Cantilever Beam Deflection and Stress Analysis with Matlab

This tutorial is related to analyzing the deflection and stress distribution of a cantilever beam. A fixed-free cantilever beam, (Young’s Modulus 𝐸, 𝑏×ℎ cross-section, length 𝐿) is supported at its left-hand end. A Force 𝐹 is applied at a distance 𝑎 from the left-hand end. We will first calculate the shape of the deflected beam and plot the result of deflection, shear and bending moment diagram. Next, we will generate a 2D contour map of the stress distribution along the beam.

A fixed-free cantilever beam

Length of the beam, L = 1000 mm

Width of the beam, b = 50 mm

Height of the beam, h = 155 mm

Distance of the force location, a = 800 mm

Modulus of elasticity, E = 200 Gpa

Applied force, F = 25 KN

The shape of the deflected beam is calculated as [2-4],













The moment of inertia is calculated as,

I = 1/12 bh^3



Cantilever beam deflection
















Cantilever beam deflection
















The following Matlab codes are used to calculate the beam deflection and plot the results.

close all
clear 
clc
 
% Parameters from the question
E = 200*10^9;
a = 0.8;
b = 0.05;
h = 0.155;
I = (1/12)*(b*h^3);
f = 25000;
 
x = 0:0.01:0.8;
delta_1 = - (f*x.^2)./(6.*E.*I).*(3.*a - x);  % for 0 <= x <= a
delta_2 = - (f*a.^2)./(6.*E.*I).*(3.*x - a);  % for a <= x <= L
 
figure(1)
plot(x, delta_1)
xlabel('Length (m)')
ylabel('Beam Deflection (m)')
title('Cantilever Beam Deflection for 0 <= x <= a')
 
figure(2)
plot(x, delta_2)
xlabel('Length (m)')
ylabel('Beam Deflection (m)')
title('Cantilever Beam Deflection for a <= x <= L')


The shear force and bending moment diagrams for the fixed free cantilever beam is shown below:

Shear force and bending moment diagram















The Matlab codes for the calculation of shear and bending moments as well as plotting these diagrams are provided below.


MATLAB Codes

close all
clear
clc
 
L = 1; % Unit in meters
n = 2; % Number of point loads on the beam
 
for i = 1:n
fprintf('Enter the point load and distance where load acts for %d node\n',i)
Wc(i)=input('Enter the load in Newton\n');
Lc(i) = input('Enter the distance of the point load from the fixed end\n');
            
if Lc(i)>L || Lc(i)<0                
    error('Please check the Length')            
end
end
 
NL = zeros(1,n);        
NW = zeros(1,n);
       
for i=1:n
    [a,b]=max(Lc);            
    NL(i)=a;           
    NW(i)=Wc(b);          
    Lc(b)=[];           
    Wc(b)=[];        
end
 
NL(n+1) = 0;
        
% Shear force diagram
figure(1)
Ra = sum(NW);
if NL(1)==L
X1 = L;
F1 = NW(1);
elseif i==1
X1 = [L NL(1) NL(1) 0];
F1 = [0 0 Ra Ra];
else
X1 = [L NL(1)];
F1 = [0 0];
end
S=[]; X=[]; F=[];
        if n>=2
            for i=2:n+1
                x = [NL(i-1) NL(i)];
                S= [S NW(i-1)];
                f = sum(S);
                X= [X x];
                F = [F f f];
            end
        end
subplot(2,1,1)
        
plot([X1 X],[F1 F],'b')
xlim([0 L]);
hline = refline(0,0);
hline.Color = 'k';
legend('Shear Force','Reference')          
title('Shear Force Diagram of Cantilever Beam')
xlabel('Length of the Beam in Meter')
ylabel('Shear Force in Newton')
  
hold on
grid on
fprintf('Reaction Force =%d N\n',Ra)
    
% Bending moment diagram
subplot(2,1,2)
if NL(1)<L
X1=NL(1): L;
M1 = zeros(size(X1));
plot(X1,M1,'r')
hold on 
grid on
end
 
Xm=[]; 
Mm=[];
    
for i=1:n
X = NL(i+1):0.1:NL(i);
M = -(NW(1,1:i)*((NL(1,1:i))'-newX(X,i))) ;
Xm =[X Xm];
Mm =[M Mm];
end    
    
plot(Xm,Mm,'r'
xlim([0 L]);
hline = refline(0,0);
hline.Color = 'k';
legend('Bending Moment','Reference')  
    
title('Bending Moment Diagram of Cantilever Beam')
xlabel('Length of the Beam in Meter')
ylabel('Bending Moment in Newton-Meter')
hold off
grid on
 
function x = newX(X,i)                                                
[~,d] = size(X);                                                 
x = zeros(i,d);                                                
for j=1:i                                                        
x(j,1:d) = X;                                                   
end                                                             
end    

The bending stress in the beam is calculated as [2-4],





Where, δ is calculated as,








In the following, the stress distributions for the cantilever beam are shown. The Matlab codes are appended after the results.

Mesh plot for cantilever beam deflection
















Contour plot for cantilever beam deflection
















Mesh plot for cantilever beam deflection
















Contour plot for cantilever beam deflection
















MATLAB CODES

close all
clear 
clc
 
% Parameters from the question
L = 1;
E = 200*10^9;
a = 0.8;
b = 0.05;
h = 0.155;
I = (1/12)*(b*h^3);
f = 25000;
M1 = f*a;
M2 = f*(L-a);
 
x = 0:0.01:0.8;
y = 0:0.01:0.8;
 
delta_1 = - (f*x.^2)./(6.*E.*I).*(3.*a - x);  % for 0 <= x <= a
delta_2 = - (f*a.^2)./(6.*E.*I).*(3.*x - a);  % for a <= x <= L
 
sigma_1 = (M1.*delta_1)./I;
 
x = rand(100,1)*4-2;
y = rand(100,1)*4-2;
z = (M1.*((f*x.^2)./(6.*E.*I).*(3.*a - x)))./I;
z1 = (M2.*((f*a.^2)./(6.*E.*I).*(3.*x - a)))./I;
 
F = TriScatteredInterp(x,y,z);
ti = -1:.25:1;
[qx,qy] = meshgrid(ti,ti);
qz = F(qx,qy);
 
figure(1)
scatter3(x,y,z)
title('Mesh Plot for 0 <= x <= a')
hold on
mesh(qx,qy,qz)
 
figure(2)
contour(qx,qy,qz)
title('Contour Plot for 0 <= x <= a')
 
F = TriScatteredInterp(x,y,z1);
ti = -1:0.25:1;
[qx,qy] = meshgrid(ti,ti);
qz = F(qx,qy);
 
figure(3)
scatter3(x,y,z)
title('Mesh Plot for a <= x <= L')
hold on
mesh(qx,qy,qz)
 
figure(4)
contour(qx,qy,qz)
title('Contour Plot for a <= x <= L')



REFERENCES

[1] MATLAB 2020b Academic Version from MathWorks (https://www.mathworks.com/)

[2] Budynas-Nisbett, "Shigley's Mechanical Engineering Design," 8th Ed.

[3] Gere, James M., "Mechanics of Materials," 6th Ed.

[4] Lindeburg, Michael R., "Mechanical Engineering Reference Manual for the PE Exam," 13th Ed.

[5] Solid Mechanics Part I: An Introduction to Solid Mechanics by ‪Piaras Kelly‬‬



#FixedFreeBeam   #CantileverBeam   #StressDistribution   #ShearForce   #BendingMoment  #Matlab

Mohr's Circle Representation for the Stress Components of a Structure using Matlab

In this tutorial, we are going to solve a fundamental problem in mechanics of materials using the very popular Mohr's circle method. Mohr's circle is a two dimensional (2D) graphical representation of the stress components in a solid material. The following diagram shows a typical 2D element having all the stress components. Here,

𝞼x = Normal stress along x axis
𝞼y = Normal stress along y axis
𝞽xy = Shear stress on the plane
𝞡 = Inclination angle

Showing a 2D element with the stresses.















Now, let's assume the following parameters to represent the stresses by a Mohr's circle with Matlab codes.

𝞼x = Normal stress along x axis = 115 Mpa
𝞼y = Normal stress along y axis = - 50 Mpa
𝞽xy = Shear stress on the plane = 25 Mpa
𝞡 = Inclination angle = 25 Degree

With the above information at hand, we can draw a Mohr's circle to represent those stresses. The following figure shows the Mohr's circle, which is plotted by Matlab codes appended at the end of this blog post.

Showing Mohr's circle for the stress components














Showing stress versus cut plane angles














We can calculate the transformed normal and shear stresses, which  may be expressed as [2, 3],















The maximum principal stress is defined as [2, 3],











The following figure illustrates the maximum normal and shear stresses, along with other nomenclatures on a typical Mohr’s circle

Showing maximum normal and shear stresses along with other nomenclatures on the Mohr’s circle


MATLAB CODES

close all
clear
clc
 
% Parameters from the question
sigma_x = 115;
sigma_y = -50;
tau_xy = 25;
gridsize = 1000;
% Calling the function "mohrs_circle" defined in another script
[sigma_mohr,tau_mohr,sigma_1,sigma_2,tau_1,tau_2,center_circle,phi] = mohrs_circle(sigma_x,sigma_y,tau_xy,gridsize);
%%
% Plotting the figures
figure;
plot(sigma_mohr,tau_mohr);
grid on;
axis equal;
xlabel('Normal Stess, MPa');
ylabel('Shear Stress, MPa');
title('Mohr Circle for 2D Stresses');
hold on;
plot(sigma_1,0,'r*',sigma_2,0,'r*',center_circle,tau_1,'ro',center_circle,tau_2,'ro',center_circle,0,'r^');
%%
figure;
plot(phi*180/pi,sigma_mohr,'b',phi*180/pi,tau_mohr,'g');
grid on;
xlabel('Cut plane angle (deg)');
ylabel('Stress, MPA');
legend('Normal Stress','Shear Stress')
title('Mohr 2D Circle');

-----------------------------------------------------------------------------------------------------------------------------
function [sigma_mohr,tau_mohr,sigma_1,sigma_2,tau_1,tau_2,center_circle,phi]=mohrs_circle(sigma_x,sigma_y,tau_xy,gridsize)
% This function implements the calculation of the principal stresses
 
phi = linspace(0,pi,gridsize); % Defining the range of the angles for plotting
% Applying the formula or expression based on the stress theory
sigma_mohr = (sigma_x+sigma_y)/2+(sigma_x-sigma_y)/2*cos(2*phi)+tau_xy*sin(2*phi);
tau_mohr = -(sigma_x-sigma_y)/2*sin(2*phi)+tau_xy*cos(2*phi);
sigma_1 = (sigma_x+sigma_y)/2-sqrt(((sigma_x-sigma_y)/2)^2+tau_xy^2);
sigma_2 = (sigma_x+sigma_y)/2+sqrt(((sigma_x-sigma_y)/2)^2+tau_xy^2);
tau_1=sqrt(((sigma_x-sigma_y)/2)^2+tau_xy^2);
tau_2 = -tau_1;
center_circle = (sigma_x+sigma_y)/2;
%phi_p=atan(2*tau_xy/(sigma_x-sigma_y))/2;
end


REFERENCES

[1] Saul K. Fenster, Ansel C. Ugural. Advanced Strength and Applied Elasticity, Fourth Edition. Published by Pearson, 2003

[2] Solid Mechanics Part I: An Introduction to Solid Mechanics by ‪Piaras Kelly‬‬



#Matlab   #SolidMechanics   #MechanicsofMaterials #MohrsCircle