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

=Problem 1: Explanation of Stiffness Matrix Assembly Process= On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Node 2 on [http://upload.wikimedia.org/wikiversity/en/1/13/Eml4500.f08.2.djvu pg. 11-4]

Find
Deduce the assembly process. Then find the equilibrium of node 2 on pg. 11-4.

Solution
Assembly process: starting from 9-1

Below is the local matrix of elements 

$$\left[\begin{array}{cccccc} K^{(1)}_{11} & K^{(1)}_{12} & K^{(1)}_{13} & K^{(1)}_{14} & 0 & 0 \\ K^{(1)}_{21} & K^{(1)}_{22} & K^{(1)}_{23} & K^{(1)}_{24} & 0 & 0 \\ K^{(1)}_{31} & K^{(1)}_{32} & (K^{(1)}_{33} + K^{(2)}_{11} & (K^{(1)}_{34} + K^{(2)}_{12})& K^{(2)}_{13} & K^{(2)}_{14} \\ K^{(1)}_{41} & K^{(1)}_{42} & (K^{(1)}_{43} + K^{(2)}_{21}) & (K^{(1)}_{44} + K^{(2)}_{22}) & K^{(2)}_{23} & K^{(2)}_{24} \\ 0 & 0 & K^{(2)}_{31} & K^{(2)}_{32} & K^{(2)}_{33} & K^{(2)}_{34} \\ 0 & 0 & K^{(2)}_{41} & K^{(2)}_{42} & K^{(2)}_{43} & K^{(2)}_{44} \end{array} \right]\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \end{Bmatrix}=\begin{Bmatrix} F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ F_{5} \\ F_{6} \end{Bmatrix}$$

Using the K values given from class,only three calculations need to be done for each element, the rest can be found from the symmetry of the matrix. They either become positive or negative variations of the values calculated below.  $$ K_{11} = \frac{9}{16} = .5625$$ $$ K_{12} = \frac{3\sqrt{3}}{16} = .3248$$ $$ K_{33} = \frac{49}{16} = 3.0623$$ $$ K_{34} = -2.175$$ $$ K_{43} = -2.175$$ $$ K_{44} = 2.6875$$<br\><br\>

Boundary conditions

Nodes 1 and 2 are fixed therefore no displacement occurs at these nodes

$$\textbf{d} = \begin{Bmatrix} 0 \\ 0 \\ d_{3} \\ d_{4} \\ 0 \\ 0 \end{Bmatrix}$$<br\><br\> $$ d_{1}=d_{2}=d_{5}=d_{6}=0 $$<br\> By eliminating the known degrees of freedom it starts to reduce the global force-displacement equation.<br\>

$$\left[ \begin{array}{cc} K_{13} & K_{14} \\ K_{23} & K_{24} \\ K_{33} & K_{34} \\ K_{43} & K_{44} \\ K_{53} & K_{54}\\K_{63} & K_{64}\end{array} \right]\begin{Bmatrix} d_{3} \\ d_{4}\end{Bmatrix}=\begin{Bmatrix} F_{3} \\ F_{4}\end{Bmatrix}$$

By canceling the correct rows in the global stiffness matrix and the force matrix results forms the following equation<br\>

$$ \left[ \begin{array}{cc} k_{33} & k_{34}\\ k_{43}& k_{44} \end{array} \right]_{2x2} $$ $$\begin{Bmatrix} d_{3}\\d_{4} \end{Bmatrix}_{2x1}=\begin{Bmatrix}F_{3}\\F_{4} \end{Bmatrix}_{2x1}$$<br\><br\>

Next find the inverse of the stiffness matrix.<br\> $$det (k)=(3.0625)(2.6875)-(-2.175)(-2.175)=3.50$$<br\><br\>

$$ k^{-1}=\frac {1}{3.50}$$$$ \left[ \begin{array}{cc}2.6875&2.175\\ 2.175&3.0625\end{array} \right]=\left[ \begin{array}{cc}.7678&.6214\\ .6214&.8750\end{array} \right]$$<br\><br\> To find the displacement matrix, the inverse matrix is then multiplied by $$\begin{Bmatrix}{0} \\ {P} \end{Bmatrix}$$, where P=7.<br\><br\> $$\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=\left[ \begin{array}{cc}.7678&.6214\\ .6214&.8750\end{array} \right]\begin{Bmatrix}{0}\\{7}\end{Bmatrix}=\begin{Bmatrix}{4.352}\\{6.1271}\end{Bmatrix}$$<br\><br\>

$$k_{4x4}^{(e)}d_{4x1}^{(e)}=f_{4x1}^{(e)}$$ where $$e=1,2$$<br\><br\>

$$d^{(1)}=\begin{Bmatrix}{0}\\{0}\\{4.352}\\{6.1271}\end{Bmatrix}$$, $$f^{(1)}=\begin{Bmatrix}{f_{1}^{(1)}}\\{f_{2}^{(1)}}\\{f_{3}^{(1)}}\\{f_{4}^{(1)}}\end{Bmatrix}$$<br\><br\>

$$d^{(2)}=\begin{Bmatrix}{4.352}\\{6.1271}\\{0}\\{0}\end{Bmatrix}$$, $$f^{(2)}=\begin{Bmatrix}{f_{1}^{(2)}}\\{f_{2}^{(2)}}\\{f_{3}^{(2)}}\\{f_{4}^{(2)}}\end{Bmatrix}$$<br\><br\>

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

Given
A L2-ODE-CC system with the following conjugate roots:

$$ \lambda{_{1}}=-0.5+i2$$

$$ \lambda{_{2}}=-0.5-i2$$

The following are the initial conditions

$$ y(0)=1, y{}'(0)=0$$

Find
Find the homogeneous L2-ODE-CC in standard form in terms of $$y{_{h}}(t) $$ , the damping ratio  $$ \zeta $$  , and the solution with the above roots and initial conditions. Plot the solution.

Using the plot of the solution, find the (approximate) amplitudes of oscillation in a cycle (i.e. two successive maximum values of $$y_{h}(t)$$ , denoted by  $$y_{h}^{(n)}$$  and  $$y_{h}^{(n+1)}$$  with  $$ n $$  being the cycle number.)

Find the log decrements for $$ n=1,2,3 $$  using:

$$\delta^{(n)}=log\frac{y_{h}^{(n)}}{y_{h}^{(n+1)}}$$

Find the average log decrement:

$$ \delta\approx\frac13[\delta^((1))+\delta^((2))+\delta^((3))]$$

Verify the value of $$\delta$$  using:

$$\delta=\frac{2\pi \zeta}{\sqrt{1-\zeta^{2}}}$$

Find the quality factor $$ Q $$  and the loss factor  $$ \eta $$  for the system.

Homogeneous L2-ODE-CC in standard form
$$y{_{h}}''+ay{_{h}}'+by{_{h}}=0$$

$$\lambda^{2}+a\lambda+b=0$$

$$\Delta=a^{2}-4b< 0$$

So this system is underdamped.

'''$$\left ( \lambda-\lambda_{1}\right )\left ( \lambda-\lambda_{2} \right )=0 $$'''

Plug in the values for '''$$\lambda_{1},\lambda_{2} $$'''

'''$$\left ( \lambda-(-0.5+i2) \right )\left ( \lambda-(-0.5-i2) \right )=0 $$'''

Simplify everything inside the parentheses:

'''$$\left ( \lambda+0.5-i2\right )\left ( \lambda+0.5+i2 \right )=0 $$'''

Foil everything to get:

'''$$\lambda^{2}+0.5\lambda+i2\lambda+0.5\lambda+0.25+0.5i2-i2\lambda-0.5i2-(i2)^{2}=0 $$'''

Reduce everything:

$$\lambda^{2}+\lambda+0.25+4=0$$

$$\lambda^{2}+\lambda+4.25=0$$

The homogeneous L2-ODE-CC in standard form is:

$$ {y_{h}}''+{y_{h}}'+4.25y{_h}=0 $$

Damping Ratio
The correlation for damping ratio is given by:

$$ 2w\zeta=a$$

From the homogeneous L2-ODE-CC in standard form, it is known that:

$$ a=1$$

$$ b=2.25$$

Solving for $$ w $$  gives:

$$ w^{2}=b$$

$$ w=1.5$$

Plug in the value into:

$$ \zeta=\frac{1}{2w}$$

$$ \zeta=\frac{1}{3}$$

Solution and plot given the initial conditions
$$ y=c_{1}e^{(\alpha x)}cos(\beta x)+c_{2}e^{(\alpha x)}sin(\beta x)$$

From the conjugate roots, we know that:

$$ \alpha=-0.5, \beta= 2$$

Those values are plugged in:

$$ y=c_{1}e^{(-0.5 x)}cos(2 x)+c_{2}e^{(-0.5 x)}sin(2 x)$$

The first derivative is found to be:

$$y{}'=-0.5c_{1}e^{(-0.5 x)}cos(2 x)-2c_{1}e^{(-0.5x)}sin(2x)-0.5c_{2}e^{(-0.5 x)}sin(2 x)+2c_{2}e^{(-0.5x)}cos(2x)$$

The initial conditions are plugged in the two above equations to solve the equations.

The first equation gives:

$$ y(0)=c_{1}e^{(0)}cos(0)+c_{2}e^{(0)}sin(0)$$

This results in:

$$c_{1}=1$$

The second equation gives:

$$y{}'(0)=-0.5c_{1}e^{(0)}cos(0)-2c_{1}e^{(0)}sin(0)-0.5c_{2}e^{(0)}sin(0)+2c_{2}e^{(0)}cos$$

This reduces to:

$$ -0.5+2c_{2}=0$$

$$ c_{2}=0.25 $$

Plugging in the coefficients, the solution is:

$$y=e^{-0.5x}cos(2x)+0.25e^{-0.5x}sin(2x)$$



Approximate amplitudes of oscillation
Looking at the graph of the solution, the approximate amplitudes of oscillation are found by locating two successive maximum values of $$y_{h}(t)$$.

The amplitudes of oscillation are:

$$y_{h}^{(1)}=1$$

$$y_{h}^{(2)}=0.2$$

$$y_{h}^{(3)}=0.09$$

$$y_{h}^{(4)}=0.02$$

Log Decrements
The log decrements for $$ n=1,2,3 $$  are found using the equation:

$$\delta^{(n)}=log\frac{y_{h}^{(n)}}{y_{h}^{(n+1)}}$$

The log decraments are:

$$\delta^{(1)}=log\frac{1}{0.2}=0.699 $$

$$\delta^{(2)}=log\frac{0.2}{0.09}=0.347$$

$$\delta^{(3)}=log\frac{0.09}{0.02}=0.653$$

Average Log Decrement
The average log decrement is found using the equation:

$$ \delta\approx\frac13[\delta^{(1)}+\delta^{(2)}+\delta^{(3)}]$$

The values found previously are plugged in:

$$ \delta\approx\frac13[0.699+0.347+0.653]$$

$$ \delta\approx 1.699 $$

To check the answer, we use another equation:

$$\delta=\frac{2\pi (1/3)}{\sqrt{1-(1/3)^{2}}}$$

$$\delta= 2.221$$

Quality Factor
The quality factor is given by the equation:

$$Q=\frac{\sqrt{1-\zeta ^{2}}}{2\zeta}$$

$$Q=\frac{\sqrt{1-(1/3 ^{2}}}{2(1/3)}$$

$$Q=1.414$$

Loss Factor
The loss factor is the inverse of the quality factor:

$$\eta=\frac{1}{Q}$$

$$\eta=\frac{1}{1.414}$$

$$\eta=0.707$$

= Problem 3: Exercise 21, p. 103 Kim and Sankar: Matlab, Statics, and CALFEM = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Truss shown in figure.

Find
Determine the global stiffness matrix, the unknown displacement degrees of freedom, the reaction forces, and the member forces for the given truss system using Matlab. Verify the reaction and member forces with those computed using a Statics approach. Verify the global stiffness matrix, the unknown displacement degrees of freedom, the reaction forces, and the member forces using CALFEM.

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

Enter the x, y, and z coordinates (this program was designed solve 1D, 2D, and 3D problems, so the z coordinates must be entered. As long as all nodes have the same z coordinates, the solution will match a 2D problem). >> x=[  0  144  144  288    0 -144]; >> y=[ 108   0  108    0    0    0]; >> z=[  0    0    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    0             4     0             5 -1200             7   400             8     0            10     0            13     0            14     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. Only the known displacements are entered into this matrix. If the force for a particular degree of freedom is unknown, the its displacement is known. Because there is no force in the z direction, it is assumed that all displacements in the z direction are equal to 0. >> Q_in=[ 3 0 6 0            9 0            11 0            12 0            15 0            16 0            17 0            18 0];

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

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=[2 3 4 1 5 5 1 6 2];

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

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> Modulus=[10000 10000 10000 10000 10000 10000 10000 10000 10000];

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> Area=[3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159];

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

The stiffness matrix k: k = 1.0e+03 * Columns 1 through 8 441.5679        0 -111.7010   83.7757 -218.1660         0         0         0            0  416.5516   83.7757  -62.8318         0         0         0         0    -111.7010   83.7757  548.0329  -83.7757         0         0 -218.1660         0      83.7757  -62.8318  -83.7757  353.7198         0 -290.8880         0         0    -218.1660         0         0         0  329.8669  -83.7757 -111.7010   83.7757            0         0         0 -290.8880  -83.7757  353.7198   83.7757  -62.8318            0         0 -218.1660         0 -111.7010   83.7757  329.8669  -83.7757            0         0         0         0   83.7757  -62.8318  -83.7757   62.8318            0         0 -218.1660         0         0         0         0         0            0 -290.8880         0         0         0         0         0         0    -111.7010  -83.7757         0         0         0         0         0         0     -83.7757  -62.8318         0         0         0         0         0         0     Columns 9 through 12 0        0 -111.7010  -83.7757            0 -290.8880  -83.7757  -62.8318    -218.1660         0         0         0            0         0         0         0            0         0         0         0            0         0         0         0            0         0         0         0            0         0         0         0     436.3319         0 -218.1660         0            0  290.8880         0         0    -218.1660         0  329.8669   83.7757            0         0   83.7757   62.8318

The displacement matrix Q: Q = 1.0000   9.9122       2.0000  -17.9909       3.0000         0       4.0000    7.3339       5.0000  -26.2033       6.0000         0       7.0000    6.2452       8.0000  -23.1093       9.0000         0      10.0000   12.8343      11.0000         0      12.0000         0      13.0000    3.6669      14.0000  -17.9909      15.0000         0      16.0000         0      17.0000         0      18.0000         0

Every 3 rows correspond to a separate global node. (ie. global node 1 experiences a positive displacement in the x direction, a negative displacement in the y direction, and no displacement in the z direction) It can be seen that global node 4 (rows 10-12) experiences only a positive displacement in the x direction.

The Force matrix, F F = 1.0e+03 * -0.0000      0.0000            0       0.0000      -1.2000            0       0.4000       0.0000            0      -0.0000       0.9000            0       0.0000            0            0      -0.4000       0.3000            0

This shows that at global node 4, a positive vertical force of 900 lb is exerted on the truss. Additionally, at global node 6 a negative horizontal force of 400 lb and a positive vertical force of 300 lb is exerted on the truss. This is in accord with a force balance of the entire truss.

The local displacement matric Q_bar

Q_bar = -21.5891  18.8618  -12.8343   17.9909    3.6669   -3.6669    9.9122         0  -26.2033     -18.7243   10.2674   -7.3339   17.9909    7.3339         0    6.2452   -2.8648  -23.1093

The internal spring force matrix F_bar F_bar = 1.0e+03 * -0.5000   1.5000   -1.2000         0   -0.8000   -0.8000    0.8000    0.5000   -0.9000       0.5000   -1.5000    1.2000         0    0.8000    0.8000   -0.8000   -0.5000    0.9000

Matlab Code

Using Statics (Hand Calculation)
Begin by solving for the reaction force at global node 4 using a moment balance about global node 6.

$$\Sigma M_{6}=-24R_{2y}-9R_{3x}+36R_{4y}$$

$$R_{4y}=(24*1200+9*400)/36=900$$

Sequentially move from node to node solving the force balance equations at each node.

The following table shows the equivalence between the reaction forces computed by Matlab and those computed from the Statics approach.

The following table shows the equivalence between the member forces computed by Matlab and those computed from the Statics approach.

CALFEM Verification
Using CALFEM the following code was created in Matlab to calculate the stiffness matrix, compute the unknown displacement degrees of freedom, the reaction forces, and the member forces.

Comparing Results

These results from both the Matlab program and the CALFEM program show that the Matlab code is working correctly. CALFEM uses the known information to compute the unknowns therefore the known information for displacement, reaction forces, and member forces are not included in their respective matrices. Also the Matlab program outputs the Z components of these matrices.

Stiffness Matrix

CALFEM

Matlab

Displacements

CALFEM

Matlab Results from both CALFEM and Matlab. d1 = 9.91 inches d2  = -17.99 inches

d3 = 7.33 inches d4  = -26.20 inches

d5 = 6.24 inches d6  = -23.10 inches

d7 = 12.8 inches d8  =      0 inches

d9 = 3.66 inches d10 = -17.99 inches

d11 =   0 inches d12 =      0 inches

Reaction Forces

CALFEM

Matlab Results from both CALFEM and Matlab. F1 =    0 lbf F2  =     0 lbf

F3 =    0 lbf F4  = -1200 lbf

F5 =  400 lbf F6  =     0 lbf

F7 =    0 lbf F8  =   900 lbf

F9 =    0 lbf F10 =     0 lbf

F11 = -400 lbf F12 =  300 lbf

Member Forces

CALFEM

Matlab Results from both CALFEM and Matlab. F1 =  500 lbf

F2 = -1500 lbf

F3 = 1200 lbf

F4 =    0 lbf

F5 =  800 lbf

F6 =  800 lbf

F7 = -800 lbf

F8 = -500 lbf

F9 =  900 lbf

= Problem 4: (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 Contidions
 * $$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)$$
First find the homogenous solution $$y_h(t)$$. $$y_h(t) = C_1 e^{(\alpha t)} sin(\beta t)+C_2 e^{(\alpha t)} cos(\beta t)$$, $$\lambda_1 = \bar \lambda_2 = \alpha + i\beta$$, where $$\alpha = -0.5, \beta = 2$$ $$y_h(t)=C_1 e^{(-0.5 t)} sin(2 t)+C_2 e^{(-0.5 t)} cos(2 t)$$

By setting $$y''(t) + y'(t) + 4.25y(t) = f_0cos(\bar \omega t)$$ and solving for the particular solution:

$$y_p(t) = \frac{(f_0 cos(t \bar \omega)(-\bar \omega^2 sin^2(2 t)-\bar \omega^2 cos^2(2 t)+4.25)+f_0 \bar \omega sin^2(2 t) sin(t \bar \omega)+f_0 \bar \omega cos^2(2 t) sin(t \bar \omega))}{(\bar \omega^4-7.5 \bar \omega^2+18.0625)}$$

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$$

$$y_p(t) = \mathbb A \frac{( cos(t 0.9\omega)(- 0.81\omega^2 sin^2(2 t)- 0.81\omega^2 cos^2(2 t)+4.25)+ 0.9\omega sin^2(2 t) sin(t 0.9\omega)+ 0.9\omega cos^2(2 t) sin(t 0.9\omega))}{( 0.6561\omega^4-6.075\omega^2+18.0625)}$$ $$\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$$

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) = \mathbb A \frac{( cos(t 0.9\omega)(- 0.81\omega^2 sin^2(2 t)- 0.81\omega^2 cos^2(2 t)+4.25)+ 0.9\omega sin^2(2 t) sin(t 0.9\omega)+ 0.9\omega cos^2(2 t) sin(t 0.9\omega))}{( 0.6561\omega^4-6.075\omega^2+18.0625)}$$ $$y(t)=C_1 e^{(-0.5 t)} sin(2 t)+C_2 e^{(-0.5 t)} cos(2 t)$$
 * $$+\mathbb A \frac{( cos(t 0.9\omega)(- 0.81\omega^2 sin^2(2 t)- 0.81\omega^2 cos^2(2 t)+4.25)+ 0.9\omega sin^2(2 t) sin(t 0.9\omega)+ 0.9\omega cos^2(2 t) sin(t 0.9\omega))}{( 0.6561\omega^4-6.075\omega^2+18.0625)}$$

Simplify with: $$\omega = \sqrt{4.25}$$ $$1 = C_2 + \frac{2.10032(-0.81 \omega^2+4.25)}{4.09456}$$

$$C_2 = 0.58579$$ $$0 = C_1(2)+0.58579(-0.5)+ \mathbb A \frac{3.4425}{4.09456}$$ $$C_1 = -0.73648$$ The complete solution is: $$y(t)=-0.73648 e^{(-0.5 t)} sin(2 t)+0.58579 e^{(-0.5 t)} cos(2 t)$$
 * $$+\mathbb A \frac{( cos(t 0.9\omega)(- 0.81\omega^2 sin^2(2 t)- 0.81\omega^2 cos^2(2 t)+4.25)+ 0.9\omega sin^2(2 t) sin(t 0.9\omega)+ 0.9\omega cos^2(2 t) sin(t 0.9\omega))}{( 0.6561\omega^4-6.075\omega^2+18.0625)}$$

Plots








= Problem 5: Finite Element Analysis Of A Spring System = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Spring arrange shown in the figure.

Find
Find the solution using matlab and hand calculations

Solution
From Hooke's law for a spring system we know:



$$\displaystyle f=kd$$

$$\displaystyle f = force$$

$$\displaystyle k = stiffness$$

$$\displaystyle d = displacement$$

In matrix form:

$$ \begin{bmatrix} f_{1x} \\ f_{2x} \\ \end{bmatrix} = \begin{bmatrix} k & -k \\ -k & k \\ \end{bmatrix} \begin{bmatrix} d_{1x} \\ d_{2x} \\ \end{bmatrix} $$

First we create a matrix to show the system's Degree of Freedoms.

The first column represents each spring(#1,#2,#3). The second and third column represent the connectivity of each node.

$$ DOF = \begin{bmatrix} 1 & 1 & 2\\ 2 & 2 & 3\\ 3 & 2 & 3\\ \end{bmatrix} $$

The following is a matrix of the Force vector

$$ F = \begin{bmatrix} F_{1x}\\ F_{2x}\\ F_{3x}\\ \end{bmatrix} $$

Next we need the element stiffness matrix for each spring:

Element stiffness matrix for Spring #1:

$$ k_{1} = \begin{bmatrix} 3000 & -3000 \\ -3000 & 3000 \\ \end{bmatrix}= \begin{bmatrix} 3000 & -3000 & 0\\ -3000 & 3000 & 0\\ 0 & 0 & 0\\ \end{bmatrix} $$

Element stiffness matrix for Spring #2:

$$ k_{2} = \begin{bmatrix} 1500 & -1500 \\ -1500 & 1500 \\ \end{bmatrix}= \begin{bmatrix} 0 & 0 & 0\\ 0 & 1500 & -1500\\ 0 & -1500 & 1500\\ \end{bmatrix} $$

Element stiffness matrix for Spring #3:

$$ k_{3} = \begin{bmatrix} 3000 & -3000 \\ -3000 & 3000 \\ \end{bmatrix}= \begin{bmatrix} 0 & 0 & 0\\ 0 & 3000 & -3000\\ 0 & -3000 & 3000\\ \end{bmatrix} $$

To obtain the global stiffness matrix, we add up each element stiffness matrix:

$$ K = \begin{bmatrix} 3000 & -3000 & 0\\ -3000 & 3000 & 0\\ 0 & 0 & 0\\ \end{bmatrix}+ \begin{bmatrix} 0 & 0 & 0\\ 0 & 1500 & -1500\\ 0 & -1500 & 1500\\ \end{bmatrix}+ \begin{bmatrix} 0 & 0 & 0\\ 0 & 3000 & -3000\\ 0 & -3000 & 3000\\ \end{bmatrix}= \begin{bmatrix} 3000 & -3000 & 0\\ -3000 & 7500 & -4500\\ 0 & -4500 & 4500\\ \end{bmatrix} $$

Plugging each matrix into Hooke's Law we obtain:

$$ \begin{bmatrix} 3000 & -3000 & 0\\ -3000 & 7500 & -4500\\ 0 & -4500 & 4500\\ \end{bmatrix} \begin{bmatrix} d_{1x} \\ d_{2x} \\ d_{3x} \\ \end{bmatrix}= \begin{bmatrix} F_{1x}\\ F_{2x}\\ F_{3x}\\ \end{bmatrix} $$

Because node one is attached to the wall, it's displacement in zero:


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

$$\displaystyle d_{1x}=0$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Now we solve the follwing equation:

$$ \begin{bmatrix} 3000 & -3000 & 0\\ -3000 & 7500 & -4500\\ 0 & -4500 & 4500\\ \end{bmatrix} \begin{bmatrix} 0 \\ d_{2x} \\ d_{3x} \\ \end{bmatrix}= \begin{bmatrix} F_{1x}\\ F_{2x}\\ F_{3x}\\ \end{bmatrix} $$

Solving the System we obtain the following equations:


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

$$\displaystyle 3000*0 -3000*d_{2x}+0*d_{3x} = F_{1x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


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

$$\displaystyle -3000*0+7500*d_{2x}-4500*d_{3x} = F_{2x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


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

$$\displaystyle 0*0-4500*d_{2x}+4500*d_{3x} = F_{3x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Simplifying the above equations we obtain:
 * {| style="width:100%" border="0"

$$\displaystyle -3000*d_{2x} = F_{1x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


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

$$\displaystyle 7500*d_{2x}-4500*d_{3x} = F_{2x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }


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

$$\displaystyle -4500*d_{2x}+4500*d_{3x} = F_{3x}$$
 * style="width:97%" |
 * style="width:97%" |
 * }

Solving we obtain:


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

$$\displaystyle F_{1x}= 40 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Plugging this value in


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

$$\displaystyle {\frac{40}{-3000}}=d_{2x}=-0.01333 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Using substitution for the remaining values, we obtain magnitudes of:

$$ d_{1x}= 0 $$

$$ d_{2x}= 0.01333 $$

$$ d_{3x}= 0.01333 $$

$$ F_{1x}= 40 $$

$$ F_{2x}= 20 $$

$$ F_{3x}= 40 $$

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

Given
Spring-rigid-body arrangement as shown. Global nodes 1 and 5 are supports. A force F=-1000 N is applied to global node 3 in the x direction. Motion is constrained to the horizontal direction. Spring constants for springs 1-6 are given as 500, 400, 600, 200, 400, and 300 N/mm respectively.

Find
Find the displacements of global nodes 2, 3, and 4, the forces in the springs, and the reactions at the walls using Matlab/by hand and verify using CALFEM.

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

Enter the x, y, and z coordinates (this program was designed solve 1D, 2D, and 3D problems, so the y and z coordinates must be entered. As long as all nodes have the same y and z coordinates, the solution will match a 1D problem). The exact x coordinates for this problem are not important. As long as the relative location of each nodes remains the same, the solution will be correct. >> x=[0 1 2 3 4]; >> y=[0 0 0 0 0]; >> z=[0 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=[4 0; 7 -1000; 10 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. Only the known displacements are entered into this matrix. If the force for a particular degree of freedom is unknown, the its displacement is known. >> 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];

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

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 2 2 1 3 4];

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

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> spring_const=[500 400 600 200 400 300];

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

The stiffness matrix k: k = 700       -500        -200           0           0        -500        1500        -600        -400           0        -200        -600        1200        -400           0           0        -400        -400        1100        -300           0           0           0        -300         300

The displacement matrix Q: Q = 1.0000        0    2.0000         0    3.0000         0    4.0000   -0.8542    5.0000         0    6.0000         0    7.0000   -1.5521    8.0000         0    9.0000         0   10.0000   -0.8750   11.0000         0   12.0000         0   13.0000         0   14.0000         0   15.0000         0

This shows that the first degree of freedom (x direction) for global node 2 experiences a negative displacement of .8542 mm, the first degree of freedom (x direction) for global node 3 experiences a negative displacement of 1.5521 mm, and the first degree of freedom (x direction) of global node 4 experiences a negative displacement of .875 mm.

The Force matrix, F F = 1.0e+03 * 0.7375        0         0   -0.0000         0         0   -1.0000         0         0         0         0         0    0.2625         0         0

This shows that the first degree of freedom (x direction) for global node 1 exerts a positive force of 737.5 N and the first degree of freedom (x direction) for global node 5 exerts a positive force of 262.5 N.

The internal spring force matrix F_bar F_bar = 427.0833   8.3333  418.7500  310.4167 -270.8333 -262.5000 -427.0833   -8.3333 -418.7500 -310.4167  270.8333  262.5000

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 2 2 4           3 2 3            4 1 3            5 3 4            6 4 5];

Construct an empty stiffness matrix K that has side lengths equal to the number of dimensions times the number of global nodes. For this problem, motion is constrained to 1 dimension and there are 5 global nodes. >> K=zeros(5,5);

Construct the force matrix F. The forces corresponding to global nodes 2, 3, and 4 are known; the forces corresponding to global nodes 1 and 5 are unknown. >> F=zeros(5,1); F(3)=-1000;

Construct the displacement matrix Q. Global nodes 1 and 5 have known displacements of 0. >> Q=[1 0 5 0];

The function spring1e generates the element stiffness matrices. Its parameter is the spring stiffness for each respective element. The stiffness are given in the problem statement.

The commands: >> Ke2=spring1e(200); >> Ke3=spring1e(300); >> Ke4=spring1e(400); >> Ke5=spring1e(500); >> Ke6=spring1e(600);

Yield the following stiffness matrices: >> Ke2=[200 -200 -200 200];   >> Ke3=[300 -300 -300 300];   >> Ke4=[400 -400 -400 400];   >> Ke5=[500 -500 -500 500];   >> Ke6=[600 -600 -600 600];

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

This results in the following stiffness matrix. >> K = 700       -500        -200           0           0        -500        1500        -600        -400           0        -200        -600        1200        -400           0           0        -400        -400        1100        -300           0           0           0        -300         300

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     -0.8542      -1.5521      -0.8750            0   r = 737.5000     -0.0000       0.0000       0.0000     262.5000

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

The function extract is used to computer the displacement of each spring, which is needed to determine the force in each spring. >> ed1=extract(Edof(1,:),q); >> ed2=extract(Edof(2,:),q); >> ed3=extract(Edof(3,:),q); >> ed4=extract(Edof(4,:),q); >> ed5=extract(Edof(5,:),q); >> ed6=extract(Edof(6,:),q);

This results in the following displacements. The number following 'ed' corresponds to the element number. The first value represents the displacement of the left side of each spring. The second column corresponds to the displacement of the right side of each spring. ed1 =      0   -0.8542 ed2 = -0.8542  -0.8750 ed3 = -0.8542  -1.5521 ed4 =      0   -1.5521 ed5 = -1.5521  -0.8750 ed6 = -0.8750        0

The function spring1s can be used to determine the force in each spring. Its input parameters are the spring stiffness and the spring displacements. >> es1=spring1s(500,ed1) >> es2=spring1s(400,ed2) >> es3=spring1s(600,ed3) >> es4=spring1s(200,ed4) >> es5=spring1s(400,ed5) >> es6=spring1s(300,ed6)

This results in the following spring forces. es1 = -427.0833 es2 =  -8.3333 es3 = -418.7500 es4 = -310.4167 es5 = 270.8333 es6 = 262.5000

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

=References=