User:Eml4500.f08.lulz.strack/hw4

MATLAB Code: Five-bar and Three-bar Trusses
The first part of the required MATLAB code constructs and loads a truss system composed of five bars. The code below represents Example 1.4 (FiveBarTrussEx.m) on page 25 of the textbook. The first part defines the specifications of each bar and node including modulus of elasticity (e), cross-sectional area (a), and loading values (P). The "nodes" matrix defines the absolute locations of the four nodes. The "conn" matrix defines the global node numbers at the ends of each bar or element. Each row represents an element and nodes are numbered by the lowest first. As shown in the code, there are 5 sets of 2 numbers for the 2 nodes corresponding to each element. The "lmm" matrix represents the global degrees of freedom (dof) for each element. Each node has 2 dof's; therefore, each element has 4 corresponding dof's. Degrees of freedom are numbered from 1 to 2*(number of nodes) where the first 2 correspond to Global Node 1, the second 2 correspond to the Global Node 2, etc. From all this data, the global stiffness matrix (K) is created. The displacement of the loaded node is calculated using the relevant parts of the global stiffness matrix, inverting them, and multiplying by the load array (R). Finally, the reactions (reactions) at each node are calculated by multiplying the displacement array (d) by the global stiffness matrix. The functions "PlaneTrussElement.m", "NodalSoln.m", and "PlaneTrussResults.m" are used again from Section 6 of Homework 2.

FiveBarTrussEx.m
% Five bar truss example e1=200*10^3; e2=70*10^3; a1=40*100; a2=30*100; a3=20*100; P = 150*10^3; nodes = 1000*[0, 0; 1.5, 3.5; 0, 5; 5, 5]; conn = [1,2; 2,4; 1,3; 3,4; 2,3]; lmm = [1,2,3,4; 3, 4, 7, 8; 1, 2, 5, 6; 5, 6, 7, 8; 3, 4, 5, 6]; K=zeros(8); % Generate stiffness matrix for each element and assemble it. for i=1:2 lm=lmm(i,:); con=conn(i,:); k=PlaneTrussElement(e1, a1, nodes(con,:)); K(lm, lm) = K(lm, lm) + k; end for i=3:4 lm=lmm(i,:); con=conn(i,:); k=PlaneTrussElement(e1, a2, nodes(con,:)); K(lm, lm) = K(lm, lm) + k; end

lm=lmm(5,:); con=conn(5,:); k=PlaneTrussElement(e2, a3, nodes(con,:)); K(lm, lm) = K(lm, lm) + k

% Define the load vector R = zeros(8,1); R(4)=-P

