clc;
clear all;
E=200e9;
p=0.3
d=1;
t=1;
L=5;
endload=-1000
nmeshx=32;
nmeshy=8;
%........%
Xspacing=(L/nmeshx)/2;
yspacing=(d/nmeshy)/2;
nel=nmeshx*nmeshy;
nnode=(nmeshx*2+1)*(nmeshy+1)+nmeshy*(nmeshx+1);
%.........coordinates
nodecoord=zeros(nnode,2);
c=0;
for i=1:(nmeshy*2+1)
if mod(i,2)==1
for j=1:(2*nmeshx+1)
c=c+1;
xcoord=Xspacing*(j-1);
ycoord=yspacing*(i-1);
nodecoord(c,:) = [xcoord ycoord];
end
else
for j=1:(nmeshx+1)
c=c+1;
xcoord=(2*Xspacing)*(j-1);
ycoord=yspacing*(i-1);
nodecoord(c,:)=[xcoord ycoord];
end
end
end
%..........conenectivy
conn=zeros(nel,4);
rdiff=nmeshx*2+nmeshx+2;
rdiffsmall=nmeshx*2+1;
c=0;
for j=0:nmeshy-1
for i=0:nmeshx-1
c=c+1;
c1=(2*i+1)+rdiff*j;
c2=c1+2;
c4=c1+rdiff;
c3=c4+2;
conn(c,:)=[c1;c2;c3;c4];
end
end
%..........load matrix
Load = zeros(nmeshy*2+1,3);
laod(1,:)=[rdiffsmall 0 endload/(nmeshy*2+1)];
for i=1:nmeshy*2
Load(i+1,:)=[Load(i,1)+round((nmeshx+1)*(abs(cos((i+1)*pi()/2)))+rdiffsmall*(abs(cos(i*pi()/2 )))) 0 endload/(nmeshy*2+1)];
end
nLoad=size(Load,1);
%..................support
s=zeros(nmeshy*2+1,3);
s(1,:)=[1 1 1];
for i=1:nmeshy*2 %s=[nodelid ux uy ] 1..fixed 0...free
s(i+1,:)=[Load(i,1)+ 1 1 1];
end
numsupport=size(s,1);
c=0;
for i=1:nel
nodecoordelem=zeros(4,2);
for j=1:4
nodecoordelem(j,:)=nodecoord(conn(i,j),:);
end
c=c+1;
nodecoordelemlist(c).nodecoordelem=nodecoordelem;
end
EM=E/(1-p^2).*[1 p 0;p 1 0; 0 0 (1-p)/2];
%shape functions
syms ks ta
n1 = 0.25*(1-ks)*(1-ta);
n2 = 0.25*(1+ks)*(1-ta);
n3 = 0.25*(1+ks)*(1+ta);
n4 = 0.25*(1-ks)*(1+ta);
N8=[n1;n2;n3;n4];
dnks = [-0.25*(1-ta) 0.25*1*(1-ta) 0.25*1*(1+ta) 0.25*-1*(1+ta)];
dnta = [-0.25*(1-ks) 0.25*-1*(1+ks) 0.25*1*(1+ks) 0.25*(1-ks)];
D8=[dnks
dnta];
%element stiffness matrix
for g=1:nel
jac=D8*nodecoordelemlist(g).nodecoordelem; %jacobian matrix
capgamma = inv(jac); %inverse jacobian matrix
B=sym(zeros(3,8));
for k=0:3
m=2*k+1;
n=m+1;
B(1,m) = capgamma(1,1)*dnks(k+1)+capgamma(1,2)*dnta(k+1);
B(2,n) = capgamma(2,1)*dnks(k+1)+capgamma(2,2)*dnta(k+1);
B(3,m) = B(2,n);
B(3,n) = B(1,m);
end
kelemlist(g).kelem = int(int(transpose(B)*EM*B*t*det(jac),ks,-1,1),ta,-1,1);
kelemlist(g).B = B;
end
%global stiffness matrix
kk = zeros(nnode*2,nnode*2);
for k=1:nel
kelem=kelemlist(k).kelem;
for j=1:4
for i=1:4
n=conn(k,i);
m=conn(k,j);
kk(2*n-1,2*m-1)=kk(2*n-1,2*m-1)+kelem(2*i-1,2*j-1);
kk(2*n-1,2*m)=kk(2*n-1,2*m)+kelem(2*i-1,2*j);
kk(2*n,2*m-1)=kk(2*n,2*m-1)+kelem(2*i,2*j-1);
kk(2*n,2*m)=kk(2*n,2*m)+kelem(2*i,2*j);
end
end
end
%apply the boundry condition
kkWosup=kk;
for i=1:numsupport
n=s(i,1);
if(s(i,2)==1)
kk(2*n-1,:)=0;
kk(:,2*n-1)=0;
kk(2*n-1,2*n-1)=1;
end
if(s(i,3)==1)
kk(2*n,:)=0;
kk(:,2*n)=0;
kk(2*n,2*n)=1;
end
end
%load vector
Loadvector = zeros(nnode*2,1);
for i=1:nLoad
n =Load(i,1);
Loadvector(2*n-1)=Load(i,2);
Loadvector(2*n) = Load(i,3);
end
%displacement computation
u = kk\Loadvector;
for i=1:nnode
unodelist(i).node = [u(2*i-1);u(2*i)];
end
for i=1:nel
temp = [];
for j=1:4
temp = [temp ; unodelist(conn(i,j)).node];
end
uelemlist(i).elem = temp;
end
beamtheory = (endload*L^3)/(3*E*(t*d^3/12)+(endload*L))/((E/2*(1-p))*(5/6)*d*t);
Stress=beamtheory/(2*Xspacing*yspacing*-1e-2);
%stress computation
scp = [-1 -1;1 -1;1 1;-1 1];
for i=1:nel
uelem = uelemlist(i).elem';
temp = [];
felemlist(i).felem = double(kelemlist(i).kelem*uelem');
sigma = EM*kelemlist(i).B*uelem';
temp = [];
for j=1:4
temp = [temp double(subs(subs(sigma,ks,scp(j,1)),ta,scp(j,2)))];
end
sigmaelemlist(i).sigma = temp;
end
fkk = kkWosup*u;
