% Torsion on stepped bar, Torq1.m
clear                       % variables
format short e              % E format
torq1 = 4e3                 % moment at node 1
torq2 = -10e3               % moment at node 2

gx = [ 11 11 ]*1e6          % shear modulus
len = [4 3]*12              % element lengths
dia = [1.25  1.5]           % element diameters
j = pi/32*dia .^4;          % vector of polar moment
k12 = j(1)*gx(1)/len(1)     % stiff const
k23 = j(2)*gx(2)/len(2)
k22 = k12 + k23

beta = k12*1e7              % zienk factor
stiff = [ k12    -k12     0    % stiffness matrix
         -k12     k22   -k23
          0      -k23    k23 ]
df =    [torq1; torq2;  0*beta]; % combined d, F
zienk = stiff;
zienk(3,3)=zienk(3,3)+beta

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

for i = 1:length(df)-1        % element displacements
  eldis(i) = rotat(i+1) - rotat(i);
end
stress = gx .* dia/2 .* eldis ./ len;   % element stresses

fprintf('\n      Torque1 ')
fprintf('\n   Your name, MENG 421, ')
disp(date)
fprintf('\n     Node   Rotation     Torque ')
fprintf('\n            Radian       Kip-in \n')
for i = 1:length(df)
  fprintf('    %3.0f   %11.2e   %6.1f \n',i,rotat(i),torque(i)/1000)
end

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

%       Torque1 
%  Your name, MENG 421, 02-Mar-2004

%      Node   Rotation     Torque 
%             Radian       Kip-in 
%       1     3.33e-002      4.0 
%       2    -3.95e-002    -10.0 
%       3    -1.09e-008      6.0 

%      Elem   Stress 
%              Ksi 
%       1    -10.43 
%       2      9.05