User:Eml4500.f08.team.foskey.ckf/hw6

MATLAB ASSIGNMENT

Electric Pylon Problem Overview
This assignment required the use of a previously created MATLAB file called electric_pylon.m. This code can be seen below. This code creates a two dimensional truss system that resembles an electric pylon. Although the code for the truss system is complex, the code only plots the truss using defined nodal coordinates.

The first part of the assignment involved scaling the height of the truss to a max height of 60m, applying a 1000N load, and calculating the stresses in the elements. The elements under the highest tensile and compressive stresses are identified. The setup of the problem is seen below. In the figure, the truss plot was created using the provided code, and the applied load P was added manually.



The second part of the assignment required the masses of the elements to be accounted for. The half of the mass of each element was lumped to each end. The eigenvalues and eigenvectors were found, and the lowest three were plotted. Also the lowest three vibrational periods were found.

Provided MATLAB Code
The following MATLAB code was provided.

An Explanation Of The Provided Code
Essentially, all the provided code does is plot the truss system. The code is long due to the fact that there are over 90 elements that make up the truss. However, the code is relatively simple.

The code begins by declaring there to be 46 global nodes, and 91 elements. Then is calculates the total number of degrees of freedom. This is accomplished by multiplying the number of nodes by two, because there are two degrees of freedom (horizontal and vertical) for each node.

Next, the code defines the global nodal locations into a position matrix using x-y-z coordinates. For example, the line position(:, 1) = [ 1.5; 0; 0]; gives the location of node 1 as being at the coordinate (1.5, 0, 0). This is repeated for all 46 global nodes.

Then, the code organizes the locations of the nodes into an easier to understand (and use) form. This is done by looping through all of the nodes defining x, y, and z variables for each node. For example, the line y(i) = position(2,i); will define the y coordinate of global node "i" as the second value of the position matrix (which is the y coordinate) in row "i" (which represents node "i"). Again, this is done for all 46 global nodes.

Once the nodal coordinates are defined, the code then assigns a node number to each end of each element. This is done by using a line of code such as node_connect(1, 89) = 39;. In this example line of code, the local node 1 of element 89 is defined as being equal to global node 39. Two lines of code correspond to each element - one line for each local node. This is done for all 91 elements.

Next, the code plots each element. This is done by looping through all of the elements. In each loop, the code obtains the global node numbers for each end of an element, defining these as "node_1" and "node_2". A vector is then defined for each direction going from one node to the other. For example in the line of code xx = [x(node_1),x(node_2)]; , the vector "xx" is defined as a vector going from the x value for node_1 to the x_value of node_2. Once the three vectors ( one each for the x, y and z directions) are defined for an element, they are plotted. This plot will represent one element. The plotting process is repeated for all 91 elements, until the entire truss is plotted.

The provided code ends by titling the plot, and the axis. It should also be noted that the code is designed to be able to plot a three dimensional truss, however all of the z coordinates for the nodes are zero. Therefore the code plots a two dimensional truss.

Additional Functions
The following functions are used in the following codes to solve for certain parts of the truss problem, particularly the loading code.

PlaneTrussElement calculates elemental stiffness matrices NodalSoln calculates the nodal displacements and the reactions of the truss PlaneTrussResults calculates the strain, stress, and axial force in each element

A more complete explanation of these functions can be found here.

PlaneTrussElement.m 

NodalSoln.m

PlaneTrussResults.m 

The function TransientPlaneTrussElement is used in the lumped mass code. The function creates an elemental stiffness matrix(like the function PlaneTrussElement ), as well as creates a mass matrix for the element. This function was taken from page 568 in the textbook.

TransientPlaneTrussElement.m

MATLAB Code For Loading The Pylon
The following MATLAB code is used to apply a 1000N downward load to the far right tip of the truss. The unloaded and loaded trusses are plotted, and the stresses in the elements are evaluated. The code is commented throughout, and explained further below.

An Explanation Of The Loading Code
The pylon generated by the originally provided code is to be loaded with a 1000N downward force on the right tip. The original, unloaded truss is to be plotted with the deformed truss, and the elements undergoing the highest stresses (tensile and compressive) are to be marked.

The code for loading the pylon can be broken down into three sub-parts: 1) converting the provided code to a usable form for analyzing the truss system, 2) plotting the trusses, and 3) determining the maximum stresses.

In the first part of the loading code (given the file name "solvingpylon1.m") takes the data from the provided code and converts it into variables that will work in a format of solving the truss similar to earlier MATLAB projects (NOTE: see The Two Bar Truss, The Six Bar Truss , and The Three Bar Space Truss) The loading code begins by running the provided code to get all of the data that that code generates. This eliminates the need to retype all of the original code material into a new code. Next, the loading code determines the highest point of the originally created truss and determines a scaling factor. The highest point is determined by setting a dummy variable "keep" equal to zero initially, and then checking the y coordinate value of each node to see if it is greater than the "keep" value. If it is, then that value becomes the new "keep" value, until all nodes have been checked. The scaling factor was found using the form $$ Scaling Factor = \frac{60}{keep} $$. Every node coordinate is then multiplied by the scaling factor, increasing the size of the pylon to a maximum height of 60m.

