Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Problem 1
syms s
Kmid=[1 -1;-1 1];
E=80*10^9; u=[];
ERR=[];
for NELEM=1:40
X=linspace(0,.2,NELEM+1);
Le=X(NELEM+1)-X(NELEM);
GSTIFF=zeros(NELEM+1,NELEM+1);
FE=zeros(NELEM+1,1);
for i=1:NELEM KSTIFF=E*(1/(2*Le))*Kmid*int(0.001*(2-
4.5*(Le*s+X(i)+X(i+1))),-1,1);
RN=[i,i+1];
GSTIFF(RN,RN)=GSTIFF(RN,RN)+KSTIFF;
f=10^3*5*(Le/4)*int((Le*s+X(i)+X(i+1))*[1-s;1+s],-1,1);
FE(RN,1)=FE(RN,1)+f; end %now i have Gstiff and FE
%now implying boundary conditions
FE(NELEM+1,1)=FE(NELEM+1,1)+5*10^3;
%imposing PENALTY to satisfy BC
penalty=10^12; %u(0)=0
GSTIFF(1,1)=GSTIFF(1,1)+penalty;
u=linsolve(GSTIFF,FE); Uact=1.628952867*10^(-5);%@last node exact
solution analytical computation ERR(NELEM)=Uact-u(NELEM+1); %@last node
end disp('displacement
is: ') disp(u) plot(ERR)
displacement is: 1.0e-
04 *
0.0001
0.0017
0.0034
0.0051
0.0069
0.0087
0.0105
0.0124
0.0144
0.0164
0.0184
0.0206
0.0228
0.0250
0.0273
0.0297
Scanned by CamScanner
0.0322
0.0348
0.0374
0.0402
0.0431
0.0461
0.0492
0.0524
0.0558
0.0594
0.0631
0.0671
0.0713
0.0757
0.0805
0.0855
0.0910
0.0969
0.1033
0.1103
0.1182
0.1270
0.1370
0.1487
0.1628
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
clear
all clc
%defination of data
E=200*10^9; %elastic constant
A=pi/4*(2*10^(-3))^2; %area of truss element
%nodal coordinates
NodeCoord=[0 0;24 10;36 15;48 20;60 15;72 10;96 0;84 15;72 32;60 37;48 42;
36 37;24 32;12 15];
%connectivity
Con=[1 2;2 3;3 4;4 5;5 6;6 7;7 8;8 9;9 10;10 11;11 12;12 13;13 14;14 1;
2 14;2 13;2 12;3 12;4 12;4 11;4 10;5 10;6 10;6 9;6 8];
FE=zeros(28,1);
Gstiff=zeros(28,28);
for i=1:length(Con)
element=Con(i,:);
%finding element dimensions
y1=NodeCoord(element(1),2);
y2=NodeCoord(element(2),2);
x1=NodeCoord(element(1),1);
x2=NodeCoord(element(2),1);
dy=(y2-y1); dx=(x2-x1);
theta=vpa(atan(dy/dx)); Le=vpa(sqrt(dy^2+dx^2)); c=vpa(cos(theta));
s=vpa(sin(theta)); Le=0.0254*Le; Kelem=E*(A/Le)*[c*c s*c -c*c -s*c;s*c
s*s -s*c -s*s;-c*c -s*c c*c s*c;-s*c -s*s %forming Global stiffness matrix
%finding row number
RN=[2*element(1)-1,2*element(1),2*element(2)-1,2*element(2)];
Gstiff(RN,RN)=Gstiff(RN,RN)+Kelem;
end
%providing boundary conditions
FE(8)=vpa(200*0.453592*9.81)+FE(8);
penalty=[10^30 10^30;10^30 10^30];
Scanned by CamScanner
%node 1 and 7 is having zero
displacement
Gstiff([1 2],[1 2])=Gstiff([1 2],[1 2])+penalty;
Gstiff([13 14],[13 14])=Gstiff([13 14],[13 14])+penalty;
%finding displacement of each node
u=vpa(linsolve(Gstiff,FE))
Warning: Matrix is close to singular or badly scaled. Results may be
inaccu RCOND = 2.016641e-26. u =
0.0000000000010674864189433048803027510484542
-0.0000000000010674864189433042744575758387172
0.0000062503924293070272855056719596423
0.0010473542522755055136918223013254
-0.000016147809394080716086355281602138
0.001505059973595197081452923271172
-0.0000000000049239017804896344873029677396343
0.0018702552798103289077857880329248
0.000016147800427363498859982191002693
0.0015050599757098037483060704033733
-0.0000062504005149379740808733736667779
0.0010473542565047190642385510628287
-0.0000000000073909439508234375667931588059372
0.0000000000073909439508234375667931588059372
0.00021496499258478158982631212037262
0.00054824420290592180537297961606669
0.00010779455993022140226477584290521
0.00085625232391710261355133315674948
0.0002456023738416395407399561712225
0.0015050599757098039651465049004742
-0.0000000000068622913984415729801553403843388
0.0015517520588309672169985420353555
-0.00024560238668513626044895281630431
0.0015050599735951972982933577682729
-0.00010779457189263202433748833897198
0.00085625231968788884616416989814525
-0.0002149650015514991492539959017094
0.00054824419656210039535071398830723
Scanned by CamScanner
clear all
E=200*10^9; I=1; %assuming M.O.I. of
beam to be 1 m4
syms s; NELEM=200; %total no
of elements
Gstiff=zeros(2*(NELEM+1),2*(NELEM+
1)); FE=zeros(2*(NELEM+1),1); for
i=1:NELEM
x=linspace(0,5,NELEM+1); Le=x(NELEM)-x(NELEM-1); Ni=[1-
3*s^2+2*s^3;s*Le-2*s^2*Le+s^3*Le;3*s^2-2*s^3;-Le*s^2+Le*s^3]; %shape function
if i<=NELEM/5
q1=(1-
(Le*s+x(i))*3*10^3+10^3);
Rq=int(q1*Ni*Le,0,1); end
if i>NELEM/5 &&
i<=4*NELEM/5 q2=10^3;
Rq=int(q2*Ni*Le,0,1);
end
if i>4*NELEM/5
q3=(Le*s+x(i)-
4)*3*10^3;
Rq=int(q3*Ni*Le,0,1);
end
Kstiff=E*I/Le^3*[12 6*Le -12 6*Le;6*Le 4*Le^2 -6*Le 2*Le^2;
-12 -6*Le 12 -6*Le;6*Le 2*Le^2 -6*Le 4*Le^2]; RN=[2*i-
1,2*i,2*(i+1)-1,2*(i+1)];
FE(RN,1)=FE(RN,1)+Rq;
Gstiff(RN,RN)=Gstiff(RN,RN)+Kstiff; end
%applying BC FE(2*(NELEM/2)+1,1)=FE(2*(NELEM/2)+1,1)-
100*10^3;
penalty=10^25;
Gstiff(1,1)=Gstiff(1,1)+penalty;
Gstiff(2*(NELEM+1)-1,2*(NELEM+1)-1)=Gstiff(2*(NELEM+1)-1,2*(NELEM+1)-
1)+penalty;
Gstiff(2*(4*NELEM/5+1)-1,2*(4*NELEM/5+1)-1)=Gstiff(2*(4*NELEM/5+1)-
1,2*(4*NELEM/5+1
Gstiff(2*(NELEM/5+1)-1,2*(NELEM/5+1)-1)=Gstiff(2*(NELEM/5+1)-1,2*(NELEM/5+1)-
1)+pen
Gstiff(2*(NELEM/2+1),2*(NELEM/2+1))=Gstiff(2*(NELEM/2+1),2*(NELEM/2+1))+penalty
;
v=vpa(linsolve(Gstiff,FE));
%finding value of deflection
only u=[]; p=1;
Scanned by CamScanner
for
i=1:2*(NELEM+1)
if mod(i,2)==1
u(p,1)=v(i,1);
p=p+1; end
end plot(u)
Scanned by CamScanner
%finding moment values Mom=[];
du=[]; for i=1:NELEM
du(i)=(u(i+1)-u(i))/Le; end
for i=1:(NELEM-1)
Mom(i)=E*I*(du(i+1)-du(i))/Le;
end figure plot(Mom)
ylabel('Moment diagram')
%finding shear force ShrFrs=[];
for i=1:(NELEM-2)
ShrFrs(i)=(Mom(i+1)-Mom(i))/Le;
end figure plot(ShrFrs)
ylabel('shear force diagram')
Scanned by CamScanner
Scanned by CamScanner
Scanned by CamScanner
