% stepped bar, thermex04.m 
% element 1, 2 steel, elements 3 aluminum
% fixed on left end, 25F temperature change
%
clear                       % variables
format short e              % E format

temp = 25                   % temperature difference
f2 = -2*4e3                 % force at node 2
f3 = -2*2e3                 % force at node 2
f4 =  10e3                  % force at node 1

ex = [ 29 29 11 ]*1e6       % elastic modulus, steel, alum
len = [2 3 4]*12            % element lengths
dia = [1.5 1.5 1.2]         % element diameters
alpha = [6.5 6.5 13.1]*1e-6 % thermal expansion
area = pi/4*dia .^2
aeat = alpha .*ex .*area*temp
k12 = ex(1) * area(1)/len(1) % EA/L
k23 = ex(2) * area(2)/len(2)
k34 = ex(3) * area(3)/len(3)
k22 = k12 + k23;
k33 = k23 + k34;
beta = k12*1e7                 % zienk factor
stiff = [ k12    -k12     0     0  % stiffness matrix
         -k12     k22   -k23    0
          0      -k23    k33  -k34
          0        0    -k34   k34 ]
df =    [ 0*beta
          f2+aeat(1)-aeat(2)
          f3+aeat(2)-aeat(3)
          f4+aeat(3)           ]   % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;

displ = inv(zienk) * df;     % displacements
force = stiff * displ;       % force vector

for i = 1:length(df)-1       % element displacements
  eldis(i) = displ(i+1) - displ(i);
end

% find element stresses but remove a*e*t
stress = ex .* (eldis ./ len - alpha*temp);
force(1) = force(1) + aeat(1);   % remove aeat
force(2) = force(2) + aeat(2) - aeat(1);
force(3) = force(3) + aeat(3) - aeat(2);
force(4) = force(4) - aeat(3);

% print the results
fprintf('\n         AxialAxE04 ')
fprintf('\n  Your name, MENG 421, ')
disp(date)

fprintf('\n   Node Displacement    Force')
fprintf('\n            Inches       Lb \n')
for i = 1:length(df)
  fprintf('    %2.0f   %11.3e %7.0f \n',i,displ(i),force(i))
end

fprintf('\n   Elem    Stress')
fprintf('\n             Psi \n')
for i = 1:length(df)-1
  fprintf('    %2.0f   %7.0f \n',i,stress(i))
end

% End of Matlab code

%          AxialAxE04 
%   Your name, MENG 421, 05-Mar-2004

%    Node Displacement    Force
%             Inches       Lb 
%      1    2.963e-010    2000 
%      2    2.963e-003   -8000 
%      3    1.303e-002   -4000 
%      4    6.733e-002   10000 

%    Elem    Stress
%              Psi 
%      1     -1132 
%      2      3395 
%      3      8842