Next, the loading code is given the parameters of the elements. Each element has a Young's Modulus of 200GPa, and a cross sectional area of 4cm2. The load P is also entered as -1000, because the load is downward (in the negative y direction). The code then creates a 91x2 matrix of the nodal coordinates, with each node represented by a row, the x coordinates in the first column, and the y coordinates in the second row. The total number of displacements is also determined by multiplying the number of nodes by 2, which is the number of degrees of freedom for each node.

The loading program then creates the connectivity array by looping through each element and filling the connectivity array conn with the associated values from the node_connect matrix of the provided code. The location matrix master array lmm is then created. This is accomplished by using the node_connect matrix form the provided code. To fill lmm, one has to realize that each node has two associated degrees of freedom. Using this, the matrix is filled in for elements 1 and 2 like so:

(Note: Only the bold values are actually put into the lmm matrix)

This is accomplished by taking the global node number, decreasing it by one, doubling it, and adding 1 for the horizontal dof, or 2 for the vertical. For example, in element 1, the global nodes associated with it are 1 and 3, and the global dofs are 1, 2, 5 and 6. Looking at the code:

d1= 2*(node_connect(1,i)-1)+1; d2= 2*(node_connect(1,i)-1)+2; d3= 2*(node_connect(2,i)-1)+1; d4= 2*(node_connect(2,i)-1)+2;

lmm(i,1)=d1; lmm(i,2)=d2; lmm(i,3)=d3; lmm(i,4)=d4;

and realizing that node_connect(1,i) will give the global node 1, and that node_connect(2,i) will give 3 (the other global node), it follows that d1, d2, d3 and d4 will equal 1, 2, 5, and 6 respectively. This makes the first row of lmm to be $$ lmm=\begin{bmatrix} 1 & 2 & 5 & 6 \\ \end{bmatrix} $$

From there, the code uses parts of older codes (see note above) to constrain nodes 1 and 2, and to set up the force matrix R. The code then calculates the elemental stiffness matrices and the global stiffness matrix by using the PlaneTrussElement function. The code then calculates the strain, stress and axial force in each element using the NodalSoln and PlaneTrussResults functions.

The second part of the loading code plots the trusses. The code begins by closing the plot created by the provided code. This is why, when running the code, a plot will appear and quickly disappear. The code then creates a matrix of the displaced nodal coordinates defnodes by adding the displacements found in part 1 to the original nodal coordinates. The code then plots the original truss by creating vectors for each element that go from the node at one end, to the node at the other end of an element. This process is repeated for the deformed truss, this time using the coordinates of the dislocated nodes.

In the third part of the loading code, the results matrix from part 1 is examined. Since the second column of the results matrix represents the stress in an element, that value is examined. The stress value of the first element is examined first, and it is checked to be a compressive force or a tensile force (a compressive stress will be negative, and a tensile stress will be positive). The value is then stored the highest stress. The next element's stress is examined, and determined if it is compressive or tensile. The value is then compared to the highest stress. If it is of greater magnitude, it becomes the highest stress, and its element number is remembered. This process continues until all elements have been checked.

The code ends by plotting the elements under the highest stresses in different colors so that they might be easily identified, using the same process explained earlier. The plot is titled, axis are labeled and a legend is created.

MATLAB Code For Lumped Mass
The following code is used to complete the second portion of the assignment. The code is commented throughout, and a more detailed explanation of the code will follow.

An Explanation Of The Lumped Mass Code
The lumped mass code (solvingpylon2.m) begins by running the loading code (solvingpylon1.m) to receive the already calculated values that will be needed. These values include nodal coordinates, and element properties such as cross sectional area and modulus of elasticity. The density of the elements is then given in g/cm3, and then converted to kg/m3. The length of each element is calculated, and the mass of each element is found by multiplying the element's length by the cross sectional area and by the density.

The mass of each element is then divided by two, and each half mass is added to the global node represented by each end of the element. This is done for all of the elements until the entire mass of the truss is concentrated at the nodes.

A mass matrix is then created by using the function TrainsientPlaneTrussElement.

A reduced mass matrix M_bar, and a reduced stiffness matrix K_bar is then created using MATLAB's setdiff function. This function eliminates the values of the matrices that correspond to the constrained points. This is valid due to the principle of virtual work. The eigenvalues and eigenvectors are then found using MATLAB's eig function.

The natural frequencies $$ \omega $$, vibrational frequencies f, and the vibrational periods T are then found using the following formulas: $$ \omega = \sqrt{\lambda} $$

$$f = \frac{\omega}{2\pi}$$

$$T = \frac{1}{f} $$

where $$ \lambda $$ is an eigenvalue.

The code then proceeds to find the three lowest eigenvalues. This is accomplished by setting variables for the lowest eigenvalues (low1, low2, and low3) equal to large numbers, and then looping through all of the eigenvalues, replacing low1, low2, and low3 with any eigenvalues that are encountered that are lower than the previous values. While doing this, the code keeps track of which rows the lowest eigenvalues were in.

