Plane Strain

clc;
clear;
close all;


[ ElementList, KTotal, KReduced, Force, ElemData ] = …
BiLinearPreProcessor( ‘Coordinates.txt’, …
‘Connectivity.txt’, ‘Supports.txt’, ‘Forces.txt’);
Supports = LoadFromFile(‘Supports.txt’ , ‘ %d ‘ , [1,Inf]);
UTemp = KReduced\Force(:, 2);
U = zeros(size(Supports,1)+size(UTemp,1),1);
SupportsTemp=sort(Supports,1);
for i=1:size(U,1)
if(isempty(SupportsTemp))
SupportsTemp=0;
end
if(SupportsTemp(1,1)==i)
U(i)=0;
SupportsTemp(1,:)=[];
else
U(i)=UTemp(1);
if(isempty(UTemp))
UTemp=0;
end
UTemp(1,:)=[];
end
end

finalForce = KTotal*U;

[ ElementList, Error ] = T4PostProcessor( ElementList, KTotal, U );

% Post-processing

% — Show max stress contour
axstress=subplot(2,2,1);
title(‘Max stress’);
DSF = 300000; %DispScaleFactor
for i=1:size(ElementList,2)
hold on;
coor = ElementList(i).Coordinates;
for j=1:4
coor(j,1)=coor(j,1)+ElementList(i).Displacements(2j-1)DSF;
coor(j,2)=coor(j,2)+ElementList(i).Displacements(2j )DSF;
end
syms x y;
strs=ElementList(i).PStress(1);
fill3(coor(:,1),coor(:,2),[1,2,3,4],…

[subs(strs,{x,y},{coor(1,:)}),…
subs(strs,{x,y},{coor(2,:)}),…
subs(strs,{x,y},{coor(3,:)}),…
subs(strs,{x,y},{coor(4,:)})]

);
colormap(axstress);
for j=1:4
disp(double(subs(ElementList(i).PStress(1),{x,y},{coor(j,:)})));
end
end

% — Show max shear stress contour
axstress=subplot(2,2,3);
title(‘Max shear stress’);
for i=1:size(ElementList,2)
hold on;
coor = ElementList(i).Coordinates;
for j=1:4
coor(j,1)=coor(j,1)+ElementList(i).Displacements(2j-1)DSF;
coor(j,2)=coor(j,2)+ElementList(i).Displacements(2j )DSF;
end
syms x y;
strs=ElementList(i).PStress(3);
fill3(coor(:,1),coor(:,2),[0,0,0,0],…

[subs(strs,{x,y},{coor(1,:)}),…
subs(strs,{x,y},{coor(2,:)}),…
subs(strs,{x,y},{coor(3,:)}),…
subs(strs,{x,y},{coor(4,:)})]

);
colormap(axstress);
for j=1:4
disp(double(subs(ElementList(i).PStress(1),{x,y},{coor(j,:)})));
end
end

% — Show displacement contour
subplot(2,2,2);
title(‘Displacement’);
for i=1:size(ElementList,2)
hold on;
coor = ElementList(i).Coordinates;
for j=1:4
coor(j,1)=coor(j,1)+ElementList(i).Displacements(2j-1)DSF;
coor(j,2)=coor(j,2)+ElementList(i).Displacements(2j )DSF;
end
syms x y;
displacement=zeros(4,1);
for j=1:4
displacement(j)=(ElementList(i).Displacements(2j-1)^2+… ElementList(i).Displacements(2j)^2)^0.5;
end
disp(displacement);
fill3(coor(:,1),coor(:,2),[1,2,3,4],displacement);
colormap jet;
end

% — Show displacement gradient
subplot(2,2,4);
title(‘Strain’);
for i=1:size(ElementList,2)
coord = ElementList(i).Coordinates;
[X,Y] = meshgrid(linspace(coord(1,1),coord(2,1),10),…
linspace(coord(2,2),coord(3,2),10));
disp([X,Y]);
dispx = double(subs(ElementList(i).Strain(1),[x,y],{X,Y}));
dispy = double(subs(ElementList(i).Strain(2),[x,y],{X,Y}));
hold on;
quiver(X,Y,dispx,dispy);
end

% — Show displacement gradient
subplot(2,2,4);
title(‘Strain’);
for i=1:size(ElementList,2)
coord = ElementList(i).Coordinates;
[X,Y] = meshgrid(linspace(coord(1,1),coord(2,1),10),…
linspace(coord(2,2),coord(3,2),10));
disp([X,Y]);
dispx = double(subs(ElementList(i).Strain(1),[x,y],{X,Y}));
dispy = double(subs(ElementList(i).Strain(2),[x,y],{X,Y}));
hold on;
quiver(X,Y,dispx,dispy);
end

Matlab output Figure
ABAQUS verification