Stiffness Matrix for Plane Stress Problem (4-node isoparametric element with 2×2 integration points)

clear
format long

% material properties
E=200e09;
nu=0.33;

D=zeros(3,3);
D(1,1)=1;
D(1,2)=nu;
D(2,1)=nu;
D(2,2)=1;
D(3,3)=(1-nu)/2;
for i=1:3
for j=1:3
D(i,j)=E/(1-nu*nu)*D(i,j);
end
end

%coordinates of the nodes (clockwise)
coord=zeros(4,2);

coord(1,1)=0;
coord(1,2)=0;
coord(2,1)=0;
coord(2,2)=0.1;
coord(3,1)=0.1;
coord(3,2)=0.1;
coord(4,1)=0.1;
coord(4,2)=0;

%number of gauss points
%samp(i,1) -> sampling points
%samp(i,2) -> weights
ngp=2;
samp=zeros(2,2);
samp(1,1)=-0.577350269;
samp(2,1)=0.577350269;
samp(1,2)=1;
samp(2,2)=1;

Area=0;
%stiffness matrix
Kel=zeros(8,8);
for i=1:ngp
for j=1:ngp
%sampling point
xi=samp(i,1);
eta=samp(j,1);

%interpolation function evaluated at the integration point
fun(1)=(1-xi)*(1-eta)/4;
fun(2)=(1-xi)*(1+eta)/4;
fun(3)=(1+xi)*(1+eta)/4;
fun(4)=(1+xi)*(1-eta)/4;
% der stores derivatives of interpolation functions
% with respect to the natural coordinates xi, eta
der(1,1)=-(1-eta)/4;
der(1,2)=-(1+eta)/4;
der(1,3)=(1+eta)/4;
der(1,4)=(1-eta)/4;
der(2,1)=-(1-xi)/4;
der(2,2)=(1-xi)/4;
der(2,3)=(1+xi)/4;
der(2,4)=-(1+xi)/4;

% jacobian of the transformation
jac=der*coord;
% inverse of the jacobian
jac1=jac^-1;
% determinant of the jacobian
det(jac)
% derxy stores derivatives of the interpolation functions
% with respect to x and y
derxy=jac1*der;
Area=Area+det(jac)*samp(i,2)*samp(j,2);

%B matrix
for i1=1:4
B(1,(i1-1)*2+1)=derxy(1,i1);
B(2,(i1-1)*2+2)=derxy(2,i1);
B(3,(i1-1)*2+1)=derxy(2,i1);
B(3,(i1-1)*2+2)=derxy(1,i1);
end
%integration
Kel=Kel+B’*D*B*det(jac)*samp(i,2)*samp(j,2);

end
end

Area
Kel

Quadrilateral Elements (4 nodes) / Plane Elasticity

clear
format long
% 2D plane stress using 4-node quad elements
% size of the problem
Lx=1;
Ly=1;
thickness=1;

% number of elements in x and y directions
nx=10;
ny=10;

% constant element size
hx=Lx/nx;
hy=Ly/ny;

% number of nodes and number of elements
nn=(nx+1)*(ny+1);
nel=nx*ny;

% material properties
E=200e09;
nu=0.33;

% coordinates of the nodes
xcoord=zeros(nn,1);
ycoord=zeros(nn,1);

for i=1:nx+1
for j=1:ny+1
i1=i+(nx+1)*(j-1);
xcoord(i1)=(i-1)*hx;
ycoord(i1)=(j-1)*hy;
end
end

% nf stores equation numbers
% first index -> node
% second index -> local dof number
% uspecified stores the essential boundary conditions
nf=zeros(nn,2);
uspecified=zeros(nn,2);
for i=1:nx+1
nf(i,1)=1;
nf(i,2)=1;
uspecified(i,1)=0.;
uspecified(i,2)=0.;
end

% number the equations and store them in the nf array
n=0;
for i=1:nn
for j=1:2
if nf(i,j)~= 1
n=n+1;
nf(i,j)=n;
else
nf(i,j)=0;
end
end
end