The code then organizes the eigenvectors. The calculated eigenvectors are contained in an 88x88 matrix that is very cumbersome. Each column of the eigenvector matrix represents one eigenvector. Since the assignment required the finding of the lowest three eigenvalue-eignevector pairs, the needed eigenvectors are found in the columns that correspond to the rows in which the lowest eigenvalues were in. The code extracts the desired eigenvectors and puts them into column matrices called modes1, modes2, and modes3.

In a similar manner as the determining of the lowest eigenvalues, the lowest three vibrational periods are found.

Since there was a number of things calculated by this code, with a majority of the values being large matrices, only a selected amount of data was deemed important for output. The next section of code outputs the important data. This data consists of the eigenvalues, the three lowest eigenvalues and their corresponding eigenvectors, and the three lowest vibrational periods.

To finish the code, three plots are created. Each plot consists of an undeformed and a deformed truss. The undeformed truss is plotted using the same lines of code from the loading code. The deformed truss is plotted using lines of code similar to the code from the loading code, but with adjustments to the coordinates of the nodes. For this situation, the "deformed" nodal positions are created by adding the values from an eigenvector to the original nodal coordinates. This is done because the eigenvectors represent the natural displacements for the truss system. By adding the natural displacements to the original nodal coordinates, the result is the naturally deformed truss. This is how the truss system will move, even under no loads, due to the natural frequencies associated with the elements. These natural displacements, however, are very small. In order for the natural displacements to be seen on the plot, a scaling factor was used. This scaling factor was 100,000.

Three plots are created from the lumped mass code. Each of these plots represents the deformations that occur as a result of one of the three lowest eigenvectors.

Results
There are no outputted results from the originally provided code.

Running the loading code (solvingpylon.m) generates the following results

Where d is the displacement matrix that shows the displacements for each degree of freedom in meters reactions are the forces at each of the constrained nodes results shows the strain, stress, and axial forces in each element.

trans_elem is the element number of the element under the highest tensile stress, which is max_tens The highest tensile stress encountered was 9.0511 x 106 Pa, or 9.0511 MPa comp_elem is the number of the element under the highest compressive force, which is called max_comp The highest compressive stress encountered was 8.6957 x 106 Pa, or 8.6957 MPa

In the results matrix, the strain, stress, and axial force values are contained in columns 1, 2 and 3 respectively. The rows correspond to the element numbers. For example, elements 1 and 2 are given as follows:

Running the lumped mass code (solvingpylon2.m) results in the following output

Where eigenvalues is a column matrix that is a diagonal matrix of the matrix of eigenvalues low1, low2,& low3 are the three lowest eigenvalues (low1 being the lowest, low3 the highest) modes1, modes2, and modes3 are the three lowest eigenvectors (modes1 being the lowest, modes3 the highest) and lowT1, lowT2, & lowT3 are the three lowest vibrational periods (lowT1 being the lowest, lowT3 being the highest)

Figure From Loading Code
The following figure is created when the M-file solvingpylon1.m is run. 

The figure shows the original truss (scaled to a height of 60m) in red dotted lines. The deformed truss is shown in blue solid lines. Arrows point out the elements under the highest stresses. The element under the highest tensile stress is shown in light blue, and the element under the highest compressive stress is shown in green.

Figures From Lumped Mass Code
The following three figures are generated when the lumped mass code is run(solvingpylon2.m). Each figure represents the displacements from one of the three lowest eigenvectors. In each figure, the red dotted lines represent the undeformed truss, and the solid blue lines represent the deformed truss. The x and y axis for all three figures are in terms of meters.   

300 M Steel
300 M steel is a type of steel that can be considered a modification of AISI 4340 steel. 300 M steel contains silicon, vanadium and has a greater carbon and molybdenum content than AISI 4340 steel does. 300 M steel is considered to have good strength, toughness, ductility and fracture strength.

The strength of 300 M steel is 280 - 305 ksi (1929 -2100 MPa).

300 M steel is commonly used for applications that require strength near 300 ksi, such as aircraft landing gear.

For more information see: 300M Alloy Steel Material Property Sheet - Metal Suppliers Online

In the electric pylon problem, the highest stress encountered was the tensile stress. It was found to be 9.0511 x 106 Pa, or 9.0511 MPa. If the pylon was constructed of 300 M steel, then the maximum stress is well below the strength of the material.

Since the strength of the material is so much greater than the maximum stress, two observation can be made. First, if the pylon was constructed of 300 M intentionally, then loads much higher than the 1000N force are expected (or combinations of low loads) that would result in much higher stresses encountered in this situation. Second, if 300 M was not specifically required and was chosen without much reason, and if the pylon experiences only low levels of stress, then a different material should have been chosen instead of 300 M steel. A lower grade of steel would have still adequately been able to withstand the applied loads, and would have been less expensive than 300 M steel.