% Plane stress analysis of plates
% Plane stress analysis of a thin plate under tension at its extremes
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Warning : On running this the workspace memory will be deleted. Save any
% if any data present before running the code !!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula |
% Senior Research Fellow |
% Structural Mechanics Laboratory |
% Indira Gandhi Center for Atomic Research |
% India |
% E-mail : allwayzitzme@gmail.com |
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% Variable descriptions
% k = element matrix for stiffness
% f = element vector
% stiffness = system matrix
% force = system vector
% displacement = system nodal displacement vector
% coordinates = coordinate values of each node
% nodes = nodal connectivity of each element
% index = a vector containing system dofs associated with each element
% gausspoint = matrix containing sampling points for bending term
% gaussweight = matrix containing weighting coefficients for bending term
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in 'bcdof'
% B = matrix for kinematic equation for plane stress
% D = matrix for material property for plane stress
%----------------------------------------------------------------------------
%--------------------------------------------------------------------------
% input data
%--------------------------------------------------------------------------
clear
clc
disp('Please wait Programme is under Run')
%--------------------------------------------------------------------------
% Input data for nodal coordinate values
%--------------------------------------------------------------------------
load coordinates.dat ;
%--------------------------------------------------------------------------
% Input data for nodal connectivity for each element
%--------------------------------------------------------------------------
load nodes.dat ;
%
nel = length(nodes) ; % number of elements
nnel=4; % number of nodes per element
ndof=2; % number of dofs per node (UX,UY)
nnode = length(coordinates) ; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
edof=nnel*ndof; % degrees of freedom per element
% Units are in SI system
a = 1 ; % Length of the plate (along X-axes)
b = 1 ; % Length of the plate (along Y-axes)
elementsalongX = 10 ; % Number of elements along X-axes
elementsalongY = 10 ; % Number of elements along Y-axes
%
PlotMesh(coordinates,nodes)
E = 2.1*10^11; % Youngs modulus
nu = 0.3; % Poisson's ratio
t = 0.0254; % plate thickness
rho = 7840. ; % Density of the plate
nglx = 2; ngly = 2; % 2x2 Gauss-Legendre quadrature
nglxy=nglx*ngly; % number of sampling points per element
%--------------------------------------------------------------------------
% Input data for boundary conditions
%--------------------------------------------------------------------------
% (0,0) and (1,0) are fixed
bcdof = [ 1 2 241 242] ;
bcval = zeros(1,length(bcdof)) ;
%--------------------------------------------------------------------------
% initialization of matrices and vectors
%--------------------------------------------------------------------------
force = zeros(sdof,1); % system force vector
stiffness = zeros(sdof,sdof); % system stiffness matrix
displacement = zeros(sdof,1); % system displacement vector
eldepl = zeros(edof,1) ; % element displacement vector
index = zeros(edof,1); % index vector
B = zeros(3,edof); % kinematic matrix for bending
D = zeros(3,3); % constitutive matrix for bending
stressGP = zeros(nglxy,3) ; % Matrix containing stress components
strainGP = zeros(nglx,3) ; % Matrix containing strain components
stress = zeros(nglx,3) ;
stressG = zeros(nel,3) ;
allstressGP = zeros(nnel*nel,3) ;
allstressEXP = zeros(nnel*nel,3) ;
allstressPR = zeros(nnel*nel,3) ;
%--------------------------------------------------------------------------
% force vector
%--------------------------------------------------------------------------
P = 1e5 ; % Load
rightedge = find(coordinates(:,1)==a);
rightdof = 2*rightedge-ones(length(rightedge),1);
force(rightdof) = -P*b/(elementsalongY+1) ;
leftedge = find(coordinates(:,1)==0);
leftdof = 2*leftedge-ones(length(leftedge),1) ;
force(leftdof) = P*b/(elementsalongY+1) ;
%--------------------------------------------------------------------------
% computation of element matrices and vectors and their assembly
%--------------------------------------------------------------------------
[Gausspoint,Gaussweight]=GaussQuadrature(nglx); % sampling points & weights
D = E/(1-nu^2)*[1 nu 0 ; nu 1 0; 0 0 (1-nu)/2] ; % Constituent Matrix for Plane stress
for iel=1:nel % loop for the total number of elements
for i=1:nnel
nd(i)=nodes(iel,i); % extract connected node for (iel)-th element
xx(i)=coordinates(nd(i),1); % extract x value of the node
yy(i)=coordinates(nd(i),2); % extract y value of the node
end
k = zeros(edof,edof); % initialization of stiffness matrix
%--------------------------------------------------------------------------
% numerical integration for stiffness matrix
%--------------------------------------------------------------------------
for intx=1:nglx
xi = Gausspoint(intx,1); % sampling point in x-axis
wtx = Gaussweight(intx,1); % weight in x-axis
for inty=1:ngly
eta = Gausspoint(inty,1); % sampling point in y-axis
wty = Gaussweight(inty,1) ; % weight in y-axis
[shape,dhdr,dhds] = shapefunctions(xi,eta); % compute shape functions and
% derivatives at sampling point
jacobian = Jacobian(nnel,dhdr,dhds,xx,yy); % compute Jacobian
detjacob=det(jacobian); % determinant of Jacobian
invjacob=inv(jacobian); % inverse of Jacobian matrix
[dhdx,dhdy]=shapefunctionderivatives(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
% physical coordinate
B=fekineps(nnel,dhdx,dhdy); % kinematic matrix for stiffness
%--------------------------------------------------------------------------
% compute element stiffness matrix
%--------------------------------------------------------------------------
k = k+B'*D*B*wtx*wty*detjacob;
end
end % end of numerical integration loop for bending term
index = elementdof(nd,nnel,ndof); % extract system dofs associated with element
stiffness = assemble(stiffness,k,index); % assemble element stiffness matrices
end
%--------------------------------------------------------------------------
% apply boundary conditions
%--------------------------------------------------------------------------
[stiffness,force] = constraints(stiffness,force,bcdof,bcval);
%--------------------------------------------------------------------------
% Solve the matrix equation
%--------------------------------------------------------------------------
displacement = stiffness\force ;
% Output of results
[UX,UY] = mytable(nnode,displacement,sdof);
disp('The maximum displacement UX')
max(UX)
disp('The maximum displacement UY')
max(UY)
%--------------------------------------------------------------------------
% Contour Plots of UX and UY
%--------------------------------------------------------------------------
PlotFieldonMesh(coordinates,nodes,UX)
title('Profile of UX on plate')
PlotFieldonMesh(coordinates,nodes,UY)
title('Profile of UY on plate')
