numberElements=2; %no. of elements in system
numberNodes=3; %no. of nodes in system
dofNode =1; %degrees of freedom per node
nodesEl =2; %nodes per element
dofEl =2; %degrees of freedom per element
elementNodes=[1 2;2 3]; %element node numbers
nodeCoordinates=[0 0;75 0;125 0]; %nodal coordinates
bcDof = [0; 125]; %boundary condition locations
A = [8;4]; %element area data
E = [25000000;25000000]; %element elasticity data
alpha = [0.0000011;0.0000011]; %element conductivity data
deltaT = [0;0]; %element temp. change data
P = [0;15000;0]; %nodal loading vector
L = [75;50]; %element length vector
xx=nodeCoordinates(:,1); %generates nodal x coordinates
yy=nodeCoordinates(:,2); %generates nodal y coordinates
sysDof = dofNode*numberNodes; %number of DOFs in the system
noConstraints = length(bcDof); %number of constrained nodes
K = zeros(sysDof,sysDof); %Initialised global stiffness
F = zeros(sysDof,1); %Initialised element force vector
%element loop
for e = 1:numberElements
AEL = A(e)*E(e)/L(e); %definition
AEalphDeltT = A(e)*E(e)*alpha(e)*deltaT(e); %definition
ke = AEL*[1 -1;-1 1 ]; %ke
fe = AEalphDeltT*[-1;1]; %element forces
gDof = elementNodes(e,:); %global DOFs for element
%system row loop
for i=1:dofEl
F(gDof(i))=F(gDof(i))+fe(i); %global element forces
%system column loop
for j=1:dofEl
fs='%s %d %s %d %s %d %s %d %s %d'; %format spec (1st input in sprintf)
sprintf(fs,'e =',e,', i =',i,', j =',j,': GDof(i) =',gDof(i),', GDof(j) =',gDof(j));
% sprintf is used to display the current values in loop
K(gDof(i),gDof(j)) = K(gDof(i),gDof(j))+ke(i,j); %global K
end
end
end
for c = 1:noConstraints
i = bcDof(c); %current boundary condition being applied
j = bcDof(c); %current boundary condition being applied
K(i,:) = 0; %replace stiffness terms in row i with zero
K(:,j) = 0; %replace stiffness terms in row j with zero
F(i) = 0; %replace force terms in row j with zero
end
PF = F+P; %add up force vectors
activeDof=setdiff(1:sysDof,bcDof); %produce a list of active DOFs
Kr = K(activeDof,activeDof); %form reduced stiffness matrix
Fr = PF(activeDof); %form reduced force vector
U = Kr\Fr; %solve for displacements
uFinal = zeros(sysDof,1); %preallocating displacement vector
for j = 1:size(activeDof,2)
uIndex = activeDof(j); %index for active DOF
uFinal(uIndex) = U(j); %allocating displacement vector according to active system DOF
end
nodalForce = K*uFinal; %nodal force
elementForce = zeros(numberElements,dofEl); %preallocating force matri
for e = 1:numberElements
sprintf('Element = %d',e); %displays current element
AEL = A(e)*E(e)/L(e); %definition
AEalphDeltT = A(e)*E(e)*alpha(e)*deltaT(e); %definition
ke = AEL*[1 -1;-1 1 ]; %element stiffness
fe = AEalphDeltT*[-1;1]; %element force
gDof = elementNodes(e,:); %call global DOFs for current element
cDofStart = gDof(1); %call 1st global DOF for current element
cDofEnd = gDof(end); %call last global DOF for current element
uCurrent = uFinal(cDofStart:cDofEnd); %displacements for current elements
elementForce(e,:) = ke*uCurrent; %element individual members force matrix
end
