i am getting two different curves for code attached below. the two figures that I have attached are for the same code, please help why am i getting this error, in actual practice I should get only one type of curve as shown in figure(i).
clc;
clear all;
% EX891.m: MATLAB program to solve a static beam deflection using
% Hermitian beam elements
% Variable descriptions
% ? = element stiffness matrix
% kk = system stiffness matrix
% ff = system force vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofe in bcdof
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=2; % number of dofs per node
nnode=(nnel-1)*nel+1; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
el=1; % elastic modulus
xi=1; % moment of inertia of cross-section
tleng=1; % length of a half of the beam
leng=tleng/nel; % element length of equal size
area=1; % cross-sectional area of the beam
rho=1; % mass density (arbitrary value for this problem because
% it is not used for the static problem)
ipt=1; % option for mass matrix (arbitrary value and not used here)
bcdof(1)=1; % first dof (deflection at left end) is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=2; % 12th dof (slope at the symmetric end) is constrained
bcval(2)=0; % whose described value is 0
displacement=zeros(sdof/2,1);
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
ff(2*nel+1)=-1; % because a half of the load is applied due to symmetry
v=1:nel;
N=nel;
p=permn(v,N);
[l r]=size(p);
mm=zeros(sdof,sdof);
force=zeros(sdof, 1);
index=zeros(nnel*ndof,1);
%%
%%
i=1:l;
kk(:,:,i)=zeros(sdof, sdof,l);
mm(:,:,i)=zeros(sdof, sdof,l);
for l=1:l
for iel=1:nel % loop for the total number of elements
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof;
%
for i=1:edof
index(i)=start+i;
end
sigma=0.55;
tleng=1;
d=tleng/nel;
i=1:nel;
j=1:nel;
ee=zeros(nel,nel);
R=zeros(nel,nel);
e(i)=tleng/2*nel+(tleng/nel)*(i-1);
for i=1:nel
for j=1:nel
ee(i,j)=e(i)-e(j);
R(i,j)=sigma^2*exp(-abs(ee(i,j)/d));
end
end
C1=chol(R,'lower');
P1=normrnd(0,1,[1,50]);
Q1=[0.18066 0.83684 -1.4424 1.3025 -0.48065 1.3233]
R1=rand(1,nel);
Z1=R1(1:nel);
alpha=C1*transpose(Z1);
E0=1;
el4=E0*(1+alpha);
el2= 1 + zeros( 1, 100 );
el=el4(1:nel);
x=p(l,iel);
c(iel)=el(iel)*xi/(leng^3);
k(:,:,iel)=c(iel)*[12 6*leng -12 6*leng; 6*leng 4*leng^2 -6*leng 2*leng^2; -12 -6*leng 12 -6*leng;6*leng 2*leng^2 -6*leng 4*leng^2];
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
kk(ii,jj,l)=kk(ii,jj,l)+k(i,j,x);
end
end
if ipt==1
%
c2(iel)=el2(iel)*rho*area*leng/420 ;
m(:,:,l)=c2(iel)*[156 22*leng 54 -13*leng;...
22*leng 4*leng^2 13*leng -3*leng^2;...
54 13*leng 156 -22*leng;...
-13*leng -3*leng^2 -22*leng 4*leng^2];
%
% lumped mass matrix
%
elseif ipt==2
m=zeros(4,4);
mass=rho*area*leng;
m=diag([mass/2 0 mass/2 0]);
%
% diagonal mass matrix
%
else
%
m=zeros(4,4);
mass=rho*area*leng;
m=mass*diag([l/2 leng^2/78 1/2 leng^2/78]);
%
end
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
mm(ii,jj,l)=mm(ii,jj,l)+m(i,j,x);
end
end
n=length(bcdof);
sdof=size(kk);
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j,l)=0;
kk(j,c,l)=0;
mm(c,j,l)=0;
mm(j,c,l)=0;
end
%
mm(c,c,l)=1;
end
end
kkk(:,:,l)=kk(3:end,3:end,l);
mmm(:,:,l)=mm(3:end,3:end,l);
[V(:,:,l) D(:,:,l)]=eig(kkk(:,:,l),mmm(:,:,l));
omega_1(:,:,l)=diag(D(:,:,l));
omega(:,:,l)=omega_1(:,:,l).^0.5;
f(:,:,l)=omega(:,:,l)/(2*pi);
v1(:,:,l)=V(:,[1 2 3],l);
end
g1=reshape(v1,2*nel,[]);
b=zeros(2,size(g1,2));
g=([b;g1]);
x=linspace(0,1);
x1=0;
x2=1/nel;
x3=1;
L1=(x.^2-x*(x2+x3)+x2*x3)/((x1-x2)*(x1-x3));
L2=(x.^2-x*(x1+x3)+x1*x3)/((x2-x1)*(x2-x3));
L3=(x.^2-x*(x1+x2)+x1*x2)/((x3-x1)*(x3-x2));
L1d=(2*x-(x2+x3))/((x1-x2)*(x1-x3));
L2d=(2*x-(x1+x3))/((x2-x1)*(x2-x3));
L3d=(2*x-(x1+x2))/((x3-x1)*(x3-x2));
L1da=(2*x1-(x2+x3))/((x1-x2)*(x1-x3));
L2db=(2*x2-(x1+x3))/((x2-x1)*(x2-x3));
L3dc=(2*x3-(x1+x2))/((x3-x1)*(x3-x2));
h1=(1-2*(x-x1)*L1da).*L1.*L1;
h2=(x-x1).*L1.*L1;
h3=(1-2*(x-x2)*L2db).*L2.*L2;
h4=(x-x2).*L2.*L2;
h5=(1-2*(x-x3)*L3dc).*L3.*L3;
h6=(x-x3).*L3.*L3;
for m=3:3:100
u=h1*g(1,m)+h2*g(2,m)+h3*g(3,m)+h4*g(4,m)+h5*g((end-1),m)+h6*g(end,m);
plot(x,u)
hold on
v=zeros(length(u),1);
plot(x,v,'r')
title('Mode shapes of the Cantilever beam ;Number of element=5')
end
figure(i) figure(ii)
