User:Eml4507.s13.team3.steiner/Team Negative Damping (3): Report 3

= Problem 1: Pb. 3.1, p. 3-8 Section 3 Notes using Matlab and CALFEM = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Plane truss as shown. A force, F, of 20000 N is applied as shown. The length of each element is 1 meter. Element 2 is vertical and the elements are spaced evenly around node 1. Each element has the same physical properties: Young's modulus of 206 GPa and cross sectional area of .0001 m. Additionally, node 2 experiences an initial displacement of 2 cm in the x direction and -1 cm in the y direction; node 3 experiences an initial displacement of -3 cm in the x direction and 5 cm in the y direction.

Find
Find the unknown displacements and forces at each global node as well as the member forces of each element using Matlab and verify using CALFEM. Plot the deformed and undeformed shape on the same plot.

Using Matlab
Enter the number of global nodes. >> G=10;

Enter the x, y, and z coordinates. >> x=[-1.5 1.5 -1.5 1.5 1.5 -1.5 -4 4 4 -4]; >> y=[0 0 1.5 1.5 -1.5 -1.5 4 4 -4 -4]; >> z=[8 8 4 4 4 4 0 0 0 0];

Enter the input force matrix. The first element of each row corresponds to the global degree of freedom; the second element corresponds to the force associated with that degree of freedom. Only the known forces are entered into this matrix. >> F_in=[1 0; 2 60000; 3 0; 4 0; 5 60000; 6 0; 7 0; 8 0; 9 0; 10 0; 11 0; 12 0; 13 0; 14 0; 15 0; 16 0; 17 0; 18 0];

Enter the input displacement matrix. The first element of each row corresponds to the global degree of freedom; the second element corresponds to the displacement associated with that degree of freedom. >> Q_in=[19 0; 20 0; 21 0; 22 0; 23 0; 24 0; 25 0; 26 0; 27 0; 28 0; 29 0; 30 0];

Enter the number of elements. >> E=25;

Enter the global node corresponding to the first local node of each element. The column number corresponds to the element number of the given node. >> l1=[1 1 2 1 2 2 2 1 1 3 4 3 5 3 6 4 5 4 3 5 6 6 3 4 5];

Enter the global node corresponding to the second local node of each element. >> l2=[2 4 3 5 6 4 5 3 6 6 5 4 6 10 7 9 8 7 8 10 9 10 7 8 9];

Enter the Young's Modulus associated with each element. The column number corresponds to the element number of the given node. >> Modulus=3e7*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];

Enter the Area associated with each element. The column number corresponds to the element number of the given node. >> Area=(pi/4)*2^2*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];

Run the displacement program. >> [k Q F Q_bar F_bar]=displacement(G, x, y, z, F_in, Q_in, E, l1, l2, Modulus, Area);

The stiffness matrix k: k = 1.0e+07 * 3.0900        0         0         0   -1.5450    0.8920   -1.5450   -0.8920            0    3.0900         0   -2.0600    0.8920   -0.5150   -0.8920   -0.5150            0         0         0         0         0         0         0         0            0   -2.0600         0    2.0600         0         0         0         0      -1.5450    0.8920         0         0    1.5450   -0.8920         0         0       0.8920   -0.5150         0         0   -0.8920    0.5150         0         0      -1.5450   -0.8920         0         0         0         0    1.5450    0.8920      -0.8920   -0.5150         0         0         0         0    0.8920    0.5150

The image to the right shows a plot of the deformed and undeformed truss. The displacement matrix Q: Q = 1.0000  -0.0290       2.0000    0.0108       3.0000         0       4.0000    0.0200       5.0000   -0.0100       6.0000         0       7.0000   -0.0300       8.0000    0.0500       9.0000         0      10.0000         0      11.0000         0      12.0000         0

This shows that the first and second degrees of freedom (x and y directions) for global node 1 experience a displacement of -2.9 cm and 1.08 cm respectively. The remaining displacements were given in the problem.

The Force matrix, F F = 1.0e+05 * 0.1414      0.1414            0            0      -4.2816            0      -3.6562       2.1109            0       3.5148       2.0293            0

This shows the forces at each node. Every 3 rows correspond to 1 node. Global node 1 (rows 1-3) shows the force that was given in the problem.

The internal spring force matrix F_bar F_bar = 1.0e+05 * 4.2219   4.2816    4.0586      -4.2219   -4.2816   -4.0586

Row 1 gives the member forces with respect to the Global nodes and row 2 gives the member forces with respect to the members.

Matlab Code
This code is equivalent to that used in R2.3 except that it uses the given spring constants rather than than the modulus and area to determine the stiffness matrix.

CALFEM Verification
First begin by constructing the Edof  matrix. Column 1 corresponds to the element number; columns 2 and 3 correspond to the global node that corresponds to the first and second local nodes respectively. >> Edof=[1 1 2 5 6 2 1 2 3 4           3 1 2 7 8];

Construct an empty stiffness matrix K that has side lengths equal to the number of degrees of freedom times the number of global nodes. For this problem, there are 4 global nodes and motion is constrained to 2 degrees of freedom. >> K=zeros(8,8);

Construct the force matrix F. The forces corresponding to global node 1 are known; the forces corresponding to global nodes 2, 3, and 4 are unknown. >> F=zeros(8,1); >> F(1)=20000/sqrt(2); >> F(2)=20000/sqrt(2);

Construct the displacement matrix Q. Global nodes 2, 3, and 4 and known displacements as shown. >> Q=[3 .02 4 -.01        5 -.03         6 .05         7 0         8 0];

