% Torsion on drilled bar, Torq4.m
clear                           % variables
format short e                  % E format

torq2 =  120                    % moment at node 1
gx = 77e9                       % shear modulus
len = [0.125 0.125]             % element lengths
dia = [20 20]*1e-3              % element diameters
diai = [ 0 16]*1e-3             % element diameters
j = pi/32*(dia .^4 - diai .^4); % polar moment
k12 = j(1)*gx/len(1)            % stiff const
k23 = j(2)*gx/len(2)
k22 = k12 + k23
beta = k12*1e7                  % zienk factor

stiff = [ k12  -k12     0       % stiffness matrix
         -k12   k22  -k23
          0    -k23   k23 ];
df =    [ 0*beta; torq2; 0*beta]  % combined d, F
zienk = stiff;
zienk(1,1)=zienk(1,1)+beta;
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       Torque4 ')
fprintf('\n  Your name, MENG 421, ')
disp(date)

fprintf('\n     Node   Rotation     Torque ')
fprintf('\n             Radian        NM \n')
for i = 1:length(df)
  fprintf('    %3.0f   %11.2e   %6.1f \n',i,rotat(i),torque(i))
end

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

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

%     Node   Rotation     Torque
%             Radian        NM
%      1     7.80e-010    -75.5
%      2     7.80e-003    120.0
%      3     4.60e-010    -44.5

%     Elem   Stress
%             MPa
%      1      48.0
%      2     -48.0