% element interpolation functions
syms x y
phi(1)=(x-hx)*(y-hy)/((0-hx)*(0-hy));
phi(2)=(x-0)*(y-hy)/((hx-0)*(0-hy));
phi(3)=x*y/(hx*hy);
phi(4)=(x-hx)*y/((0-hx)*hy);

% elasticity matrix for plane stress problems
D=zeros(3,3);
D(1,1)=1;
D(1,2)=nu;
D(2,1)=nu;
D(2,2)=1;
D(3,3)=(1-nu)/2;
for i=1:3
for j=1:3
D(i,j)=E/(1-nu*nu)*D(i,j);
end
end

for i=1:4
B(1,(i-1)*2+1)=diff(phi(i),x);
B(2,(i-1)*2+2)=diff(phi(i),y);
B(3,(i-1)*2+1)=diff(phi(i),y);
B(3,(i-1)*2+2)=diff(phi(i),x);
end

% element stiffness matrix
K=int(int(thickness*B’*D*B,x,0,hx),y,0,hy);
K=double(K);

% elnode stores element node connectivity
elnode=zeros(nel,4);
m=0;
for j=1:ny
for i=1:nx
m=m+1;
elnode(m,1)=i+(nx+1)*(j-1);
elnode(m,2)=elnode(m,1)+1;
elnode(m,3)=elnode(m,2)+nx+1;
elnode(m,4)=elnode(m,3)-1;
end
end

% KG is the global stiffness matrix
KG=zeros(n,n);
% rhs is the global rhs vector
rhs=zeros(n,1);

for i=1:nel

% g is the steering vector
% g stores equation numbers
% uspecifiedel is a vector storing essential boundary conditions on
% the primary variable
jjj=0;
for j=1:4
for jj=1:2
jjj=jjj+1;
g(jjj)=nf(elnode(i,j),jj);
uspecifiedel(jjj)=uspecified(elnode(i,j),jj);
end
end

%assemble the global stiffness matrix
for i1=1:8
for j1=1:8
if g(i1)~=0 && g(j1)~=0
KG(g(i1),g(j1))=KG(g(i1),g(j1))+K(i1,j1);
end
end
end

% effect of essential boundary conditions
for i1=1:8
for j1=1:8
if g(i1) ==0 && g(j1)>0
rhs(g(j1))=rhs(g(j1))-K(i1,j1)*uspecifiedel(i1);
end
end
end
end

% external forces
j=ny+1;
for i=2:nx
i1=i+(nx+1)*(j-1);
rhs(nf(i1,2))=rhs(nf(i1,2))-200e06*thickness*hx;
end
i=1;
i1=i+(nx+1)*(j-1);
rhs(nf(i1,2))=rhs(nf(i1,2))-200e06*thickness*hx/2;

i=nx+1;
i1=i+(nx+1)*(j-1);
rhs(nf(i1,2))=rhs(nf(i1,2))-200e06*thickness*hx/2;

% solve for the primary variables
u=linsolve(KG,rhs);

% create ux and uy -> displacements in x- and y-directions
j1=0;
for i=1:nn
for j=1:2
j1=j1+1;
if nf(i,j)==0
if j==1
ux(i)=uspecified(i,j);
else
uy(i)=uspecified(i,j);
end
else
if j==1
ux(i)=u(nf(i,j));
else
uy(i)=u(nf(i,j));
end
end
end
end

fileID1 = fopen(‘ux.txt’,’w’);
fileID2 = fopen(‘uy.txt’,’w’);
for i=1:nn
fprintf(fileID1,’%12.8f %12.8f %12.8f’,[xcoord(i) ycoord(i) ux(i)] );
fprintf(fileID1,’\n’);
fprintf(fileID2,’%12.8f %12.8f %12.8f’,[xcoord(i) ycoord(i) uy(i)] );
fprintf(fileID2,’\n’);
end
fclose(fileID1);
fclose(fileID2);

