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

Electric Pylon Frame System

Problem Statement
This problem uses the electric pylon code electric_pylon.mprovided earlier, and solves the system as a frame system. This problem is identical in concept to the earlier electric pylon problem, except that the earlier problem was a truss system and did not account for bending moments in the bars.

The assignment again 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,as well as the elements under the highest shear forces. The setup of the problem is seen below.



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

Additional Functions
This loading code uses three functions: PlaneFrameElement, PlaneFrameResults, and NodalSoln.

PlaneFrameElement
This function works by receiving an element's modulus of elasticity, cross sectional area, moment of inertia, and the coordinates of the elements endpoints. The length of the element is calculated based on the endpoint coordinates, and the director cosines are calculated using the length and coordinates.

Then the local element stiffness matrix k_ is found. This matrix is converted into its global coordinate system counterpart k by using the transformation matrix T. All three of these matrices are the matrices used for a frame system, which accounts for the bending moment degrees of freedom.

PlaneFrameResults
This function works by receiving an element's modulus of elasticity, cross sectional area, moment of inertia, the coordinates of the elements endpoints, and the displacements of the element. The length of the element is calculated based on the endpoint coordinates, and the director cosines are calculated using the length and coordinates. The displacement matrix d is created in local coordinates using the transformation matrix T on the global displacements. The strain, eps, is found by dividing the total element elongation by the element's original length. The stress, sigma, is found by multiplying the strain by the modulus of elasticity e. The axial force in the element is then found by multiplying the stress by the cross sectional area A of the element.

NodalSoln
The function NodalSoln reads in the global stiffness matrix, the applied loads, a matrix ,debc, containing the degrees of freedom with specified values (in this case the constrained degrees of freedom), and a matrix, ebcVals ,with the values of the nodes (in this case zeros). The function first calculates the number of global degrees of freedom. Through several steps, the function ultimately calculates the global displacements, and the reactions (in a global coordinate system) of the bars. The displacement matrix, and reaction matrix is returned to the original call function. This function is the same as the function used for the truss system because it is not dependent on the number of degrees of freedom at each node, only the overall total number of degrees of freedom.

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 frame is to be plotted with the deformed frame, 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 frame system, 2) plotting the frames, and 3) determining the maximum stresses.

In the first part of the loading code (given the file name "electricframe.m") takes the data from the provided code and converts it into variables that will work in a format of solving the frame similar to earlier MATLAB projects involving truss systems (see the original Electric Pylon Truss Problem) 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

$$ScalingFactor = \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. Since the problem is a frame system, a moment of inertia value is also needed. The bars were assumed to have square cross sections. This leads to calculating the moment of inertia to be

$$I = \frac{b*h^3}{12} = \frac{A^2}{12}$$

Where I is the moment of inertia, b is the base of the cross section, and h is the height. Since the cross section is a square, b and h are equal. The end result is one twelfth the square of the cross sectional area.

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 3, which is the number of degrees of freedom for each node( horizontal, vertical, and a rotation).

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 three 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, multiplying it by three, and adding 1 for the horizontal dof, 2 for the vertical, or 3 for the rotation. For example, in element 1, the global nodes associated with it are 1 and 3, and the global dofs are 1, 2, 3, 7, 8 and 9. Looking at the code:

%create location matrix master array using data from original code lmm=zeros(n_elem,6); for i=1:n_elem d1= 3*(node_connect(1,i)-1)+1; d2= 3*(node_connect(1,i)-1)+2; d3 = 3*(node_connect(1,i)-1)+3; d4= 3*(node_connect(2,i)-1)+1; d5= 3*(node_connect(2,i)-1)+2; d6= 3*(node_connect(2,i)-1)+3; lmm(i,1)=d1; lmm(i,2)=d2; lmm(i,3)=d3; lmm(i,4)=d4; lmm(i,5)=d5; lmm(i,6)=d6; end

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, d4, d5, and d6 will equal 1, 2, 3, 7, 8 and 9 respectively. This makes the first row of lmm to be

$$ lmm=\begin{bmatrix} 1 & 2 & 3 & 7 & 8 & 9 \\ \end{bmatrix} $$

From there, the code uses parts of older codes 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 PlaneFrameElement function. The code then calculates the strain, stress and axial force in each element using the NodalSoln and PlaneFrameResults functions.

The second part of the loading code plots the frames. 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 frame 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 frame, 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.

To find the highest nodal rotation, a loop check every third displacement from d, until the highest rotation and its correspond in node is found. Then the conn matrix is used to find the elements that touch that node. The displacements for those elements are found and converted to their axial equivalents. The displacement equations are differentiated twice to obtain the shear equations, which are then used to find the maximum shear in the elements touching the node with the highest rotation.

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.

Results Output
Running the loading code "electricframe.m" generates the following output in MATLAB

Comments On Output
Due to the fact that there are 91 elements in the frame system, the elemental stiffness matrices and the global stiffness matrix are not shown in the output. Displaying these matrices would have made the output far too long, and the global stiffness matrix would be too large to display on a computer screen in a logical manner. The matrix would be broken up and displayed in segments.

The output does display three matrices, the displacement matrix d, the reactions matrix reactions and the results matrix results.

In the displacement matrix, every third value represents a rotational displacement. The first six values of d are zero because nodes 1 and 2 were constrained.

The reactions matrix shows the forces acting at the constrained nodes. This information can be understood more easily in table format

The results matrix shows the strain, stress and axial forces in each element, and is again easier to understand in table form

The final few lines of output identify the elements under the highest stresses and roations. Based upon the output, the element under the highest tensile stress is element 77, which is under 2123 Pa of tensile stress. The element under the highest compressive stress is element 32, which is under 1.632 Pa of compressive stress.

"rotate_node" is the node with the highest rotational displacement (node 33), which is 0.14978 radians counter clockwise. The two elements touching that node are element 55 and 77. The shear force in element 55 is 1.5257e+009 N/m2 and the shear force in element 77 is 1.6221e+009 N/m2.

Resulting Figure
Running the above MATLAB code generates the following figure



In the figure, the red dotted lines represent the undeformed frame and the blue solid lines represent the deformed frame. The light blue element is the element under the highest tensile stress, and the green element is the element under the highest compressive stress. Node 33 is identified, as it is the node under the highest rotational displacement. Element 55 in pink is identified as the element with one of the highest shear stresses. Element 77 is identified by the arrow, as it is element under the highest tensile stress, as well as the highest shear stress.

Notes: The arrow and node 33 identifier in the figure were added manually, and element 77 will be pink when running the MATLAB code above. These things were changed manually to better convey the message of the problem.

It should also be noted that this figure is an "incomplete" figure. A complete figure would show the elements as curved lines, showing the bending the elements occur. This would have been cumbersome to show in the figure as there are so many elements, and the coding for this would have been considerably more difficult.

Frame System vs. Truss System
Compared to the same load in a frame system, as seen in the Original Electric Pylon Problem one notices from the figures that the frame system as a whole deforms less than the truss system. This is because the frame system accounts for bending of the elements and rotations at the nodes. Also it is interesting to note that different elements are now under the highest stresses. In the frame system, the element with the highest compressive stress is now located farther from the actual applied load.

Also it should be noted that the stresses and forces in the elements is lower in the frame system. This is because the forces are now distributed as axial loads and as bending moments. There are also no shear forces in the truss problem, as there are no bending moments in the truss.