% Axial force and temperature by Zienkiewicz
% Problem from Beer, prob 1.2, page 13
% bar constrained at right end,  Delta T = 35
%
clear                       % variables
format short e              % E format
% format long               % double precision

temp = 35                   % temperature difference
f1 = -12e3                  % force at node 1
f2 = 2*25e3                 % force at node 2

ex = [11 6.5]*1e6           % elastic modulus, alum, magnesium
len = [3 4]*12              % element lengths
dia = [1.5 2.5]             % element diameters
alpha = [13 14 ]*1e-6       % thermal expansion
a = pi/4*dia .^2            % area
aeat = alpha .*ex .*a*temp
k12 = ex(1) * a(1)/len(1)   % EA/L
k23 = ex(2) * a(2)/len(2)

k22 = k12 + k23;
beta = k12*1e7                 % zienk factor
stiff = [ k12    -k12     0    % stiffness matrix
         -k12     k22   -k23
          0      -k23    k23 ]
df =   [ f1-aeat(1)
         f2+aeat(1)-aeat(2)
         0*beta            ] % combined d, F
zienk = stiff;
zienk(3,3)=zienk(3,3)+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(2);


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

fprintf('\n   Node Displacement  Force')
fprintf('\n            inches     Kip \n')
for i = 1:length(df)
  fprintf('    %2.0f   %10.2e %7.1f \n',i,displ(i),force(i)/1000)
end

fprintf('\n   Elem    Stress')
fprintf('\n            Ksi \n')
for i = 1:length(df)-1
  fprintf('    %2.0f   %7.2f \n',i,stress(i)/1000)
end

%        Therm6 
%  Your name, MENG 421, 05-Feb-2004

%    Node Displacement  Force
%             inches     Kip 
%      1   -4.96e-003   -12.0 
%      2    3.36e-002    50.0 
%      3    4.14e-009   -38.0 

%    Elem    Stress
%             Ksi 
%      1      6.79 
%      2     -7.74