fileID1 = fopen(‘sxx.txt’,’w’);
fileID2 = fopen(‘syy.txt’,’w’);
fileID3 = fopen(‘sxy.txt’,’w’);
for i=1:nel
xaverage=0.;
yaverage=0.;
for j=1:4
xaverage=xaverage+xcoord(elnode(i,j))*0.25;
yaverage=yaverage+ycoord(elnode(i,j))*0.25;
end
for j=1:4
uelement((j-1)*2+1)=ux(elnode(i,j));
uelement((j-1)*2+2)=uy(elnode(i,j));
end

stress=D*B*uelement’;
stress=subs(stress,x,hx/2);
stress=subs(stress,y,hy/2);
stress=double(stress);
fprintf(fileID1,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(1)] );
fprintf(fileID1,’\n’);
fprintf(fileID2,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(2)] );
fprintf(fileID2,’\n’);
fprintf(fileID3,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(3)] );
fprintf(fileID3,’\n’);

end
fclose(fileID1);
fclose(fileID2);
fclose(fileID3);

Triangular Elements (3 nodes) / Plane Elasticity

clear
format long

nn=5;
nel=4;

thickness=1;
E=200e09;
nu=0.33;

xcoord=zeros(4,1);
ycoord=zeros(4,1);

xcoord(1)=0.;
ycoord(1)=0.;
xcoord(2)=1.;
ycoord(2)=0.;
xcoord(3)=1.;
ycoord(3)=1.;
xcoord(4)=0.;
ycoord(4)=1.;
xcoord(5)=0.5;
ycoord(5)=0.5;

elnode=zeros(nel,3);
elnode(1,1)=1;
elnode(1,2)=2;
elnode(1,3)=5;
elnode(2,1)=2;
elnode(2,2)=3;
elnode(2,3)=5;
elnode(3,1)=5;
elnode(3,2)=3;
elnode(3,3)=4;
elnode(4,1)=1;
elnode(4,2)=5;
elnode(4,3)=4;

nf=zeros(nn,2);
uspecified=zeros(nn,2);

nf(1,1)=1;
nf(1,2)=1;
nf(2,2)=1;
uspecified(1,1)=0.;
uspecified(1,2)=0.;
uspecified(2,2)=0.;

n=0;
for i=1:nn
for j=1:2
if nf(i,j)~= 1
n=n+1;
nf(i,j)=n;
else
nf(i,j)=0;
end
end
end

D=zeros(3,3);
D(1,1)=1;
D(1,2)=nu;
D(2,1)=nu;
D(2,2)=1;
D(3,3)=(1-nu)/2;
for i=1:3
for j=1:3
D(i,j)=E/(1-nu*nu)*D(i,j);
end
end

KG=zeros(n,n);
rhs=zeros(n,1);
for i=1:nel
jjj=0;
for j=1:3
for jj=1:2
jjj=jjj+1;
g(jjj)=nf(elnode(i,j),jj);
uspecifiedel(jjj)=uspecified(elnode(i,j),jj);
end
end

for j=1:3
x(j)=xcoord(elnode(i,j));
y(j)=ycoord(elnode(i,j));
end
Area=(x(2)*y(3)-x(3)*y(2)+x(3)*y(1)-x(1)*y(3)+x(1)*y(2)-x(2)*y(1))*0.5;
B=zeros(3,6);

B(1,1)=(y(2)-y(3))/(2*Area);
B(1,3)=(y(3)-y(1))/(2*Area);
B(1,5)=(y(1)-y(2))/(2*Area);
B(2,2)=(x(3)-x(2))/(2*Area);
B(2,4)=(x(1)-x(3))/(2*Area);
B(2,6)=(x(2)-x(1))/(2*Area);
B(3,1)=B(2,2);
B(3,2)=B(1,1);
B(3,3)=B(2,4);
B(3,4)=B(1,3);
B(3,5)=B(2,6);
B(3,6)=B(1,5);

K=thickness*B’*D*B*Area;

for i1=1:6
for j1=1:6
if g(i1)~=0 && g(j1)~=0
KG(g(i1),g(j1))=KG(g(i1),g(j1))+K(i1,j1);
end
end
end