% Nodal solution and reactions [d, reactions] = NodalSoln(K, R, [1,2,7,8], zeros(4,1)) results=[]; for i=1:2 results = [results; PlaneTrussResults(e1, a1, ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end for i=3:4 results = [results; PlaneTrussResults(e1, a2, ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results = [results; PlaneTrussResults(e2, a3, ...       nodes(conn(5,:),:), d(lmm(5,:)))]

Output to the Command Window

>> FiveBarTrussEx

K =

32600       76067       -32600       -76067            0            0            0            0        76067  2.9749e+005       -76067 -1.7749e+005            0    -1.2e+005            0            0 -32600      -76067  2.4309e+005  1.1914e+005       -32998        32998 -1.7749e+005       -76067 -76067 -1.7749e+005 1.1914e+005  2.4309e+005        32998       -32998       -76067       -32600 0           0       -32998        32998    1.53e+005       -32998    -1.2e+005            0 0   -1.2e+005        32998       -32998       -32998    1.53e+005            0            0 0           0 -1.7749e+005       -76067    -1.2e+005            0  2.9749e+005        76067 0           0       -76067       -32600            0            0        76067        32600

R =

0          0           0     -150000           0           0           0           0

d =

0           0      0.53895     -0.95306       0.2647      -0.2647            0            0

reactions =

54927 1.5993e+005 -54927     -9926.7

results =

-0.0001743     -34.859 -1.3944e+005 -3.15e-005     -6.2999       -25200 -5.2941e-005     -10.588       -31764 -5.2941e-005     -10.588       -31764 0.00032087      22.461        44922

The results show the strain, stress, and axial force in each element corresponding to the three columns in the "results" matrix.

To visualize the loaded and unloaded system, the following code graphs the system in its loaded and unloaded state. In order to be able to realize the direction of node movement, the unfixed node displacements have been magnified by a factor of two (using the variable "m"). As shown, the original values from the displacement array (d3-d6) above were used below to get the right proportions of the nodal displacements.

FiveBarTrussDisplay.m
%Displaying the five bar truss clear; close;

m = 2;

d1 = 0;         %disp of node 1 in the x direction d2 = 0;         %disp of node 1 in the y direction d3 = 0.53895*m; %disp of node 2 in the x direction d4 = -0.95306*m; %disp of node 2 in the y direction d5 = 0.2647*m;  %disp of node 3 in the x direction d6 = -0.2647*m; %disp of node 3 in the y direction d7 = 0;         %disp of node 4 in the x direction d8 = 0;         %disp of node 4 in the y direction

%input nodal coordinates num_nodes = 4;              %total number of nodes num_elems = 5;              %total number of elements tot_dofs = 2 * num_nodes;   %total dofs of the system

%undeformed positions pos(1,:)=[0 1.5 0 5]; pos(2,:)=[0 3.5 5 5];

%undeformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

%element connectivity array %follows: %node_conn(local node number, element number) = global node number

node_conn(1,1)=1;    %element 1 node_conn(2,1)=2;

node_conn(1,2)=2;    %element 2 node_conn(2,2)=4;

node_conn(1,3)=1;    %element 3 node_conn(2,3)=3;

node_conn(1,4)=3;    %element 4 node_conn(2,4)=4;

node_conn(1,5)=2;    %element 5 node_conn(2,5)=3;

%plot the undeformed 2-bar truss system for i=1:num_elems node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; axis([-1 6 -1 6]) plot(xx,yy,'--rs') hold on end

n = -0.5; o = -0.25; text(0+n,0,'GN1','HorizontalAlignment','center') text(1.5+n,3.5,'GN2','HorizontalAlignment','center') text(0+n,5,'GN3','HorizontalAlignment','center') text(5+n,5,'GN4','HorizontalAlignment','center') text(0.75+o,1.75,'E1','HorizontalAlignment','center') text(3.25+o,4.25,'E2','HorizontalAlignment','center') text(0+o,2.5,'E3','HorizontalAlignment','center') text(2.5+o,5,'E4','HorizontalAlignment','center') text(0.75+o,4.25,'E5','HorizontalAlignment','center')

%deformed positions pos(1,:)=[0+d1 1.5+d3 0+d5 5+d7]; pos(2,:)=[0+d2 3.5+d4 5+d6 5+d8];

%deformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

text(1.5+d3+n,3.5+d4,'GN2','HorizontalAlignment','center') text(0+d5+n,5+d6,'GN3','HorizontalAlignment','center') text(pos(1,1)+(pos(1,2)-pos(1,1))/2+o,pos(2,1)+(pos(2,2)-pos(2,1))/2,'E1','HorizontalAlignment','center') text(pos(1,2)+(pos(1,4)-pos(1,2))/2+o,pos(2,2)+(pos(2,4)-pos(2,2))/2,'E2','HorizontalAlignment','center') text(pos(1,1)+(pos(1,3)-pos(1,1))/2+o,pos(2,1)+(pos(2,3)-pos(2,1))/2,'E3','HorizontalAlignment','center') text(pos(1,3)+(pos(1,4)-pos(1,3))/2+o,pos(2,3)+(pos(2,4)-pos(2,3))/2,'E4','HorizontalAlignment','center') text(pos(1,2)+(pos(1,3)-pos(1,2))/2+o,pos(2,2)+(pos(2,3)-pos(2,2))/2,'E5','HorizontalAlignment','center')

%plot the deformed 2-bar truss system for i=1:num_elems hold on   node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; plot(xx,yy,'-s') end

title('Deformed and Undeformed Five-bar Truss System') xlabel('X') ylabel('Y')



The dotted line represents the undeformed system shape. The solid line represents the deformed system shape. The abbreviation "GN" and "E" stand for Global Node and Element, respectively. Note only nodes 2 and 3 are displaced because they are unfixed.

For the three bar truss created in the class notes, the two bar truss system code (Two Bar Truss System.m) was modified from Section 6 of Homework 2. One node and element each were added to the system, and e, A, P, and L values were changed to coincide with the problem.

ThreeBarTrussEx.m
% Three bar truss example clear all; e = [2 4 3]; A = [3 1 2]; P = 30; L = [5 5 10]; alpha = pi/6; beta = -pi/6; gamma = -3*pi/4; nodes = [0, 0; L(1)*cos(alpha), L(1)*sin(alpha); L(1)*cos(alpha) + L(2)*cos(beta), L(1)*sin(alpha) + L(2)*sin(beta); L(1)*cos(alpha) + L(3)*cos(gamma), L(1)*sin(alpha) + L(3)*sin(gamma)]; dof = 2*length(nodes); conn = [1 2; 2 3; 2 4]; lmm = [1 2 3 4; 3 4 5 6; 3 4 7 8]; elems = size(lmm,1); K = zeros(dof); R = zeros(dof,1); debc = [1 2 5 6 7 8]; ebcVals = zeros(length(debc),1);

%load vector R = zeros(dof,1); R(4) = P;

% Assemble global stiffness matrix K = zeros(dof); for i = 1:elems lm = lmm(i,:); con = conn(i,:); k_local = e(i)*A(i)/L(i)*[1 -1; -1 1] k = PlaneTrussElement(e(i), A(i), nodes(con,:)) K(lm, lm) = K(lm, lm) + k; end K R % Nodal solution and reactions [d, reactions] = NodalSoln(K, R, debc, ebcVals) results = []; for i = 1:elems results = [results; PlaneTrussResults(e, A, nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results

Output to the Command Window

>> ThreeBarTrussEx

k_local =

1.2000  -1.2000   -1.2000    1.2000

k =

0.9000   0.5196   -0.9000   -0.5196    0.5196    0.3000   -0.5196   -0.3000   -0.9000   -0.5196    0.9000    0.5196   -0.5196   -0.3000    0.5196    0.3000

k_local =

0.8000  -0.8000   -0.8000    0.8000

k =

0.6000  -0.3464   -0.6000    0.3464   -0.3464    0.2000    0.3464   -0.2000   -0.6000    0.3464    0.6000   -0.3464    0.3464   -0.2000   -0.3464    0.2000

k_local =

0.6000  -0.6000   -0.6000    0.6000

k =

0.3000   0.3000   -0.3000   -0.3000    0.3000    0.3000   -0.3000   -0.3000   -0.3000   -0.3000    0.3000    0.3000   -0.3000   -0.3000    0.3000    0.3000

K =

0.9000   0.5196   -0.9000   -0.5196         0         0         0         0    0.5196    0.3000   -0.5196   -0.3000         0         0         0         0   -0.9000   -0.5196    1.8000    0.4732   -0.6000    0.3464   -0.3000   -0.3000   -0.5196   -0.3000    0.4732    0.8000    0.3464   -0.2000   -0.3000   -0.3000         0         0   -0.6000    0.3464    0.6000   -0.3464         0         0         0         0    0.3464   -0.2000   -0.3464    0.2000         0         0         0         0   -0.3000   -0.3000         0         0    0.3000    0.3000         0         0   -0.3000   -0.3000         0         0    0.3000    0.3000

R =

0    0     0    30     0     0     0     0

d =

0        0  -11.6737   44.4051         0         0         0         0

reactions =

-12.5672  -7.2557   22.3866  -12.9249   -9.8194   -9.8194

results =

2.4186      4.8371       9.6742       7.2557       14.511       9.6742       14.511       6.4625       12.925        25.85       19.387       38.775        25.85       38.775       2.3145       4.6289       9.2578       6.9434       13.887       9.2578       13.887

After running this code, the output produces displacement of node 2 to be -11.6737 in the x-direction and 44.4051 in the y-direction. On the other hand, the "results" matrix has 7 columns instead of 4. This is due to the fact that vectors for the modulus of elasticity (e), cross-sectional area (A), and element length (L) were passed as variables instead of the scalar values. The first column is the strain array for elements 1-3. This is followed by two 3x3 matrices created for stresses and axial forces of the elements instead of two arrays of length 3 each. The real values for stress and axial force lie along the diagonals of each of these matrices.

The code for displaying the three bar truss was similar to the code displaying the five bar truss. However, a displacement magnification factor (m) of 0.25 was used to ease visualization of the deformation.

ThreeBarTrussDisplay.m
%Displaying the three bar truss clear; close;

m = 0.25;

d1 = 0;         %disp of node 1 in the x direction d2 = 0;         %disp of node 1 in the y direction d3 = -11.6737*m; %disp of node 2 in the x direction d4 = 44.4051*m; %disp of node 2 in the y direction d5 = 0;         %disp of node 3 in the x direction d6 = 0;         %disp of node 3 in the y direction d7 = 0;         %disp of node 4 in the x direction d8 = 0;         %disp of node 4 in the y direction

%input nodal coordinates num_nodes = 4;              %total number of nodes num_elems = 3;              %total number of elements tot_dofs = 2 * num_nodes;   %total dofs of the system

%element 1 details L1 = 5; E1 = 2; A1 = 3; theta1 = 30; %deg thetar1 = (pi * theta1)/180;

%element 2 details L2 = 5; E2 = 4; A2 = 1; theta2 = -30; %deg thetar2 = (pi * theta2)/180;

%element 3 details L3 = 10; E3 = 3; A3 = 2; theta3 = -135; %deg thetar3 = (pi * theta3)/180;

%undeformed positions pos(1,:)=[0,L1*cos(thetar1),L1*cos(thetar1)+L2*cos(thetar2),L1*cos(thetar1)+L3*cos(thetar3)]; pos(2,:)=[0,L1*sin(thetar1),L1*sin(thetar1)+L2*sin(thetar2),L1*sin(thetar1)+L3*sin(thetar3)];

%undeformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

%element connectivity array %follows: %node_conn(local node number, element number) = global node number

node_conn(1,1)=1;    %element 1 node_conn(2,1)=2;

node_conn(1,2)=2;    %element 2 node_conn(2,2)=3;

node_conn(1,3)=2;    %element 3 node_conn(2,3)=4;

%plot the undeformed 2-bar truss system for i=1:num_elems node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; axis([-5 10 -10 20]) plot(xx,yy,'--rs') hold on end

n = -1; o = -0.5; text(x(node_conn(1,1))+n,y(node_conn(1,1)),'GN1','HorizontalAlignment','center') text(x(node_conn(2,1))+n,y(node_conn(2,1)),'GN2','HorizontalAlignment','center') text(x(node_conn(2,2))+n,y(node_conn(2,2)),'GN3','HorizontalAlignment','center') text(x(node_conn(2,3))+n,y(node_conn(2,3)),'GN4','HorizontalAlignment','center') text(x(node_conn(2,1))/2+o,y(node_conn(2,1))/2,'E1','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,2))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,2))-y(node_conn(2,1)))/2,... 'E2','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,3))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,3))-y(node_conn(2,1)))/2,... 'E3','HorizontalAlignment','center')

%deformed positions pos(1,:)=[0+d1,L1*cos(thetar1)+d3,L1*cos(thetar1)+L2*cos(thetar2)+d5,L1*cos(thetar1)+L3*cos(thetar3)+d7]; pos(2,:)=[0+d2,L1*sin(thetar1)+d4,L1*sin(thetar1)+L2*sin(thetar2)+d6,L1*sin(thetar1)+L3*sin(thetar3)+d8];

%deformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

text(x(node_conn(2,1))+n,y(node_conn(2,1)),'GN2','HorizontalAlignment','center') text(x(node_conn(2,1))/2+o,y(node_conn(2,1))/2,'E1','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,2))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,2))-y(node_conn(2,1)))/2,... 'E2','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,3))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,3))-y(node_conn(2,1)))/2,... 'E3','HorizontalAlignment','center')

%plot the deformed 2-bar truss system for i=1:num_elems hold on   node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; plot(xx,yy,'-s') end

title('Deformed and Undeformed Three-bar Truss System') xlabel('X') ylabel('Y')



The dotted line represents the undeformed system shape. The solid line represents the deformed system shape. The abbreviation "GN" and "E" stand for Global Node and Element, respectively. Note only node 2 is displaced because it is unfixed.