The function bar2e generates the element stiffness matrices for a bar with 2 degrees of freedom. Its parameters are the x and y coordinates of the element and a vector containing the Young's modulus and cross sectional area. This information is given in the problem statement.

The Vector containing the Young's modulus and cross sectional area. >> E=206e9; >> A=1e-4; >> ep=[E A];

The x and y coordinates for each element. >> ex1=[0 sqrt(3)/2]; >> ey1=[0 -.5]; >> ex2=[0 0]; >> ey2=[0 1]; >> ex3=[0 -sqrt(3)/2]; >> ey3=[0 -.5];

The commands: >> Ke1=bar2e(ex1, ey1, ep); >> Ke2=bar2e(ex2, ey2, ep); >> Ke3=bar2e(ex3, ey3, ep);

Yield the following stiffness matrices: >> Ke1 = 1.0e+07 * 1.5450  -0.8920   -1.5450    0.8920          -0.8920    0.5150    0.8920   -0.5150          -1.5450    0.8920    1.5450   -0.8920           0.8920   -0.5150   -0.8920    0.5150      Ke2 = 0          0           0           0                0    20600000           0   -20600000                0           0           0           0                0   -20600000           0    20600000   >> Ke3 = 1.0e+07 * 1.5450   0.8920   -1.5450   -0.8920          0.8920    0.5150   -0.8920   -0.5150         -1.5450   -0.8920    1.5450    0.8920         -0.8920   -0.5150    0.8920    0.5150

The global stiffness matrix can be assembled with the element stiffness matrices using the following command. >> K=assem(Edof(1,:),K,Ke1); >> K=assem(Edof(2,:),K,Ke2); >> K=assem(Edof(3,:),K,Ke3);

This results in the following stiffness matrix. >> K = 1.0e+07 * 3.0900        0         0         0   -1.5450    0.8920   -1.5450   -0.8920               0    3.0900         0   -2.0600    0.8920   -0.5150   -0.8920   -0.5150               0         0         0         0         0         0         0         0               0   -2.0600         0    2.0600         0         0         0         0         -1.5450    0.8920         0         0    1.5450   -0.8920         0         0          0.8920   -0.5150         0         0   -0.8920    0.5150         0         0         -1.5450   -0.8920         0         0         0         0    1.5450    0.8920         -0.8920   -0.5150         0         0         0         0    0.8920    0.5150

This stiffness matrix is in accord with that computed in Matlab in the previous section.

The function solveq takes the stiffness matrix, known force matrix, and displacement matrix and returns the complete force and displacement matrices. >> [q r]=solveq(K,F,Q) q = -0.0290      0.0108       0.0200      -0.0100      -0.0300       0.0500            0            0   r = 1.0e+05 * 0.0000      0.0000            0      -4.2816      -3.6562       2.1109       3.5148       2.0293

The displacement and force matrices are in accord with those computed in Matlab in the previous section.

The function extract is used to evaluate the displacements of the local nodes for each beam, which is needed to determine the force in each beam. >> ed1=extract(Edof(1,:),q); >> ed2=extract(Edof(2,:),q); >> ed3=extract(Edof(3,:),q);

This results in the following displacements. The number following 'ed' corresponds to the element number. Each column corresponds to the respective degree of freedom for the beam. ed1 = -0.0290   0.0108   -0.0300    0.0500 ed2 = -0.0290   0.0108    0.0200   -0.0100 ed3 = -0.0290   0.0108         0         0

