transmission_tower.m

LU decomposition routines in matlab

% transmission_tower.m
% set up and solve the equations for a symmetric transmission tower
%
% Kaveh Ghayour, Matthias Heinkenschloss
% Jan 26, 2001



% set number of bar (m) and number of node (n)
m = 40;
n = 19;

% bar is an m times 2 array. The i-th row contains
% information for bar i in the following form:
%  [ lower-left-node   upper-right-node ]

bar ...
= [ 1   2;
    1   4;
    2   3;
    3   4;
    2   4;
    1   3;
    3   5;
    3   6;
    4   5;
    4   6;
    5   6;
    5   8;
    6   7;
    7   8;
    5   7;
    6   8;
    8   9;
    7   10;
    7   9;
    8   10;
    9   10;
    10  12;
    9   13;
    9   12;
    10  13;
    10  14;
    13  14;
    11  12;
    9   11;
    12  13;
    12  15;
    13  15;
    16  9;
    16  11;
    17  16;
    17  11;
    10  18;
    18  14;
    18  19;
    19  14];
% node is an n times 2 array. The i-th row contains
% information for node i in the following form:
%  [x-coordinate-of-node   y-coordinate-of-node]
% Coordinates are given in meters so that our units
% are compatible.

base       = 10;                      % half the length of the base
side_angle = 80.0*(pi/180.0);         % angle between base and sides
height     = base * tan(side_angle);  % height of the transmission tower
temp       = 1.0/ tan(side_angle);    % cotan(side_angle)
arm        = base;                    % length of the top horizontal side of the hangover 
node ...
= [ -base   0;
     base   0;
     -temp*5.*height/6.               1.*height/6.;
      temp*5.*height/6.               1.*height/6.;
     -temp*4.*height/6.               2.*height/6.;
      temp*4.*height/6.               2.*height/6.;
     -temp*3.*height/6.               3.*height/6.;
      temp*3.*height/6.               3.*height/6.;
     -temp*2.*height/6.               4.*height/6.;
      temp*2.*height/6.               4.*height/6.;
     -temp*1.*height/6.-arm           5.*height/6.
     -temp*1.*height/6.               5.*height/6.;
      temp*1.*height/6.               5.*height/6.;
      temp*1.*height/6.+arm           5.*height/6;
      0.0                             height;
     -temp*1.*height/6.-arm           4.*height/6;
     -temp*1.*height/6.-2.*arm        4.*height/6;
      temp*1.*height/6.+arm           4.*height/6;
      temp*1.*height/6.+2.*arm        4.*height/6];
          
     
% plot the undeformed truss
truss_plot(bar, node, 1)

% area is an m array. The i-th element contains the cross sectional
% area of bar i. We assume that all bars have rectangular cross
% section of size 0.01 [m] times 0.05 [m], i.e., area = 0.0005 [m^2].
area = 0.0005*ones(m,1);

% young is an m array. The i-th element contains the Young's
% modulus for bar i. We assume that all bars are made of the same 
% material with Young's modulus is 195 GPa 
young = 195.e9*ones(m,1);  

% set indices of fixed displacements
fixed = [1 2 3 4 ];

% determine indices of free displacements
free = [];
for i = 1:2*n
    if (~any(fixed==i))
       free = [free; i];
    end
end  

% determine the stiffness matrix
[K] = stiff( bar, node, area, young, free);

% factor the stiffness matrix
[K, ipivt, iflag] = lu_pp( K );



disp('Hit return to apply a load.......'); pause;

% right hand side
b     = zeros(2*n,1);
b(34)  = 80000*9.80665;
b(38)  = 80000*9.80665;

u                 = b(free);
[u, iflag]        = lu_pp_sl( K, u, ipivt );
if( iflag ~= 0 )
    error( ['    lu_pp_sl returned with iflag = ', int2str(iflag)])
end

% determine the locations of the node of the deformed truss
dnode = node;
for i = 1:size(free(:),1)
    j = free(i);
    if( mod(j,2) == 0 ) 
        % j is even, i.e., u(i) is vertical displacement of node j/2
        dnode(j/2,2) = node(j/2,2) - u(i);
    else
        % j is odd, i.e., u(i) is horizontal displacement of node (j+1)/2
        dnode((j+1)/2,1) = node((j+1)/2,1) + u(i);
    end
end

% plot the undeformed and the deformed truss
truss_plot(bar, node, 0, dnode)

disp('Hit return to apply another load.......'); pause;
% TRY OUT A NEW LOADING NOW!
b     = zeros(2*n,1);
b(34)  = 80000*9.80665;
b(38)  = 80000*9.80665;
for i=3:n
   b(2*i-1) = 20000;
end

u                 = b(free);
[u, iflag]        = lu_pp_sl( K, u, ipivt );
if( iflag ~= 0 )
    error( ['    lu_pp_sl returned with iflag = ', int2str(iflag)])
end

% determine the locations of the node of the deformed truss
dnode = node;
for i = 1:size(free(:),1)
   j = free(i);
  if( mod(j,2) == 0 ) 
      % j is even, i.e., u(i) is vertical displacement of node j/2
     dnode(j/2,2) = node(j/2,2) - u(i);
 else
    % j is odd, i.e., u(i) is horizontal displacement of node (j+1)/2
    dnode((j+1)/2,1) = node((j+1)/2,1) + u(i);
 end
end

% plot the undeformed and the deformed truss
truss_plot(bar, node, 0, dnode)