Code
%code solves 1D elasticity
problems with up to 4 elements (5 nodes)
clc;
n=input('Press
one for one, two for two, three for three and four for four elements in the
system '); %elements in the system
clc;
if n==1
A1=input('Enter the area of the element in
mm^2 '); %area input
E1=input('Enter the Young''s modulus for
the element in MPa '); %Young's modulus input
L1=input('Enter the length of the element
in mm '); %length input
DT=input('Enter the change in temperature
in degree Celsius ');
TA1=input('Enter the coefficient of thermal
expansion for the element in per degree Celsius ');
Px1i=input('Enter the pressure acting on
the face at node i '); Px1j=input('Enter the pressure acting on the face at
node j '); %pressure input
X=input('Enter the body force acting on the
part under examination in N/mm^3 '); %body force input
P1i=input('Enter the nodal force on node i
'); P1j=input('Enter the nodal force on node j '); %nodal force input
clc;
I=[1 -1;-1 1]; kg=((A1*E1)/L1)*I; %global
stiffness matrix calculation
FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; %thermal
force vector calculation
FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1];
fp1j=A1*Px1j*FPj.'; %pressure force vector
FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; %body force
vector
FN1=[P1i P1j].'; %nodal force vector
fg=ft1+fp1i+fp1j+bf1+FN1; %global force
vector calculation
u=input('Press one for u1=0, two for u2=0
'); %boundary conditions
clc;
if u==1
fgu=zeros(1,1); fgu(1,1)=fg(2,1);
%global force vector when u1=0
kgu=zeros(1,1); kgu(1,1)=kg(2,2);
%global stiffness matrix when u1=0
U=kgu\fgu; %displacement calculation
e1=(U(1,1)-0)/L1; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; %strain due to forces
S1=ef1*E1; %stress calculation
R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction
force calculation
uU=['Displacement at node j is
',num2str(U),' mm']; disp(uU); %displacement output
rR=['Reaction force at node i is
',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element is
',num2str(ef1)]; disp(eE); %strain output
sS=['Stress in the element is
',num2str(S1),' MPa']; disp(sS); %stress output
end
if u==2
fgu=zeros(1,1); fgu(1,1)=fg(1,1);
%global force vector when u1=0
kgu=zeros(1,1); kgu(1,1)=kg(1,1);
%global stiffness matrix when u1=0
U=kgu\fgu; %displacement calculation
e1=(0-U(1,1))/L1; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; %strain due to forces
S1=ef1*E1; %stress calculation
R1=(kg(2,1)*U(1,1))-fg(2,1); %reaction
force calculation
uU=['Displacement at node i is
',num2str(U),' mm']; disp(uU); %displacement output
rR=['Reaction force at node j is
',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element is
',num2str(ef1)]; disp(eE); %strain output
sS=['Stress in the element is
',num2str(S1),' MPa']; disp(sS); %stress output
end
end
if n==2
A1=input('Enter the area of element number
one in mm^2 '); A2=input('Enter the area of element number two in mm^2 ');
%area input
E1=input('Enter the Young''s modulus for
element number one in MPa '); E2=input('Enter the Young''s modulus for element
number two in MPa '); %Young's modulus input
L1=input('Enter the length of element one
in mm '); L2=input('Enter the length of element two in mm '); %length input
DT=input('Enter the change in temperature
in degree Celsius ');
TA1=input('Enter the coefficient of thermal
expansion for element number one in per degree Celsius '); TA2=input('Enter the
coefficient of thermal expansion for element number two in per degree Celsius
'); %thermal expansion coefficient input
Px1i=input('Enter the pressure acting on
the face at node i for element number one '); Px2j=input('Enter the pressure
acting on the face at node j for element number two '); %pressure input
X=input('Enter the body force acting on the
part under examination in N/mm^3 '); %body force input
P1i=input('Enter the nodal force on node i
for element number one '); P1j=input('Enter the nodal force on node j for
element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for
element number two '); %nodal force input
clc;
I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I;
ke2=((A2*E2)/L2)*I; %element stiffness matrix calculation
kg=zeros(n+1,n+1); %initialize global
stiffness matrix
kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2);
kg(1,3)=0; kg(2,1)=ke1(2,1); kg(2,2)=ke1(1,1)+ke2(2,2); kg(2,3)=ke2(2,1);
kg(3,1)=0; kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2); %global stiffness matrix
formulation
FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.';
ft2=E2*A2*TA2*DT*FT.'; %thermal force vector
FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1];
fp2j=A2*Px2j*FPj.'; %pressure force vector
FB=[1 1]; bf1=(X*A1*L1/2)*FB.';
bf2=(X*A2*L2/2)*FB.'; %body force vector
FN1=[P1i P1j].'; FN2=[P2i P2j].'; %nodal
force vector
f1=ft1+fp1i+bf1+FN1; f2=ft2+fp2j+bf2+FN2;
%force vector calculation
fg=zeros(3,1); %global force vector
initialize
fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1);
fg(3,1)=f2(2,1); %global force vector calculation
u=input('Press one for u1=0, two for u3=0,
three for u1=u3=0 '); %boundary conditions
clc;
if u==1
fgu=zeros(2,1); fgu(1,1)=fg(2,1);
fgu(2,1)=fg(3,1); %global force vector when u1=0
kgu=zeros(2,2); kgu(1,1)=kg(2,2);
kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix
when u1=0
U=kgu\fgu; %displacement calculation
e1=(U(1,1)-0)/L1;
e2=(U(2,1)-U(1,1))/L2; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; %strain due to
forces
S1=ef1*E1; S2=ef2*E2; %stress
calculation
R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction
force calculation
uU=['Displacement at node two and three
is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement
output
rR=['Reaction force at element one,
node i is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one and two
is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output
sS=['Stress in the element one and two
is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output
end
if u==2
fgu=zeros(2,1); fgu(1,1)=fg(1,1);
fgu(2,1)=fg(2,1); %global force vector when u3=0
kgu=zeros(2,2); kgu(1,1)=kg(1,1);
kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); %global stiffness matrix
when u3=0
U=kgu\fgu; %displacement calculation
R1=(kg(3,2)*U(2,1))-fg(3,1); %reaction
force calculation
e1=(U(2,1)-U(1,1))/L1;
e2=(0-U(2,1))/L2; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; %strain due to
forces
S1=ef1*E1; S2=ef2*E2; %stress
calculation
uU=['Displacement at node one and two
is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement
output
rR=['Reaction force at element two,
node j is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one and two
is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output
sS=['Stress in the element one and two
is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output
end
if u==3
fgu=zeros(1,1); fgu(1,1)=fg(2,1);
%global force vector when u3=0
kgu=zeros(1,1); kgu(1,1)=kg(2,2);
%global stiffness matrix when u1=u3=0
U=kgu\fgu; %displacement calculation
R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction
force calculation
R3=(kg(3,2)*U(1,1))-fg(3,1); %reaction
force calculation
e1=(U(1,1)-0)/L1; e2=(0-U(1,1))/L2;
%strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; %strain due to
forces
S1=ef1*E1; S2=ef2*E2;%stress
calculation
uU=['Displacement at node two is
',num2str(U(1,1)),' mm']; disp(uU); %displacement output
rR1=['Reaction force at element one,
node i is ',num2str(R1),' N']; disp(rR1); %reaction force output
rR3=['Reaction force at element two,
node j is ',num2str(R3),' N']; disp(rR3); %reaction force output
eE=['Strain in the element one and two
is ',num2str(ef1), ' and ',num2str(ef2)]; disp(eE); %strain output
sS=['Stress in the element one and two
is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output
end
end
if n==3
A1=input('Enter the area of element number
one in mm^2 '); A2=input('Enter the area of element number two in mm^2 ');
A3=input('Enter the area of element number three in mm^2 '); %area input
E1=input('Enter the Young''s modulus for
element number one in MPa '); E2=input('Enter the Young''s modulus for element
number two in MPa '); E3=input('Enter the Young''s modulus for element number
three in MPa '); %Young's modulus input
L1=input('Enter the length of element one
in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter
the length of element three in mm '); %length input
DT=input('Enter the change in temperature
in degree Celsius ');
TA1=input('Enter the coefficient of thermal
expansion for element number one in per degree Celsius '); TA2=input('Enter the
coefficient of thermal expansion for element number two in per degree Celsius
'); TA3=input('Enter the coefficient of thermal expansion for element number
three in per degree Celsius '); %thermal expansion coefficient input
Px1i=input('Enter the pressure acting on
the face at node i for element number one '); Px3j=input('Enter the pressure
acting on the face at node j for element number three '); %pressure input
X=input('Enter the body force acting on the
part under examination in N/mm^3 '); %body force input
P1i=input('Enter the nodal force on node i
for element number one '); P1j=input('Enter the nodal force on node j for
element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for
element number two '); P3i=P2j; P3j=input('Enter the nodal force on node j for
element number three '); %nodal force input
clc;
I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I;
ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; %element stiffness matrix calculation
kg=zeros(n+1,n+1); %initialize global stiffness
matrix
kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2);
kg(2,1)=ke1(2,1); kg(2,2)=ke1(2,2)+ke2(1,1); kg(2,3)=ke2(1,2);
kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2)+ke3(1,1); kg(3,4)=ke3(1,2);
kg(4,3)=ke3(2,1); kg(4,4)=ke3(2,2); %global stiffness matrix calculation
FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.';
ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; %thermal force vector
FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1];
fp3j=A3*Px3j*FPj.'; %pressure force vector
FB=[1 1]; bf1=(X*A1*L1/2)*FB.';
bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; %body force vector
FN1=[P1i P1j].'; FN2=[P2i P2j].'; FN3=[P3i P3j].'; %nodal force vector
f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2;
f3=ft3+fp3j+bf3+FN3; %force vector calculation
fg=zeros(4,1); %global force vector
initialize
fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1);
fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1); %global force vector formulation
u=input('Press one for u1=0, two for u4=0,
three for u1=u4=0 '); %boundary conditions
clc;
if u==1
fgu=zeros(3,1); fgu(1,1)=fg(2,1);
fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); %global force vector when u1=0
kgu=zeros(3,3); kgu(1,1)=kg(2,2);
kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4);
kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); %global stiffness matrix when u1=0
U=kgu\fgu; %displacement calculation
R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction
force calculation
e1=(U(1,1)-0)/L1;
e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0;
%strain due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
%stress calculation
uU=['Displacement at node two, three
and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and
',num2str(U(3,1)),' mm']; disp(uU); %displacement output
rR=['Reaction force at element one,
node i is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one, two and
three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE);
%strain output
sS=['Stress in the element one, two and
three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa'];
disp(sS); %stress output
end
if u==2
fgu=zeros(3,1); fgu(1,1)=fg(1,1);
fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u4=0
kgu=zeros(3,3); kgu(1,1)=kg(1,1);
kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3);
kgu(3,2)=kg(3,2); kgu(3,3)=kg(3,3); %global stiffness matrix when u4=0
U=kgu\fgu; %displacement calculation
R1=(kg(4,3)*U(3,1))-fg(4,1); %reaction
force calculation
e1=(0-U(3,1))/L1;
e2=(U(3,1)-U(2,1))/L2; e3=(U(2,1)-U(1,1))/L3; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain
due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
%stress calculation
uU=['Displacement at node one, two and
three is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and ',num2str(U(3,1)),'
mm']; disp(uU); %displacement output
rR=['Reaction force at element three,
node j is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one, two and
three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE);
%strain output
sS=['Stress in the element one, two and
three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa'];
disp(sS); %stress output
end
if u==3
fgu=zeros(2,1); fgu(1,1)=fg(2,1);
fgu(2,1)=fg(3,1); %global force vector when u1=u4=0
kgu=zeros(2,2); kgu(1,1)=kg(2,2);
kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix
when u4=u1=0
U=kgu\fgu; %displacement calculation
R1=(kg(1,2)*U(1,1))-fg(1,1);
R2=(kg(4,3)*U(2,1))-fg(4,1); %reaction force calculation
e1=(U(1,1)-0)/L1;
e2=(U(2,1)-U(1,1))/L2; e3=(0-U(2,1))/L3; %strain calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0;
%strain due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
%stress calculation
uU=['Displacement at node two and three
is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement
output
rR=['Reaction force at element one,
node i and element three, node j is ',num2str(R1),' N and ',num2str(R2),' N'];
disp(rR); %reaction force output
eE=['Strain in the element one, two and
three is ',num2str(ef1),' , ',num2str(ef2),' and ',num2str(ef3)]; disp(eE);
%strain output
sS=['Stress in the element one, two and
three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa'];
disp(sS); %stress output
end
end
if n==4
A1=input('Enter the area of element number
one in mm^2 '); A2=input('Enter the area of element number two in mm^2 ');
A3=input('Enter the area of element number three in mm^2 '); A4=input('Enter
the area of element number four in mm^2 '); %area input
E1=input('Enter the Young''s modulus for
element number one in MPa '); E2=input('Enter the Young''s modulus for element
number two in MPa '); E3=input('Enter the Young''s modulus for element number
three in MPa '); E4=input('Enter the Young''s modulus for element number four
in MPa '); %Young's modulus input
L1=input('Enter the length of element one
in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter
the length of element three in mm '); L4=input('Enter the length of element
four in mm '); %length input
DT=input('Enter the change in temperature
in degree Celsius ');
TA1=input('Enter the coefficient of thermal
expansion for element number one in per degree Celsius '); TA2=input('Enter the
coefficient of thermal expansion for element number two in per degree Celsius
'); TA3=input('Enter the coefficient of thermal expansion for element number
three in per degree Celsius '); TA4=input('Enter the coefficient of thermal
expansion for element number four in per degree Celsius '); %thermal expansion
coefficient input
Px1i=input('Enter the pressure acting on
the face at node i for element number one '); Px4j=input('Enter the pressure
acting on the face at node j for element number four '); %pressure input
X=input('Enter the body force acting on the
part under examination in N/mm^3 '); %body force input
P1i=input('Enter the nodal force on node i
for element number one '); P1j=input('Enter the nodal force on node j for
element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for
element number two '); P3i=P2j; P3j=input('Enter the nodal force on node j for
element number three '); P4i=P3j; P4j=input('Enter the nodal force on node j
for element number four '); %nodal force input
clc;
I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I;
ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; ke4=((A4*E4)/L4)*I; %element stiffness
matrix calculation
kg=zeros(n+1,n+1); %initialize global stiffness
matrix
kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2);
kg(2,1)=ke1(2,1); kg(2,2)=ke1(2,2)+ke2(1,1); kg(2,3)=ke2(1,2);
kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2)+ke3(1,1); kg(3,4)=ke3(1,2);
kg(4,3)=ke3(2,1); kg(4,4)=ke3(2,2)+ke4(1,1); kg(4,5)=ke4(1,2); kg(5,4)=ke4(2,1);
kg(5,5)=ke4(2,2); %global stiffness matrix calculation
FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.';
ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; ft4=E4*A4*TA4*DT*FT.'; %thermal
force vector
FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1];
fp4j=A4*Px4j*FPj.'; %pressure force vector
FB=[1 1]; bf1=(X*A1*L1/2)*FB.';
bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; bf4=(X*A4*L4/2)*FB.'; %body force
vector
FN1=[P1i P1j].'; FN2=[P2i P2j].'; FN3=[P3i P3j].'; FN4=[P4i P4j].';%nodal force
vector
f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2;
f3=ft3+bf3+FN3; f4=ft4+fp4j+bf4+FN4; %force vector calculation
fg=zeros(5,1); %global force vector
initialize
fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1);
fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1)+f4(1,1); fg(5,1)=f4(2,1); %global force
vector formulation
u=input('Press one for u1=0, two for u5=0,
three for u1=u5=0 '); %boundary conditions0
clc;
if u==1
fgu=zeros(4,1); fgu(1,1)=fg(2,1);
fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); fgu(4,1)=fg(5,1); %global force vector when
u1=0
kgu=zeros(4,4);
kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3);
kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3);
kgu(3,3)=kg(4,4); kgu(3,4)=kg(4,5); kgu(4,3)=kg(5,4); kgu(4,4)=kg(5,5);%global
stiffness matrix when u1=0
U=kgu\fgu; %displacement calculation
R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction
force calculation
e1=(U(1,1)-0)/L1;
e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(U(4,1)-U(3,1))/L4; %strain
calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0;
ef4=e4-e0; %strain due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
S4=ef4*E4; %stress calculation
uU=['Displacement at node two, three,
four and five is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm,
',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement
output
rR=['Reaction force at element one,
node i is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one, two,
three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and
',num2str(ef4)]; disp(eE); %strain output
sS=['Stress in the element one, two,
three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa
and ',num2str(S4),' MPa']; disp(sS); %stress output
end
if u==2
fgu=zeros(4,1); fgu(1,1)=fg(1,1);
fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); fgu(4,1)=fg(4,1); %global force vector when
u5=0
kgu=zeros(4,4);
kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2);
kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3); kgu(3,2)=kg(3,2);
kgu(3,3)=kg(3,3); kgu(3,4)=kg(3,4); kgu(4,3)=kg(4,3); kgu(4,4)=kg(4,4); %global
stiffness matrix when u5=0
U=kgu\fgu; %displacement calculation
R1=(kg(5,4)*U(4,1))-fg(5,1); %reaction force calculation
e1=(U(2,1)-U(1,1))/L1;
e2=(U(3,1)-U(2,1))/L2; e3=(U(4,1)-U(3,1))/L3; e4=(0-U(4,1))/L4; %strain
calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0;
%strain due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
S4=ef4*E4; %stress calculation
uU=['Displacement at node one, two,
three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm,
',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement
output
rR=['Reaction force at element four,
node j is ',num2str(R1),' N']; disp(rR); %reaction force output
eE=['Strain in the element one, two,
three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and
',num2str(ef4)]; disp(eE); %strain output
sS=['Stress in the element one, two,
three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa
and ',num2str(S4),' MPa']; disp(sS); %stress output
end
if u==3
fgu=zeros(3,1); fgu(1,1)=fg(1,1);
fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u1=u5=0
kgu=zeros(3,3);
kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3);
kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3);
kgu(3,3)=kg(4,4); %global stiffness matrix when u5=0
U=kgu\fgu; %displacement calculation
R1=(kg(5,4)*U(3,1))-fg(5,1);
R2=(kg(1,2)*U(1,1))-fg(5,1); %reaction force calculation
e1=(U(1,1)-0)/L1;
e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(0-U(3,1))/L4; %strain
calculation
e0=TA1*DT; %thermal strain
ef1=e1-e0; ef2=e2-e0; ef3=e3-e0;
ef4=e4-e0; %strain due to forces
S1=ef1*E1; S2=ef2*E2; S3=ef3*E3;
S4=ef4*E4; %stress calculation
uU=['Displacement at node two, three
and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),'
mm']; disp(uU); %displacement output
rR=['Reaction force at element one,
node i and element four, node j is ',num2str(R1),' N and ',num2str(R2),' N'];
disp(rR); %reaction force output
eE=['Strain in the element one, two,
three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and
',num2str(ef4)]; disp(eE); %strain output
sS=['Stress in the element one, two,
three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa
and ',num2str(S4),' MPa']; disp(sS); %stress output
end
end
Sample calculations for problem with one element
With first boundary condition (u1=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 1
Enter the area of the element in mm^2 20
Enter the Young's modulus for the element in MPa 200000
Enter the length of the element in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for the element
in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i 0
Enter the pressure acting on the face at node j 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i 0
Enter the nodal force on node j 10000
Press one for u1=0, two for u2=0 1
Result
Displacement at node j is 0.26 mm
Reaction force at node i is -10000 N
Strain in the element is 0.0025
Stress in the element is 500 MPa
With second boundary condition (u2=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 1
Enter the area of the element in mm^2 20
Enter the Young's modulus for the element in MPa 200000
Enter the length of the element in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for the element
in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i 0
Enter the pressure acting on the face at node j 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i -10000
Enter the nodal force on node j for element number one 0
Press one for u1=0, two for u2=0 2
Result
Displacement at node i is -0.26 mm
Reaction force at node j is 10000 N
Strain in the element is 0.0025
Stress in the element is 500 MPa
Sample calculations for problem
with two elements
Problem 5.1a, Fagan pg. 94
With first boundary condition (u1=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 2
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 10
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number two 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 10000
Press one for u1=0, two for u3=0, three for u1=u3=0 1
Result
Displacement at node two and three is 0.26 mm and 0.77 mm
Reaction force at element one, node i is -10000 N
Strain in the element one and two is 0.0025 and 0.005
Stress in the element one and two is 500 MPa and 1000 MPa
With second boundary condition (u3=0)
Problem 5.1a, Fagan pg. 94 mirrored
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 2
Enter the area of element number one in mm^2 10
Enter the area of element number two in mm^2 20
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number two 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one
-10000
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 0
Press one for u1=0, two for u3=0, three for u1=u3=0 2
Result
Displacement at node one and two is -0.77 mm and -0.26 mm
Reaction force at element two, node j is 10000 N
Strain in the element one and two is 0.005 and 0.0025
Stress in the element one and two is 1000 MPa and 500 MPa
With third boundary condition (u3=u1=0)
Problem 5.1b, Fagan pg. 94
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 2
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 10
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number two 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 0
Press one for u1=0, two for u3=0, three for u1=u3=0 3
Result
Displacement at node two is 0.0033333 mm
Reaction force at element one, node i is 266.6667 N
Reaction force at element two, node j is -266.6667 N
Strain in the element one and two is -6.6667e-05 and
-0.00013333
Stress in the element one and two is -13.3333 MPa and
-26.6667 MPa
Sample calculations for problem
with three elements
Problem 5.1f, Fagan pg. 95
With first boundary condition (u1=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 3
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 20
Enter the area of element number three in mm^2 10
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the change in temperature in degree Celsius 0
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number three in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number three 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 20000
Enter the nodal force on node j for element number two 0
Enter the nodal force on node j for element number three
10000
Press one for u1=0, two for u4=0, three for u1=u4=0 1
Result
Displacement at node two, three and four is 1.25 mm, 1.5 mm
and 2 mm
Reaction force at element one, node i is -50000 N
Strain in the element one, two and three is 0.0125, 0.0025
and 0.005
Stress in the element one, two and three is 2500 MPa, 500
MPa and 1000 MPa
With second boundary condition (u4=0)
Problem 5.1f, Fagan pg. 95 mirrored
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 3
Enter the area of element number one in mm^2 10
Enter the area of element number two in mm^2 20
Enter the area of element number three in mm^2 20
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the change in temperature in degree Celsius 0
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number three in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number three 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one
-10000
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two
-20000
Enter the nodal force on node j for element number three 0
Press one for u1=0, two for u4=0, three for u1=u4=0 2
Result
Displacement at node one, two and three is -2 mm, -1.5 mm
and -1.25 mm
Reaction force at element three, node j is 50000 N
Strain in the element one, two and three is 0.0125, 0.0025
and 0.005
Stress in the element one, two and three is 2500 MPa, 500
MPa and 1000 MPa
With third boundary condition (u4=u1=0)
Problem 5.1f, Fagan pg. 95 with temperature load only, no
nodal forces
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 3
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 20
Enter the area of element number three in mm^2 10
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number three in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number three 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 0
Enter the nodal force on node j for element number three 0
Press one for u1=0, two for u4=0, three for u1=u4=0 3
Result
Displacement at node two and three is 0.0025 mm and 0.005 mm
Reaction force at element one, node i and element three,
node j is 300 N and -300 N
Strain in the element one, two and three is -7.5e-05 ,
-7.5e-05 and -0.00015
Stress in the element one, two and three is -15 MPa, -15 MPa
and -30 MPa
Sample calculations for problem
with four elements
With first boundary condition (u1=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 4
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 15
Enter the area of element number three in mm^2 10
Enter the area of element number four in mm^2 5
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the Young's modulus for element number four in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the length of element four in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number three in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number four in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number four 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 400
Enter the nodal force on node j for element number one 200
Enter the nodal force on node j for element number two 400
Enter the nodal force on node j for element number three 200
Enter the nodal force on node j for element number four 400
Press one for u1=0, two for u5=0, three for u1=u5=0 1
Results
Displacement at node two, three, four and five is 0.06 mm,
0.12333 mm, 0.17333 mm and 0.22333 mm
Reaction force at element one, node i is -2400 N
Strain in the element one, two, three and four is 0.0005,
0.00053333, 0.0004 and 0.0004
Stress in the element one, two, three and four is 100 MPa,
106.6667 MPa 80 MPa and 80 MPa
With second boundary condition (u5=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 4
Enter the area of element number one in mm^2 5
Enter the area of element number two in mm^2 10
Enter the area of element number three in mm^2 15
Enter the area of element number four in mm^2 20
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the Young's modulus for element number four in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the length of element four in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number three in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number four in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number four 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one -400
Enter the nodal force on node j for element number one -200
Enter the nodal force on node j for element number two -400
Enter the nodal force on node j for element number three
-200
Enter the nodal force on node j for element number four -400
Press one for u1=0, two for u5=0, three for u1=u5=0 2
Results
Displacement at node one, two, three and four is -0.22333
mm, -0.17333 mm, -0.12333 mm and -0.06 mm
Reaction force at element four, node j is 2400 N
Strain in the element one, two, three and four is 0.0004,
0.0004, 0.00053333 and 0.0005
Stress in the element one, two, three and four is 80 MPa, 80
MPa 106.6667 MPa and 100 MPa
With third boundary condition (u5=u1=0)
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 4
Enter the area of element number one in mm^2 20
Enter the area of element number two in mm^2 15
Enter the area of element number three in mm^2 10
Enter the area of element number four in mm^2 5
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the Young's modulus for element number three in MPa
200000
Enter the Young's modulus for element number four in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the length of element three in mm 100
Enter the length of element four in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element number
three in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number four in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number four 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 0
Enter the nodal force on node j for element number three 0
Enter the nodal force on node j for element number four 0
Press one for u1=0, two for u5=0, three for u1=u5=0 3
Results
Displacement at node two, three and four is -0.0058 mm,
-0.0002 mm, 0.0032 mm
Reaction force at element one, node i and element four, node
j is -132 N and 132 N
Strain in the element one, two, three and four is -0.000158,
-4.4e-05, -6.6e-05 and -0.000132
Stress in the element one, two, three and four is -31.6 MPa,
-8.8 MPa -13.2 MPa and -26.4 MPa
Sample Calculations for problems involving
variable cross-sectional area
Problem 5.1d, Fagan pg. 95
The command window at various steps
Press one for one, two for two, three for three and four for
four elements in the system 2
Enter the area of element number one in mm^2 25
Enter the area of element number two in mm^2 15
Enter the Young's modulus for element number one in MPa
200000
Enter the Young's modulus for element number two in MPa
200000
Enter the length of element one in mm 100
Enter the length of element two in mm 100
Enter the change in temperature in degree Celsius 50
Enter the coefficient of thermal expansion for element
number one in per degree Celsius 2/1000000
Enter the coefficient of thermal expansion for element
number two in per degree Celsius 2/1000000
Enter the pressure acting on the face at node i for element
number one 0
Enter the pressure acting on the face at node j for element
number two 0
Enter the body force acting on the part under examination in
N/mm^3 0
Enter the nodal force on node i for element number one 0
Enter the nodal force on node j for element number one 0
Enter the nodal force on node j for element number two 10000
Press one for u1=0, two for u3=0, three for u1=u3=0 1
Results
Displacement at node two and three is 0.21 mm and 0.55333 mm
Reaction force at element one, node i is -10000 N
Strain in the element one and two is 0.002 and 0.0033333
Stress in the element one and two is 400 MPa and
666.6667 MPa