% Axial force and temperature by Zienkiewicz, therm7.m
% Problem from Hibler 4-37, 38, page146
% bar constrained at each end,  Delta T = 90
%
clear                       % variables
format short e              % E format
ex = [70 200 70]*1e9      % elastic modulus, steel, alum
len = [0.25 0.5 0.25]       % element lengths
dia = [50 20 50]*1e-3       % element diameters
alpha = [23 11.7 23]*1e-6     % thermal expansion
temp = 65                   % temperature difference
f1 = -2*75e3                % force at node 2
f2 = 2*100e3                % force at node 2

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
         f1+aeat(1)-aeat(2)
         f2+aeat(2)-aeat(3)
		 0*beta ] % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
zienk(4,4)=zienk(4,4)+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       Therm7 ')
fprintf('\n      Your name, ES 421, ')
disp(date)

fprintf('\n   Node  Displacement   Force')
fprintf('\n            meters       kN \n')
for i = 1:length(df)
  fprintf('    %2.0f   %11.3e %6.0f \n',i,displ(i),force(i)/1000)
end

fprintf('\n   Elem    Stress')
fprintf('\n            MPa \n')
for i = 1:length(df)-1
  fprintf('    %2.0f   %7.1f \n',i,stress(i)/1e6)
end

%      ThermE01
%      Your name, ES 421, 28-Feb-2001

%   Node  Displacement   Force
%            meters       kN
%     1    2.387e-012    192
%     2    2.387e-005   -150
%     3    6.707e-005    200
%     4    6.707e-012   -242

%   Elem    Stress
%            MPa
%     1     -98.0
%     2    -134.8
%     3    -123.4