User:Eml4500.f08.bottle.loschak/HW4

HW 4
The code used to solve the two-bar truss system was altered to solve the three-bar truss system. The altered code is displayed below, followed by the functions required to solve the system, along with a diagram near the end portraying the deformation in the system. % Three bar truss example taken from 2 bar truss example % modified for three bars instead clear all; % Matrices e, A, and L represent Young's Modulus, cross sectional area, and % length of the three elements e = [2 4 3]; A = [3 1 2]; L=[5 5 10]; % Verical force exerted on node 2 P = 30; % Angles between each element and the horizontal alpha = pi/6; beta = pi/6; ceta = pi/4;  % note: "ceta" is not actually a greek letter % Coordinates of each of the 4 nodes calculated with trigonometry nodes = [0, 0; L(1)*cos(alpha), L(1)*sin(alpha); 2*L(1)*cos(alpha), 0; L(1)*cos(alpha)-L(3)*cos(ceta),L(1)*sin(alpha)-L(3)*sin(ceta)]; % Since the problem contains 4 nodes, there are 8 global degrees of freedom dof=2*length(nodes); % The connectivity matrix describes which nodes are connected to one % another by listing each row as an element and the global nodes it % connects conn=[1,2; 2,3; 2,4]; % The location matrix master array details which global degrees of freedom % are related to matching local degrees of freedom. lmm = [1, 2, 3, 4; 3, 4, 5, 6; 3, 4, 7, 8]; % The number of elements is determined by the size of the lmm matrix elems=size(lmm,1); % The K (square) and R (columnar) matrices are initialized with zeros 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 The above code uses the following functions obtained from the website highlighted in the Fundamental Finite Element Analysis and Applications textbook by Bhatti. function k = PlaneTrussElement(e, A, coord) % PlaneTrussElement(e, A, coord) % Generates stiffness matrix for a plane truss element % e = modulus of elasticity % A = area of cross-section % coord = coordinates at the element ends x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2); L=sqrt((x2-x1)^2+(y2-y1)^2); ls=(x2-x1)/L; ms=(y2-y1)/L; k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms; ls*ms, ms^2,-ls*ms,-ms^2; -ls^2,-ls*ms,ls^2,ls*ms; -ls*ms,-ms^2,ls*ms,ms^2]; function [d, rf] = NodalSoln(K, R, debc, ebcVals) % [nd, rf] = NodalSoln(K, R, debc, ebcVals) % Computes nodal solution and reactions % K = global coefficient matrix % R = global right hand side vector % debc = list of degrees of freedom with specified values % ebcVals = specified values dof = length(R); df = setdiff(1:dof, debc); Kf = K(df, df); Rf = R(df) - K(df, debc)*ebcVals; dfVals = Kf\Rf; d = zeros(dof,1); d(debc) = ebcVals; d(df) = dfVals; rf = K(debc,:)*d - R(debc); function [eps, sigma, force] = PlaneTrussResults(e, A, coord, disps) % [eps, sigma, force] = PlaneTrussResults(e, A, coord, disps) % Compute plane truss element results % e = modulus of elasticity % A = Area of cross-section % coord = coordinates at the element ends % disps = displacements at element ends % The output quantities are eps = axial strain % sigma = axial stress and force = axial force. x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2); L=sqrt((x2-x1)^2+(y2-y1)^2); ls=(x2-x1)/L; ms=(y2-y1)/L; T=[ls,ms,0,0; 0,0,ls,ms]; d = T*disps; eps= (d(2)-d(1))/L; sigma = e.*eps; force = sigma.*A; Results of the solved 3-bar truss problem are as follows: k_local = 1.2        -1.2         -1.2          1.2 k = 0.9     0.51962         -0.9     -0.51962      0.51962          0.3     -0.51962         -0.3         -0.9     -0.51962          0.9      0.51962     -0.51962         -0.3      0.51962          0.3 k_local = 0.8        -0.8         -0.8          0.8 k = 0.6    -0.34641         -0.6      0.34641     -0.34641          0.2      0.34641         -0.2         -0.6      0.34641          0.6     -0.34641      0.34641         -0.2     -0.34641          0.2 k_local = 0.6        -0.6         -0.6          0.6 k = 0.3         0.3         -0.3         -0.3          0.3          0.3         -0.3         -0.3         -0.3         -0.3          0.3          0.3         -0.3         -0.3          0.3          0.3 K = 0.9     0.51962         -0.9     -0.51962            0            0            0            0      0.51962          0.3     -0.51962         -0.3            0            0            0            0         -0.9     -0.51962          1.8      0.47321         -0.6      0.34641         -0.3         -0.3     -0.51962         -0.3      0.47321          0.8      0.34641         -0.2         -0.3         -0.3            0            0         -0.6      0.34641          0.6     -0.34641            0            0            0            0      0.34641         -0.2     -0.34641          0.2            0            0            0            0         -0.3         -0.3            0            0          0.3          0.3            0            0         -0.3         -0.3            0            0          0.3          0.3 R = 0    0     0    30     0     0     0     0 d = 0           0      -11.674       44.405            0            0            0            0 reactions = -12.567     -7.2557       22.387      -12.925      -9.8194      -9.8194 results = 2.4186      6.4625       2.3145