The function bar2s can be used to determine the force in each beam. Its input parameters are the x and y coordinates for the beam (ex1 and ey1), the physical properties of the beam (Young's modulus and cross sectional area - ep) and the displacement (ed). >> es1=bar2s(ex1, ey1,ep,ed1) >> es1=bar2s(ex2, ey2,ep,ed2) >> es1=bar2s(ex3, ey3,ep,ed3)

This results in the following spring forces. es1 = -4.2219e+05 es2 = -4.2816e+05 es3 = -4.0586e+05

This is in accord with the forces computed in Matlab in the previous section.

= Problem 2: Dampened and Undamped System= On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Spring-damper-body arrangement as shown. Two separate forces applied to masses.

$$ M= \begin{bmatrix} m_{1} & 0\\ 0 & m_{2}\\ \end{bmatrix} $$ $$ d= \begin{bmatrix} d_{1}\\ d_{2}\\ \end{bmatrix} $$ $$ C= \begin{bmatrix} C_{1}+C_{2} & -C_{2}\\ -C_{2} & C_{2}+C_{3}\\ \end{bmatrix}$$ $$ K= \begin{bmatrix} (k_{1}+k_{2}) & -k_{2}\\ -k_{2} & (k_{2}+k_{3})\\ \end{bmatrix}

$$

Find
1.Draw the FBDs of all components in the spring-mass-damper system on p.53-13, with the known disp dofs being at the left and right supports:

$$ d_{3}=0 $$

$$  d_{4}=0 $$

2. Derive (1) p.53-12, and (1)-(3) p. 53-13.

Free Body Diagrams (Question 1)
The free body diagrams for a dampened system and undamped system can be seen to the right.

DAMPENED
Analyzing the forces on mass one, we obtain:
 * {| style="width:100%" border="0"

$$\displaystyle m_{1}d''_{1}+ k_{1}d_{1} +C_{1}d'_{1} - k_{2}d_{2}+ k_{2}d_{1}-C_{2}d'_{2}+C_{2}d'_{1} = F_{1}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


 * {| style="width:100%" border="0"

$$\displaystyle m_{1}d''_{1} + C_{1}d'_{1} + C_{2}(d'_{1}- d'_{2}) +k_{1}d_{1}+ k_{2}(d_{1}- d_{2})= F_{1}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Analyzing the forces on mass two, we obtain:
 * {| style="width:100%" border="0"

$$\displaystyle m_{2}d''_{2} + C_{2}(d'_{2}- d'_{1}) + k_{2}(d_{2}- d_{1}) +C_{3}d'_{2} +k_{3}d_{2} = F_{2}$$
 * style="width:97%" |
 * style="width:97%" |
 * }
 * {| style="width:100%" border="0"

$$\displaystyle m_{2}d''_{2} + C_{2}(d'_{2}- d'_{1}) +C_{3}d'_{2}+ k_{2}(d_{2}- d_{1}) +k_{3}d_{2} = F_{2}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Plugging equations into matrix form we obtain:

$$

\begin{bmatrix} m_{1} & 0\\ 0 & m_{2}\\ \end{bmatrix} \begin{bmatrix} d''_{1}\\ d''_{2}\\ \end{bmatrix}+ \begin{bmatrix} C_{1}+C_{2} & -C_{2}\\ -C_{2} & C_{2}+C_{3}\\ \end{bmatrix} \begin{bmatrix} (k_{1}+k_{2}) & -k_{2}\\ -k_{2} & (k_{2}+k_{3})\\ \end{bmatrix} \begin{bmatrix} d'_{1}\\ d'_{2}\\ \end{bmatrix}+ \begin{bmatrix} d_{1}\\ d_{2}\\ \end{bmatrix} $$\

Notice how this form matches each matrix given in the problem statement. Plugging in the original yields and proves the following equation of motion for a damped system:

$$ Md''+ Cd'+kd=0 $$

UNDAMPED
Analyzing the forces on mass one, we obtain:
 * {| style="width:100%" border="0"

$$\displaystyle m_{1}d''_{1}+ k_{1}d_{1} - k_{2}d_{2}+ k_{2}d_{1} = F_{1}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


 * {| style="width:100%" border="0"

$$\displaystyle m_{1}d''_{1} - k_{2}d_{2} +(k_{1}+k_{2})d_{1} = F_{1}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Analyzing the forces on mass two, we obtain:


 * {| style="width:100%" border="0"

$$\displaystyle m_{2}d''_{2}+ k_{2}d_{2} - k_{2}d_{1}+ k_{3}d_{2} = F_{2}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


 * {| style="width:100%" border="0"

$$\displaystyle m_{2}d''_{2}+ (k_{2}+k_{3})d_{2} - k_{2}d_{1}= F_{2}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Plugging equations into matrix form we obtain:

$$

\begin{bmatrix} m_{1} & 0\\ 0 & m_{2}\\ \end{bmatrix} \begin{bmatrix} d''_{1}\\ d''_{2}\\ \end{bmatrix}+ \begin{bmatrix} (k_{1}+k_{2}) & -k_{2}\\ -k_{2} & (k_{2}+k_{3})\\ \end{bmatrix} \begin{bmatrix} d_{1}\\ d_{2}\\ \end{bmatrix} $$\

Notice how this form matches each matrix given in the problem statement. Plugging in the original yields and proves the following equation of motion for an undamped system:

$$ Md''+kd=0 $$

= Problem 3: Pb. 53.6, sec53b =

Given the following values, use them to integrate the governing system of L2-ODEs-CC.

$$ m_1=3 \ m_2=2 $$ $$ c_1=1/2 \ c_2=1/4 \ c_3=1/3 $$ $$ k_1=10 \ k_2=20 \ k_3=15$$ $$ F_1(t)=0 \ F_2(t)=0 $$ $$ d_1(0)=-1 \ d_2(0)=2 $$ $$ d'_1(0)=0 \ d'_2(0)=0 $$

the engin command helps to find the time step size by finding the enginvalues of the equation λ=L $$ Kx=LMx $$

Then the smaller of the two enginvalues found

$$ 0<L_1<L_2 $$

Which is equivalent to the frequency squared

$$ L_1=(1/T_1)^2 $$

Then using $$ dt=T_1/10 $$ Using the time interval with "step 2" to integrate the governing L2-ODEs-CC $$ 0<t<T=5T_1 $$

Matlab shown below

.... Finding the eigenvales of the L2-ODEs-CC... L=eigen(K,M); [L,X]=eigen(K,M); T=L(1,1:)^2
 * K-LM| = 0

....Using step 2 to integrate the governing system of L2-ODEs_CC... Dsnap=step2(K,C,M,d0,v0,ip,f,pdisp); [Dsnap,D,V,A]=step2(K,C,M,d0,v0,ip,f,pdisp); ip = [dt T α δ [nsnap nhist timei dofi]] α=1/4 δ=1/2 dt=T/10; Plot(ip),d_1(t),d_2(t)

= Problem 4: Pb. 10, p. 100 Kim and Sankar using Matlab and CALFEM = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
The stepped bar shown in the figure above is subjected to a force at the center. Assume: E = 100 GPa, area of cross sections of the three portions shown are, respectively, 10-4m2, 2x10-4m2, and 10-4m2, and F = 10,000 N.

Find
Use FEM to determine the displacement at the center and reactions RL and RR.

Matlab
Using the Matlab Code from Problem 1 and the following inputs the displacement at the center and reactions RL and RR were determined.

The Matlab Code first takes the number of global nodes and their x, y, and z coordinates. Below is that information formatted for use in the code. This information is loaded as variables in the Matlab work space. G=5; x=[0 0.3 0.5 0.7 1]; y=[0 0 0 0 0]; z=[0 0 0 0 0];

The next parameters the code needs are the number of elements in the structure along with their corresponding local nodes. These variables are also stored in the Matlab work space. E=4; l1=[1 2 3 4]; l2=[2 3 4 5];

This is the force matrix used in the Matlab Code. The given forces are the known forces in the structure. The code does not need the known forces associated with the y and z coordinates because this is a 1-D problem. F_in=[4 0; 7 10000; 10 0];

The next parameter the Matlab Code needs is the known initial displacements. This array does include the initial y and z corresponding displacements. Not listed are the unknowns in the x direction at global node 2, 3, and 4. Q_in=[1 0; 2 0; 3 0; 5 0; 6 0; 8 0; 9 0; 11 0;... ...12 0; 13 0; 14 0; 15 0];

The final information needed by the code are the Young's Modulus and Area of each element. Modulus=[1e11 1e11 1e11 1e11]; Area=[1e-4 2e-4 2e-4 1e-4];

Running the code results in the top right figure which shows the initial structure as a blue dashed line and the displaced structure as a solid green line. The displacement in this structure is very small, an enlarged figure of global node 2 is shown in the bottom right figure to verify that the nodes have displaced.

$$Q_{x3}=0.2 m$$ $$R_{L}=5000 N$$ $$R_{R}=5000 N$$

CALFEM
Since the code for this problem was the same at the code in Problem 1 then the CALFEM Verification confirms that the code is working properly.

=Problem 3.5: 2 Finding Displacement and Forces in 2 Element Truss System= On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
A truss has two elements, as shown in the figure. The members are made of aluminum hollow square cross-section. The outer dimension of the square is 12 mm, and the inner dimension is 9 mm. (The wall thickness is 1.5 mm in all four sides.) Assume Young's Modulus E= 70 GPa and yield strenght $$\sigma _{Y}$$  = 70 MPa. The magnitude of the force at Global Node 3 (F) is equal to 1,000 N.

Find
Use FEM to determine the displacement at Global Node 3 and axial forces in Elements 1 and 2. Use von Mises' yield theory to determine if the elements will yield or not. Use Euler buckling load to determine if the elements under compressive loads will buckle.

Redesign the truss so that both the stress and buckling constraints are satisfied with a Safety Factor of N not less than 2 for stresses and N not less than 1.2 for buckling. The design should reduce the weight as much as possible.

Draw the truss designed and provide the nodal coordinates and element connectivity in the form of a table. Include the nodal displacements, forces in each member, and safety factors for stresses and buckling for each element in the form of a table.

Solution
The Matlab Code from Problem 1 is used and the specific inputs for this problem were used.

The first thing that is inputed in the Matlab code is the total number of global nodes. Next, the position of the nodes is input by putting the value of the x, y, and z coordinates for each node. G=3; x=[0 0 2]; y=[0 1 0]; z=[0 0 0];

Next the number of elements in the system is input. Next we put which global node corresponds to each local node is input for l1 and l2. E=2; l1=[1 2]; l2=[3 3];

The known forces are then input. There is only one force so we input the components of the force at the degrees of freedom 7 and 8, which is at global node 3. F_in=[7 0; 8 -1000];

Next we input the known displacements. These displacements include the forces at the reactions which are the only known forces. Q_in=[1 0; 2 0; 3 0; 4 0; 5 0; 6 0; 9 0];

Lastly, we input the Young's modulus and the cross-sectional area. The members have hollow square cross-sections. Modulus=[7e10 7e10 7e10]; Area=[6.3e-5 6.3e-5 6.3e-5];

The following image represents the initial system with a blue dashed line and shows the displaced system with a solid green line.

The displacement at Global Node 3 is: Q_out = 7.0000  -0.0009    8.0000   -0.0043

The axial forces in Elements 1 and 2 are: F = 0        0   -0.0009    0.0011

The von Mises theory states that a material will yield when the von Mises stress exceeds the yield stress for the material. $$\sigma _{VM}\geqslant \sigma _{Y}$$

The von Mises stress in the Elements of the truss is found to be 15 MPa and the yield stress is given to be 70 MPa. This means that $$\sigma _{VM}< \sigma _{Y}$$  so the elements will NOT yield. The Euler buckling load is equal to $$P_{cr}=\pi ^{2}EI/L^{^{2}}$$ where the moment of inertia of the cross section is given by $$I=\frac{a_{0}^{4}-a_{i}^{4}}{12}$$

For element 1, the Euler buckling load is 0.8161 N. For element 2, the Euler buckling is 0.1632 N.

Element 1: $$P_{cr}=\frac{\pi ^{2}(70 GPa)\frac{(12 mm)^{4}-(9 mm)^{4}}{12}}{1 m^{2}}=0.8161 N$$

Element 2: $$P_{cr}=\frac{\pi ^{2}(70 GPa)\frac{(12 mm)^{4}-(9 mm)^{4}}{12}}{\sqrt{5 m}^{2}}=0.1632 N$$

For Elements 1 and 2: $$F< P_{cr}$$  so the elements will NOT buckle under compressive force.

After conducting a brief analysis of different truss designs, the best one was found to be the same design as provided in the problem. The different designs that were looked at were a rectangular truss system with vertical, horizontal, and diagonal elements, an isoceles trianglular truss system with the elements being connected to the wall at the same angle, and the square triangular truss system from the problem. The right triangular system was chosen because it has the best mass to stress handling ratio. The figure shows the truss system.

Once the design was chosen, different dimensions were tested. Here is a list of the tests and the results. The original weight of the truss was 74.72 kg. The new weight of the truss is: $$W= (6.3*10^{-3}m^{2})(1.25 m+1.458 m)(2800\frac{kg}{m^{3}})=47.77 kg$$

The nodal coordinates, element connectivity and nodal displacements are included in the following table.

The displacements and forces were found using the same Matlab program used in the previous section. The displacement matrix is: Q_out = 7.0000  -0.0000    8.0000   -0.0020 The force matrix is: F = 1.0e-003 * 0        0 -0.4724   0.6425

= Problem 6: Project 2.3, p. 97 Kim and Sankar using Matlab and CALFEM = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Consider the plane truss in the figure above. The horizontal and vertical members have length L, while inclined members have length 1.414*L. Assume the Young's modulus E = 100 GPa, cross-sectional area A = 1.0 cm2, and L = 0.3 m.

Find
1. Determine the deflections and element forces for the following three load cases.

-Load Case A: Fx13 = Fx14 = 10,000 N

-Load Case B: Fy13 = Fy14 = 10,000 N

-Load Case C: Fx13 = 10,000 N Fx14 = -10,000 N

2. Find the three beam properties, axial rigidity(EA)eq, flexural rigidity (EI)eq, and shear rigidity (GA)eq.

3. Verify the beam model by adding two more bays to the truss (l = 8*0.3 = 2.4 m). Compute the tip deflections of the extended truss for the three load cases A-C using the Matlab program and compare the results with the deflections obtained from the equivalent beam model.

Using Matlab
1. For Load Case A Using the Matlab Code from Problem 1 and the following inputs the displacement at the center and reactions RL and RR were determined.

The Matlab Code first takes the number of global nodes and their x, y, and z coordinates. Below is that information formatted for use in the code. This information is loaded as variables in the Matlab work space. G=14; x=[0 0 0.3 0.3 0.6 0.6 0.9 0.9 1.2 1.2 1.5 1.5 1.8 1.8]; y=[0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3]; z=[0 0 0 0 0 0 0 0 0 0 0 0 0 0];

The next parameters the code needs are the number of elements in the structure along with their corresponding local nodes. These variables are also stored in the Matlab work space. E=25; l1=[1 3 5 7 9 11 2 4 6  8 10 12 1 3 5 7  9 11 13 1 3...      ...5  7  9 11]; l2=[3 5 7 9 11 13 4 6 8 10 12 14 2 4 6 8 10 12 14 4 6...    ...8 10 12 14];

This is the force matrix used in the Matlab Code. The given forces are the known forces in the structure. F_in=[2 0; 7 0; 8 0; 10 0; 11 0; 13 0; 14 0; 16 0; 17 0; 19 0; 20 0; 22 0; 23 0; 25 0; 26 0; 28 0; 29 0; 31 0; 32 0; 34 0;...    ...35 0; 37 10000; 38 0; 40 10000; 41 0];

The next parameter the Matlab Code needs is the known initial displacements. Q_in=[1 0; 3 0; 4 0; 5 0; 6 0; 9 0; 12 0; 15 0; 18 0; 21 0; 24 0; 27 0; 30 0; 33 0; 36 0; 39 0; 42 0];

The final information needed by the code are the Young's Modulus and Area of each element. Modulus=[1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11... ...1e11 1e11]; Area=[1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4... ...1e-4 1e-4];

Running the code results in the top right figure which shows the initial structure as a blue dashed line and the displaced structure as a solid green line.

For Load Case B Using the Matlab Code from Problem 1 and the following inputs the displacement at the center and reactions RL and RR were determined.

The Matlab Code first takes the number of global nodes and their x, y, and z coordinates. Below is that information formatted for use in the code. This information is loaded as variables in the Matlab work space. G=14; x=[0 0 0.3 0.3 0.6 0.6 0.9 0.9 1.2 1.2 1.5 1.5 1.8 1.8]; y=[0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3]; z=[0 0 0 0 0 0 0 0 0 0 0 0 0 0];

The next parameters the code needs are the number of elements in the structure along with their corresponding local nodes. These variables are also stored in the Matlab work space. E=25; l1=[1 3 5 7 9 11 2 4 6  8 10 12 1 3 5 7  9 11 13 1 3...      ...5  7  9 11]; l2=[3 5 7 9 11 13 4 6 8 10 12 14 2 4 6 8 10 12 14 4 6...    ...8 10 12 14];

This is the force matrix used in the Matlab Code. The given forces are the known forces in the structure. F_in=[2 0; 7 0; 8 0; 10 0; 11 0; 13 0; 14 0; 16 0; 17 0; 19 0; 20 0; 22 0; 23 0; 25 0; 26 0; 28 0; 29 0; 31 0; 32 0; 34 0;...   ...35 0; 37 0; 38 10000; 40 0; 41 10000];

The next parameter the Matlab Code needs is the known initial displacements. Q_in=[1 0; 3 0; 4 0; 5 0; 6 0; 9 0; 12 0; 15 0; 18 0; 21 0; 24 0; 27 0; 30 0; 33 0; 36 0; 39 0; 42 0];

The final information needed by the code are the Young's Modulus and Area of each element. Modulus=[1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11... ...1e11 1e11]; Area=[1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4... ...1e-4 1e-4];

Running the code results in the top right figure which shows the initial structure as a blue dashed line and the displaced structure as a solid green line.

For Load Case C Using the Matlab Code from Problem 1 and the following inputs the displacement at the center and reactions RL and RR were determined.

The Matlab Code first takes the number of global nodes and their x, y, and z coordinates. Below is that information formatted for use in the code. This information is loaded as variables in the Matlab work space. G=14; x=[0 0 0.3 0.3 0.6 0.6 0.9 0.9 1.2 1.2 1.5 1.5 1.8 1.8]; y=[0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3 0 0.3]; z=[0 0 0 0 0 0 0 0 0 0 0 0 0 0];

The next parameters the code needs are the number of elements in the structure along with their corresponding local nodes. These variables are also stored in the Matlab work space. E=25; l1=[1 3 5 7 9 11 2 4 6  8 10 12 1 3 5 7  9 11 13 1 3...      ...5  7  9 11]; l2=[3 5 7 9 11 13 4 6 8 10 12 14 2 4 6 8 10 12 14 4 6...    ...8 10 12 14];

This is the force matrix used in the Matlab Code. The given forces are the known forces in the structure. F_in=[2 0; 7 0; 8 0; 10 0; 11 0; 13 0; 14 0; 16 0; 17 0; 19 0; 20 0; 22 0; 23 0; 25 0; 26 0; 28 0; 29 0; 31 0; 32 0; 34 0;...   ...35 0; 37 10000; 38 0; 40 -10000; 41 0];

The next parameter the Matlab Code needs is the known initial displacements. Q_in=[1 0; 3 0; 4 0; 5 0; 6 0; 9 0; 12 0; 15 0; 18 0; 21 0; 24 0; 27 0; 30 0; 33 0; 36 0; 39 0; 42 0];

The final information needed by the code are the Young's Modulus and Area of each element. Modulus=[1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11 1e11... ...1e11 1e11]; Area=[1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4... ...1e-4 1e-4];

Running the code results in the top right figure which shows the initial structure as a blue dashed line and the displaced structure as a solid green line.

Displacements For the Three Cases in millimeters

2.

Using these equations and the displacements at nodes 13 and 14 the axial rigidity, flexural rigidity, and shear rigidity for the equivalent beam was found.

$$u_{tip}=\frac{Fl}{(EA)_{eq}}$$

$$v_{tip}=\frac{Fl^{3}}{3(EI)_{eq}}+\frac{Fl}{(GA)_{eq}}$$

$$v_{tip}=\frac{Cl^{2}}{2(EI)_{eq}}$$

$$(EA)_{eq}=20 MN$$

$$(EI)_{eq}=574.6 MNm$$

$$(EA)_{eq}=1.064 MN$$

3.

From the Matlab code the beam deflection was $$u_{tip}=0.0022$$ and $$v_{tip}=0.1021$$ and from the parameters found in 2. the beam delection was $$u_{tip}=0.0024$$ and $$v_{tip}=0.1023$$. These values verify the use of the equations found in 2.

CALFEM Verification
Since the code for this problem was the same at the code in Problem 1 then the CALFEM Verification confirms that the code is working properly.

=Problem 7 (fead.f08 p.14-3, HW)= On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Bar element stiffness matrix: $$\mathbf k^{(e)} = k^{(e)} \begin{bmatrix} (l^{(e)})^2 & l^{(e)}m^{(e)} & -(l^{(e)})^2 & -l^{(e)}m^{(e)} \\l^{(e)}m^{(e)} & (m^{(e)})^2 & -l^{(e)}m^{(e)} & -(m^{(e)})^2 \\-(l^{(e)})^2 & -l^{(e)}m^{(e)} & (l^{(e)})^2 & l^{(e)}m^{(e)} \\-l^{(e)}m^{(e)} & -(m^{(e)})^2 & l^{(e)}m^{(e)} & (m^{(e)})^2 \end{bmatrix}$$,

Where $$k^{(e)} = \frac{EA}{L}$$, $$l^{(e)} = \cos \theta^{(e)}$$, and $$m^{(e)} = \sin \theta^{(e)}$$

Find

 * Verify that the bar element stiffness matrix $$\mathbf k^{(e)}$$ in global coordinates can be obtained by $$\mathbf k^{(e)} = {\mathbf T}^{{(e)}T} \, \hat{\mathbf k}^{(e)} \, {\mathbf T}^{(e)}$$.
 * Verify that you also get the same expression for the bar element stiffness matrix in global coordinates with $$\mathbf k^{(e)} = \widetilde{\mathbf T}^{{(e)}T} \, \widetilde{\mathbf k}^{(e)} \, \widetilde{\mathbf T}^{(e)}$$, where $$\widetilde{\mathbf T}^{(e)}$$ is the 4x4 transformation matrix in K.2008.2 p.32, and $$\widetilde{\mathbf k}^{(e)}$$ is the 4x4 bar element stiffness matrix in local coordinates in K.2008.2 p.31.

Given
$$\mathbf \hat k^{(e)} = k^{(e)} \begin{bmatrix} 1 & -l \\-l & 1 \end{bmatrix}$$, $$\mathbf T^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}$$, $$\mathbf T^{(e)T} = \begin{bmatrix} l^{(e)} & 0 \\m^{(e)} & 0 \\0 & l^{(e)} \\0 & m^{(e)} \end{bmatrix}$$

Check that given matrix dimensions correspond
$$\mathbf k^{(e)} = {\mathbf T}^{{(e)}T} \, \hat{\mathbf k}^{(e)} \, {\mathbf T}^{(e)} \Rightarrow \begin{bmatrix} 4 \times 4 \end{bmatrix} = \begin{bmatrix} 4 \times 2 \end{bmatrix} \begin{bmatrix} 2 \times 2 \end{bmatrix} \begin{bmatrix} 2 \times 4 \end{bmatrix} $$

Multiply matrices
$$\mathbf T^{(e)T} \hat{\mathbf k}^{(e)} = k^{(e)} \begin{bmatrix} l^{(e)} & 0 \\m^{(e)} & 0 \\0 & l^{(e)} \\0 & m^{(e)} \end{bmatrix} \begin{bmatrix} 1 & -l \\-l & 1 \end{bmatrix}$$


 * $$= k^{(e)} \begin{bmatrix} l^{(e)} & -l^{(e)} \\m^{(e)} & -m^{(e)} \\-l^{(e)} & l^{(e)} \\-m^{(e)} & m^{(e)} \end{bmatrix}$$, where $$k^{(e)}$$ can be moved in front because it is a scalar.

$${\mathbf T}^{{(e)}T} \, \hat{\mathbf k}^{(e)} \, {\mathbf T}^{(e)} = k^{(e)} \begin{bmatrix} l^{(e)} & -l^{(e)} \\m^{(e)} & -m^{(e)} \\-l^{(e)} & l^{(e)} \\-m^{(e)} & m^{(e)} \end{bmatrix} \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix} $$


 * $$= k^{(e)} \begin{bmatrix} (l^{(e)})^2 & l^{(e)}m^{(e)} & -(l^{(e)})^2 & -l^{(e)}m^{(e)} \\l^{(e)}m^{(e)} & (m^{(e)})^2 & -l^{(e)}m^{(e)} & -(m^{(e)})^2 \\-(l^{(e)})^2 & -l^{(e)}m^{(e)} & (l^{(e)})^2 & l^{(e)}m^{(e)} \\-l^{(e)}m^{(e)} & -(m^{(e)})^2 & l^{(e)}m^{(e)} & (m^{(e)})^2 \end{bmatrix} = \mathbf k^{(e)}$$

Given
4x4 bar element stiffness matrix in local coordinates: $$\widetilde{\mathbf k}^{(e)} = k^{(e)} \begin{bmatrix} 1 & 0 & -1 & 0 \\0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$

Transformation matrix: $$\widetilde{\mathbf T}^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\-m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \\ 0 & 0 & -m^{(e)} & l^{(e)} \end{bmatrix}$$

Transposed transformation matrix: $$\widetilde{\mathbf T}^{(e)T} = \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix}$$

Check that given matrix dimensions correspond
$$\mathbf k^{(e)} = \widetilde{\mathbf T}^{{(e)}T} \, \widetilde{\mathbf k}^{(e)} \, \widetilde{\mathbf T}^{(e)} \Rightarrow \begin{bmatrix} 4 \times 4 \end{bmatrix} = \begin{bmatrix} 4 \times 4 \end{bmatrix} \begin{bmatrix} 4 \times 4 \end{bmatrix} \begin{bmatrix} 4 \times 4 \end{bmatrix} $$

Multiply matrices
$$\widetilde{\mathbf T}^{(e)T} \widetilde{\mathbf k}^{(e)} = k^{(e)} \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix} \begin{bmatrix} 1 & 0 & -1 & 0 \\0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} $$


 * $$= k^{(e)} \begin{bmatrix} l^{(e)} & 0 & -l^{(e)} & 0 \\m^{(e)} & 0 & -m^{(e)} & 0 \\ -l^{(e)} & 0 & l^{(e)} & 0 \\ -m^{(e)} & 0 & m^{(e)} & 0 \end{bmatrix}$$

$$\widetilde{\mathbf T}^{{(e)}T} \, \widetilde{\mathbf k}^{(e)} \, \widetilde{\mathbf T}^{(e)} = k^{(e)} \begin{bmatrix} l^{(e)} & 0 & -l^{(e)} & 0 \\m^{(e)} & 0 & -m^{(e)} & 0 \\ -l^{(e)} & 0 & l^{(e)} & 0 \\ -m^{(e)} & 0 & m^{(e)} & 0 \end{bmatrix} \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\-m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \\ 0 & 0 & -m^{(e)} & l^{(e)} \end{bmatrix} $$


 * $$= k^{(e)} \begin{bmatrix} (l^{(e)})^2 & l^{(e)}m^{(e)} & -(l^{(e)})^2 & -l^{(e)}m^{(e)} \\l^{(e)}m^{(e)} & (m^{(e)})^2 & -l^{(e)}m^{(e)} & -(m^{(e)})^2 \\-(l^{(e)})^2 & -l^{(e)}m^{(e)} & (l^{(e)})^2 & l^{(e)}m^{(e)} \\-l^{(e)}m^{(e)} & -(m^{(e)})^2 & l^{(e)}m^{(e)} & (m^{(e)})^2 \end{bmatrix} = \mathbf k^{(e)}$$

=Problem 8 (sec.53a, p.53-2b, Pb-53.4)= On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given

 * Consider an L2-ODE-CC with the following roots: $$\lambda_1,_2 = -0.5+i2$$. Find the system natural frequency $$\omega$$.


 * Let the system be excited by a periodic force of $$f(t) = f_0cos(\bar \omega t)$$. Find the expression for the particular solution $$y_p(t)$$.


 * Consider:
 * $$f_0/m = r_0 = 1$$


 * $$\bar \omega = 0.9\omega$$


 * Initial Conditions
 * $$y(0) = 1$$
 * $$y'(0) = 0$$

Find

 * Find the particular solution $$y_p(t)$$ and the amplification factor $$\mathbb A$$ for the given input.
 * Find the complete solution $$y(t)$$ satisfying the given initial conditions.
 * Plot in separate figures: $$y_h(t), y_p(t), y(t)$$
 * Plot in the same figure: $$y_h(t), y_p(t), y(t)$$

System natural frequency $$\omega$$
From R1.5 (Case 3) :


 * $$ (\lambda-\lambda_1)(\lambda-\lambda_2)$$
 * $$ (\lambda+.5+i2)(\lambda+.5-i2)$$
 * $$ \lambda^2 + .5\lambda - i2\lambda + .5\lambda + .25 + .5i2 +i2\lambda - .5i2 + 4$$
 * $$ \lambda^2 + \lambda + 4.25 = 0$$

The standard L2-ODE-CC is then
 * $$y'' + y' + 4.25y = 0 $$

Relate L2-ODE-CC to standard form of equation of motion:


 * $$y + y' + 4.25y = y + ay' + by$$
 * $$y'' + 2\zeta \omega y' + \omega^2 y$$
 * $$a = 2\zeta \omega$$
 * $$b = \omega^2$$
 * $$\omega = \sqrt{4.25}$$
 * $$\zeta = \frac1{2\sqrt{4.25}}$$

The system natural frequency: $$\omega = \sqrt{4.25}$$

Expression for particular solution $$y_p(t)$$
To find $$y_p(t)$$, use the formula for trial particular solution:


 * $$y_p(t) = de^{i(\bar \omega t - \phi)}$$,

Where $$d$$ is defined as:
 * $$d=\frac{f_0/m}{\sqrt{(\omega^2 - \bar\omega^2)^2 + (2 \zeta \omega \bar\omega)^2}}$$

And $$\phi$$ is defined by:
 * $$\tan\phi=\frac{2\zeta\omega\bar\omega}{\omega^2 - \bar\omega^2}$$

And the damped circular frequency $$\bar \omega$$
 * $$\bar \omega = \omega \sqrt{1 - \zeta^2}$$

Combining the previous equations, the entire expression is written as:
 * $$y_p(t) = \frac{f_0/m}{\sqrt{(\omega^2 - \bar\omega^2)^2 + (2 \zeta \omega \bar\omega)^2}}e^{i(\bar \omega t - \phi)}$$

But since this gives a complex solution, $$e^{i(\bar \omega t - \phi)}$$ is redefined with Euler's formula:
 * $$e^{i(\bar \omega t - \phi)} = \cos(\bar \omega t - \phi) + i \sin(\bar \omega t - \phi)$$

Since $$f(t) = f_0 \cos \bar \omega t$$ is a real forcing function, $$i \sin(\bar \omega t - \phi)$$ can be removed.

So the final expression for the particular solution $$y_p(t)$$ is:
 * $$y_p(t) = \frac{f_0/m}{\sqrt{(\omega^2 - \bar\omega^2)^2 + (2 \zeta \omega \bar\omega)^2}}\cos(\bar \omega t - \phi)$$

Which can be simplified to:
 * $$y_p(t) = d \cos(\bar \omega t - \phi)$$

Particular solution $$y_p(t)$$ for the case and amplification factor $$\mathbb A$$
Consider the case:
 * $$f_0/m = r_0 = 1$$


 * $$\bar \omega = 0.9\omega$$

The amplification factor $$\mathbb A$$ is given as:
 * $$\mathbb A := \frac1{\left[(1 - \rho^2)^2+(2 \zeta \rho)^2 \right]^{1/2}}$$,

Where $$\rho = \frac {\bar \omega}{\omega} = 0.9$$
 * $$\mathbb A = 2.10032$$

The particular solution $$y_p(t)$$ for this case is:
 * $$y_p(t) = d \cos(\bar \omega t - \phi)$$

Where $$d$$ is redefined in terms of $$\mathbb A$$ as:
 * $$d = \mathbb A \frac{f_0}{k}$$, where $$k = bm = \omega^2 m$$


 * $$d = \mathbb A \frac{f_0}{\omega^2 m}$$

Since $$f_0/m = r_0 = 1$$,


 * $$d = \frac{\mathbb A }{\omega^2}$$


 * $$y_p(t) = \frac{\mathbb A }{\omega^2} \cos(\bar \omega t - \phi)$$

Plugging in the values and simplifying results in the particular solution $$y_p(t)$$:
 * $$y_p(t) = 0.494193 \cos(1.8554 t - 1.1603)$$

Complete solution $$y(t)$$ satisfying $$y(0) = 1, y'(0) = 0$$
By superposition: $$y(t)=y_h(t)+y_p(t)$$
 * $$y_h(t)=C_1 e^{(-0.5 t)} sin(2 t)+C_2 e^{(-0.5 t)} cos(2 t)$$


 * $$y_p(t) = 0.494193 \cos(1.8554 t - 1.1603)$$


 * $$y(t)=C_1 e^{(-0.5 t)} sin(2 t)+C_2 e^{(-0.5 t)} cos(2 t)+ 0.494193 \cos(1.8554 t - 1.1603)$$

Using the conditions $$y(0) = 1, y'(0) = 0$$ with the complete solution and the derivative of the complete solution:


 * $$y(t)=C_1 e^{(-0.5 t)} sin(2 t)+C_2 e^{(-0.5 t)} cos(2 t)+ 0.494193 \cos(1.8554 t - 1.1603)$$
 * $$y'(t)=e^{-0.5t}[(-0.5 C_1 - 0.5 C_2) \sin(2t) + (2 C_1 + 2 C_2) \cos(2t)] - 0.916929 \sin(1.8554t - 1.1603)$$


 * $$y(0)=1=0+C_2+\cos(1.1603)$$
 * $$y(0)=1=0+C_2+0.197213$$
 * $$C_2 = 0.802787$$


 * $$y'(0)=0=-0.916926 \sin(1.1603) + 2 C_1 + 2 C_2$$


 * $$y'(0)=0=0.840752 + 2 C_1 + 1.60557$$


 * $$C_1 = -1.22316$$

The complete solution is:
 * $$y(t)=-1.22316 e^{(-0.5 t)} sin(2 t)+0.802787 e^{(-0.5 t)} cos(2 t)+ 0.494193 \cos(1.8554 t - 1.1603)$$

Plots








=References=