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.

The following image depicts the implementation of the schematic geometry in the 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.

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.

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.

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. 

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

The following result shows the material orientations plots on the deformed shape of the 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.

The following image shows the simulation result when a Z-pin is tied to the 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.

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. 

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.

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

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:

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.

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

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.

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

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

App Development with the Matlab App Designer - Part 1: A Simple Mechanics Calculator

In this tutorial, I am going to walk you through the basics of the app development with the relatively new Matlab App Designer platform. Mathworks has introduced the App designer in the 2016 Matlab edition, which helps you to build your own apps with the default graphical user interface (GUI) and programming environment. I have come to know this handy tool lately, and I am very eager to share my knowledge with you, especially if you are interested in object-oriented programming. 

Let's see, first how we can access the App Designer option. Open the Matlab and click on the "APPS" from the top menus.

Matlab > APPS > Design App

Next, after clicking the "Design App" icon, the following window will appear. As you can see, there are more options to choose. But, to begin with, we should select a "Blank App".

Once, you select the "Blank App", you will see the following window. On the left hand, you have the default 'drag and drop' items for your app that you can simply select and bring it to the middle window "Design View", where the development works are carried out. On your right, you have the options to customize the default app items of the left side. For example, if you want to increase the font size of your display, you can do it from the right side menu bar.

Now, if you click on the "Code View", you will see default codes already done by Matlab for the basic interface. You can edit codes here or enter your own codes that you would like to incorporate with the simulator.

Now, let's try a very simple example, let's make the very first app by the Matlab App Designer. We will make a calculator for basic calculation in mechanics where it will calculate the force based on the two input data: mass and acceleration. As we know from Newton's second law of motion that the force is a product of the mass and acceleration, this formula is simply implemented in the development of our very first app or simulator by Matlab. Let's follow the following screenshots closely.

In the above screenshot, as you can see, from the left side menu bar, objects are dragged and placed in the middle window or the design view section. To make a calculator with two inputs, I chose a simple display for the mass input and for the acceleration, I selected a knob, which can be rotated to fix a value. To show the output that is "Force", I dragged the text objects similarly from the left side. You can also format the style of the texts, sizes, etc. from the right side menu bar.

Now, it's time to code the objects that we just placed in the design window a while ago. This is the part that we may call object-oriented programming. So, what are the objects that have been selected here - an execution button, text displays, and a knob for the acceleration parameter. If we press the execution button after providing the inputs, we should have the output displayed on the right side. This is the goal for our first app development.

We will now add codes to the "Execute" button for our goal achievement. In the following image, we see that how we can access the "Callbacks", which will eventually take us to the code environment that is shown in the next screenshot.

We need to add the formula "Force = Mass × Acceleration". To do so, as you see below, we need to copy the values of the functions under each object (e.g., "Object Mass (kg)", "Choose Acceleration", and "Force (N)") that are selected and write or organize them as per the given formula. For example, the output display is the "Force (N)" and its function value is "app.ForceNEditField.Value", which is the left side of the equation. Then, the right side is the mass times acceleration that is represented by the "app.ObjectMasskgEditField.Value" and "app.ChooseAccelerationKnob.Value" subsequently under the respective functions' values.

So, that is all! We have achieved our goal - made a very simple app that calculates the force given the mass and acceleration inputs.

#Matlab   #MatlabAppDesigner   #MechanicsCalculator  #ObjectOrientedProgramming