for i1=1:6
for j1=1:6
if g(i1) ==0 && g(j1)>0
rhs(g(j1))=rhs(g(j1))-K(i1,j1)*uspecifiedel(i1);
end
end
end

end

rhs(nf(3,2))=rhs(nf(3,2))-200e06*thickness/2;
rhs(nf(4,2))=rhs(nf(4,2))-200e06*thickness/2;

u=linsolve(KG,rhs);

j1=0;
for i=1:nn
for j=1:2
j1=j1+1;
if nf(i,j)==0
if j==1
ux(i)=uspecified(i,j);
else
uy(i)=uspecified(i,j);
end
else
if j==1
ux(i)=u(nf(i,j));
else
uy(i)=u(nf(i,j));
end
end
end
end

fileID1 = fopen(‘ux.txt’,’w’);
fileID2 = fopen(‘uy.txt’,’w’);
for i=1:nn
fprintf(fileID1,’%12.8f %12.8f %12.8f’,[xcoord(i) ycoord(i) ux(i)] );
fprintf(fileID1,’\n’);
fprintf(fileID2,’%12.8f %12.8f %12.8f’,[xcoord(i) ycoord(i) uy(i)] );
fprintf(fileID2,’\n’);
end
fclose(fileID1);
fclose(fileID2);

fileID1 = fopen(‘sxx.txt’,’w’);
fileID2 = fopen(‘syy.txt’,’w’);
fileID3 = fopen(‘sxy.txt’,’w’);
for i=1:nel
xaverage=0.;
yaverage=0.;
for j=1:3
xaverage=xaverage+xcoord(elnode(i,j))*0.25;
yaverage=yaverage+ycoord(elnode(i,j))*0.25;
end
for j=1:3
uelement((j-1)*2+1)=ux(elnode(i,j));
uelement((j-1)*2+2)=uy(elnode(i,j));
end
for j=1:3
x(j)=xcoord(elnode(i,j));
y(j)=ycoord(elnode(i,j));
end

Area=(x(2)*y(3)-x(3)*y(2)+x(3)*y(1)-x(1)*y(3)+x(1)*y(2)-x(2)*y(1))*0.5;
B=zeros(3,6);

B(1,1)=(y(2)-y(3))/(2*Area);
B(1,3)=(y(3)-y(1))/(2*Area);
B(1,5)=(y(1)-y(2))/(2*Area);
B(2,2)=(x(3)-x(2))/(2*Area);
B(2,4)=(x(1)-x(3))/(2*Area);
B(2,6)=(x(2)-x(1))/(2*Area);
B(3,1)=B(2,2);
B(3,2)=B(1,1);
B(3,3)=B(2,4);
B(3,4)=B(1,3);
B(3,5)=B(2,6);
B(3,6)=B(1,5);

stress=D*B*uelement’;

fprintf(fileID1,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(1)] );
fprintf(fileID1,’\n’);
fprintf(fileID2,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(2)] );
fprintf(fileID2,’\n’);
fprintf(fileID3,’%12.8f %12.8f %12.8f’,[xaverage yaverage stress(3)] );
fprintf(fileID3,’\n’);

end
fclose(fileID1);
fclose(fileID2);
fclose(fileID3);

Information Regarding the Final Project

As I have mentioned in the class before, I believe, it is more instructive to have homeworks/projects than exams in a finite element course. Therefore, instead of the final exam, I will ask you to form teams and work on a final project. Information regarding the final project is as follows:

  • Teams: Each team will consist of two students. You might choose to work on a final project alone. But I highly recommend you to form teams.
  • Topics: Each team will propose a final project topic. I recommend you to do a preliminary research on the topic and share your thoughts with me before proposing a final project topic.
  • Expected work: Each team is expected to prepare a report which includes the background information about the final project topic, finite element formulation, information regarding the finite element implementation, verification of the finite element code (with another numerical model/analytical solution/experimental data), the source code and the necessary input files, results and discussion. The deadline for the report is 11/06/2018. On the same day, each team will have a very short (10-15 minutes) presentation (schedule and location will be announced later). My expectations on the amount of work will depend on the difficulty level of the final project topic.