Displaying the three-bar truss problem can be done with the code written below. clear; close; % The 3-bar truss problem has 8 connection points between lines that will % be displayed n_node = 8; % there are 9 dotted lines to be displayed, showing the 3-bar truss system % before deformation and showing the lines AB, AC, AD n_elem = 9; % P values calculated from the results above based on the equation % P = sqrt[(f1_1)^2 + (f2_1)^2] P2_1 = 14.512; P1_2 = 25.849; P1_3 = 13.88676; k_1 = 1.2; k_2 = 0.8; k_3 = 0.6; % Deformation in each bar calculated by P(e)/k(e) AB = P2_1/k_1; AC = P1_2/k_2; AD = P1_3/k_3; % Positions of all 8 nodes based on the origin at global node 4 position(:,1) = [10*cos(pi/4)-5*cos(pi/6); 10*sin(pi/4)-5*sin(pi/6); 0];  % global node 1 position(:,2) = [10*cos(pi/4); 10*sin(pi/4); 0];    % Point A, global node 2 position(:,3) = [10*cos(pi/4)+5*cos(pi/6); 10*sin(pi/4)-5*sin(pi/6); 0];     % global node 3 position(:,4) = [0; 0; 0]; % global node 4 position(:,5) = [10*cos(pi/4)+AB*cos(pi/6); 10*sin(pi/4)+AB*sin(pi/6); 0]; % point B position(:,6) = [10*cos(pi/4)-AC*cos(pi/6); 10*sin(pi/4)+AC*sin(pi/6); 0];  % point C position(:,7) = [10*cos(pi/4)+AD*cos(pi/4); 10*sin(pi/4)+AD*sin(pi/4); 0];  % point D position(:,8) = [-4.60; 51.48; 0];  % point E, or global node 2 after deformation % Creates x,y,z coordinates for the position of each node above for i = 1 : n_node; x(i) = position(1,i); y(i) = position(2,i); z(i) = position(3,i); end % node_connect(local node number, element number) = global node number node_connect(1,1) = 1; node_connect(2,1) = 2; node_connect(1,2) = 2; node_connect(2,2) = 3; node_connect(1,3) = 2; node_connect(2,3) = 4; node_connect(1,4) = 2; node_connect(2,4) = 5; node_connect(1,5) = 2; node_connect(2,5) = 6; node_connect(1,6) = 2; node_connect(2,6) = 7; node_connect(1,7) = 5; node_connect(2,7) = 8; node_connect(1,8) = 6; node_connect(2,8) = 8; node_connect(1,9) = 7; node_connect(2,9) = 8; % Plotting a dotted line between nodes for i = 1 : n_elem; node_1 = node_connect(1,i); node_2 = node_connect(2,i); xx = [x(node_1),x(node_2)]; yy = [y(node_1),y(node_2)]; zz = [z(node_1),z(node_2)]; plot3(xx,yy,zz,':') hold on end % Plotting a solid line to represent the deformed 3-bar truss system % Element 1 node_1 = 1; node_2 = 8; xx = [x(node_1),x(node_2)]; yy = [y(node_1),y(node_2)]; zz = [z(node_1),z(node_2)]; % Element 2 node_1 = 3; node_2 = 8; xxx = [x(node_1),x(node_2)]; yyy = [y(node_1),y(node_2)]; zzz = [z(node_1),z(node_2)]; % Element 3 node_1 = 4; node_2 = 8; xxxx = [x(node_1),x(node_2)]; yyyy = [y(node_1),y(node_2)]; zzzz = [z(node_1),z(node_2)]; plot3(xx,yy,zz,'-') plot3(xxx,yyy,zzz,'-') plot3(xxxx,yyyy,zzzz,'-') title('Deformed and undeformed three-bar truss system') xlabel('x') ylabel('y') zlabel('z') view([0 0 1])     % xy plane view

''Note: this MATLAB code was created using the following codes as examples:

- The MATLAB script for analyzing the Two Bar Truss System was obtained from Professor Vu-Quoc's website http://clesm.mae.ufl.edu/~vql/courses/fead/2008.fall/codes/twoBarTrussEx.txt

- The MATLAB script for the plot of a crab-leg structure created by X.G. Tan (2003) and obtained from Professor Vu-Quoc's website http://apollo.mae.ufl.edu/~vql/courses/fead/2006.fall/codes/matlab_plot_example2.m

- The MATLAB script for the plot of an electric pylon (undeformed five-bar truss structure) created by Team MeshWorks (2003) and obtained from Professor Vu-Quoc's website http://apollo.mae.ufl.edu/~vql/courses/fead/2006.fall/code.web.pages/electric_pylon.m''