User:EML4500.f08.A-team.kirley

Loops in MatLab
A loop is a programming tool in Matlab that is used to do repetitive work. It saves the user time. In the tutorial, matrix a represents the probabilities that two island populations will move or stay put on their side of the island. x is a column vector representing the initial populations.

a=    0.8  0.1 0.2 0.9

x=    1 Subscript text             0

Successive population states can be predicted by multiplying the two matrices. Using a loop, the operation can be repeated a defined number of times using the following form:

>> a = [ 0.8 0.1; 0.2 0.9 ] >> x = [ 1; 0 ] >> for i = 1:20, x = a*x, end

This operation will be repeated 20 times, and the populations of the island can be predicted 20 time units away.

Meeting 8, September 12, 2008
In compact notation:

[Kij]6x6[dj]6x1=[Fi]6x1 (more generally n x n)

ΣKijdj = Fi     i = 1, 2, ... , 6

Kn x n = [Kij} = global stiffness matrix

dn x 1 = {dj}n x 1 = global displacement matrix

Fn x 1 = {Fi}n x 1 = global force matrix

Recall elem FD relation:

k(l)4x4 d(l)4x1 = f(l)4x1

Where k, d, and f are the element stiffness, displacement, and force matrices, respectively.

Assembly of global matrices from element matrices:

Identify the correspondence between element displacement degrees of freedom and global displacement degrees of freedom

Global level: {d1, -d2, ... , d6}

Elem level:

*Elem 1: {d(1)1, d(1)2, d(1)3, d(1)3, d(1)4

*Elem 2: {d(2)1, d(2)2, d(2)3, d(2)3, d(2)4

d1 = d(1)1

d2 = d(1)2

d3 = d(1)3 = d(2)1

d4 = d(1)4 = d(2)2

d5 = d(2)3

d6 = d(2)4

Visualization of assembly:



Meeting 9, September 15, 2008


Kundefinedk(1)11 = 9/16

Kundefinedk(1)12 = .3247595

K33 = k(1)33 + k(2)11 = 9/16 + 5/2 = 3.0625

K34 = K34 = k(1)34 + k(2)12 = .3247595 + 5/2

The rest of the stiffness matrix can be computed similarly

Elimination of known dofs:

reduce global FD rel.

d1 = d2 = d5 = d6 = 0

EML4500.f08.A-team.kirley 20:25, 26 September 2008 (UTC)

Infinitesimal displacement


2 unknowns (XD YD)





PQ = (PQ)i^{~} = (X-X_{P})i + (Y-Y_{P})j $$



(PQ)i^{~} = PQ[ \cos{\theta}i + \sin{\theta}j ] $$

Therefore



X - X_{P} = PQ\cdot \cos{\theta} $$



Y - Y_{P} = PQ\cdot \sin{\theta} $$



\frac {Y-Y_{P}}{X-X_{P}} = \tan{\theta} $$



Y-Y{P} = (\tan{\theta})\cdot(X-X_{P}) $$

Equation for a line perpendicular to above line, passing P:



Y-Y_{P} = \tan{(\theta + \frac{\pi}{2})}(X-X{P}) $$

 Three-bar Truss System 



Local node numbering by convenience in assembly of K:



EML4500.f08.A-team.kirley 20:54, 8 October 2008 (UTC)

Connectivity Array "conn"
Consider the two-bar truss system (5-6):



conn(e,j) = global node number of local node j of element e

Location Matrix Master Array "lmm"


lmm(i,j) = eq. number (global dof #) for the element stiffness coefficient corresponding to the jth local dof #

Meeting 20 10/10/08
Note: consider case: $$\tilde{d}^{(e)}_4 \neq 0$$

$$ \tilde{d}^{(e)}_1 {=} \tilde{d}^{(e)}_2 {=} \tilde{d}^{(e)}_3 {=} 0$$

$$ \underline{\tilde{f}}^{(e)}_{4x1} {=} \underline{\tilde{k}}^{(e)}_{4x4} \cdot \underline{\tilde{d}}^{(e)}_{4x1} {=} \underline{0}_{4x1} \leftarrow 4^{th}$$ column of $$\underline{\tilde{k}}^{(e)}$$

Interpretation of trans. dofs

p. 19-3: $$ \underline{\tilde{d}}^{(e)} {=} \underline{\tilde{T}}^{(e)} \cdot \underline{d}^{(e)}$$

Similarly  $$ \underline{\tilde{f}}^{(e)} {=} \underline{\tilde{T}}^{(e)} \cdot \underline{f}^{(e)}$$

Also:  $$ \underline{\tilde{k}}^{(e)} \cdot \underline{\tilde{d}}^{(e)} {=} \underline{\tilde{f}}^{(e)} $$

$$ \rightarrow \underline{\tilde{k}}^{(e)} \cdot \underline{\tilde{T}}^{(e)} \cdot \underline{\tilde{d}}^{(e)} {=} \underline{\tilde{T}}^{(e)} \cdot \underline{\tilde{f}}^{(e)}$$

If $$\underline{\tilde{T}}$$ is invertible, then:

$$ [\underline{\tilde{T}}^{(e)-1} \cdot \underline{\tilde{k}}^{(e)} \cdot \underline{\tilde{T}}^{(e)}] \cdot \underline{d}^{(e)} {=} \underline{f}^{(e)} $$

$$ \underline{\tilde{T}}^{(e)}$$ block diag. matrix (p. 19-3)

Consider a general block-diag. matrix:



Question: What is A-1

Simpler example: Diag. matrix:



$$ {=} diag[d_{11}, d_{22} , d_{33} , ... , d_{nn}] $$

$$ \underline{B}^{-1} = Diag[\frac{1}{d_{11}}, \frac{1}{d_{22}} , ... , \frac{1}{d_{nn}}]$$

Assuming $$ d_{ii} \neq 0 $$ for $$ i {=} 1, ... , n$$

For a block-diag matrix A:

$$ \underline{A} {=} Diag[\underline{D}_1, ... , \underline{D}_s]$$

$$ \underline{A}^{-1} {=} Diag[\underline{D}_1^{-1}, ... , \underline{D}^{-1}_{s}]$$

$$ \underline{\tilde{T}}^{(e)-1} {=} Diag[\underline{R}^{(e) -1}, ... , \underline{R}^{(e) -1}] $$

p. 19-2: $$ \underline{R}^{(e) T} {=} \begin{vmatrix} l^{(e)} & -m^{(e)} \\ m^{(e)} & l^{(e)} \end{vmatrix} $$

$$ \underline{R}^{(e)T}_{2x2} \cdot \underline{R}^{(e)}_{2x2} {=} \begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix}_{2x2} {=} \underline{I}_{2x2} $$

$$ \rightarrow \underline{R}^{(e)-1} = \underline{R}^{(e)T} $$

$$ \rightarrow \underline{\tilde{T}}^{(e)-1} = Diag[\underline{R}^{(e)T}, \underline{R}^{(e)T}] $$

$$ {=} \underbrace{ (Diag[ \underline{R}^{(e)}, \underline{R}^{(e)}]) }_{\underline{\tilde{T}}^{(e)}} $$

$$ \underline{\tilde{T}}^{(e)-1} {=} \underline{\tilde{T}}^{(e)T} $$

p. 20-1: $$ \underbrace{ [\underline{\tilde{T}}^{(e)T} \cdot \underline{\tilde{k}}^{(e)} \cdot \underline{\tilde{T}}^{(e)}]}_{\underline{k}^{(e)}} \cdot \underline{d}^{(e)} {=} \underline{f}^{(e)} $$

HW: Verify: $$ \underline{k}^{(e)} {=} \underline{\tilde{T}}^{(e)T} \cdot \underline{\tilde{k}}^{(e)} \cdot \underline{\tilde{T}}^{(e)}$$

EML4500.f08.A-team.kirley 01:19, 30 October 2008 (UTC)

Meeting 27 10/29/08
$$ \hat{\underline{w}}_{2x1} {=}$$ virtual axial displacement, corresponding to $$ q^{(e)}_{2x1} $$

$$ \underline{w}_{4x1} {=} $$ virtual displacement in global coordinate system, corresponding to $$ d^{(e)}_{4x1} $$

Replace Equations (3) and (4) into (2):

$$ \hat{\underline{w}}_{2x1} \cdot (\hat{\underline{k}}^{(e)} \underline{q}^{(e)} - \underline{p}^{(e)}) {=} 0_{1x1} $$ (2)

$$ \Rightarrow (\underline{T}^{(e)}\underline{w}) \cdot [\hat{\underline{k}}^{(e)}(\underline{T}^{(e)}\underline{d}^{(e)}) - \underline{p}^{(e)}] {=} 0 $$ for all $$ w_{4x1} $$ (5)

Recall: $$ (\underline{A}\underline{B})^T {=} \underline{B}^T\underline{A}^T $$ (6)

Recall: $$ \underline{a}_{nx1} \cdot \underline{b}_{nx1} {=} \underbrace{ \underline{a}^T_{1xn}\underline{b}_{nx1} }_{1x1 (scalar)} $$ (7)

Apply (7) and (6) in (5):

$$ \underbrace{(\underline{T}^{(e)} \underline{w})^T}_{(7)}[\underline{\hat{k}}^{(e)}(\underline{T}^{(e)}\underline{d}^{(e)}) - \underline{p}^{(e)}] {=} 0_{1x1} $$ for all $$ \underline{w}_{4x1} $$

$$ \Rightarrow \underbrace{\underline{w}^T\underline{T}^{(e)T}}_{(6)}[-] {=} 0 $$

$$ \Rightarrow \underbrace{\underline{w}}_{(7)} \cdot [(\underbrace{\underline{T}^{(e)T} \underline{\hat{k}}^{(e)} \underline{T}^{(e)} \underline{d}^{(e)}}_{\underline{k^{(e)}}}) - \underbrace{\underline{T}^{(e)T} \underline{p}^{(e)}}_{\underline{f}^{(e)}}] {=} 0 $$ for all $$ \underline{w}_{4x1} $$

$$ \Rightarrow \underline{w} \cdot [\underline{k}^{(e)}\underline{d}^{(e)} - \underline{f}^{(e)}] {=} 0 $$ for all $$ \underline{w} $$

$$ \Rightarrow \underline{k}^{(e)}\underline{d}^{(e)} {=} \underline{f}^{(e)} $$

So far, only discrete cases have been discussed. Now continuous cases will be discussed.

Motivational model problem: Elastic bar with varying A(x), E(x), subjected to varying axial load (distributed) and concentrated load and inertia force (dynamics). Concentrated load and inertia force are time dependent.



Meeting 33: Wednesday, 12 November 2008. EML4500
FEM vs. PVW (cont.)

$$ u(x_{i+1}) {=} d_{i+1} $$

p. 31-1: interpolation for u(x)

Apply same interpolation for w(x), i.e.

$$ w(x) {=} N_i(x)w_i + N_{i+1}(x)w_{i+1} $$

Element stiffness matrix for element i:

$$ \int_{x_i}^{x_{i+1}} [N_i^'w_i + N_{i+1}^'w_{i+1}](EA)[N_i^'d_i + N_{i+1}^'d_{i+1}]\, dx $$

$$ N_i^' :{=} \frac{dN_i(x)}{dx} $$

Likewise for $$ N_{i+1}^' $$:

$$ N_{i+1}^' :{=} \frac{dN_{i+1}(x)}{dx} $$

Note: $$ u(x) {=} \underbrace{ \begin{vmatrix} N_i(x) & N_{i+1}(x) \end{vmatrix}}_{N(x)_{1x2}} \begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix}_{2x1} $$

$$ \frac{du(x)}{dx} = \underbrace{ \begin{vmatrix} N_i^'(x) & N_{i+1}^'(x) \end{vmatrix}}_{\underline{B}(x)_{1x2}} \begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix}_{2x1} $$

Similarly: $$ w(x) {=} \underline{N}(x) \begin{Bmatrix} W_i \\ W_{i+1} \end{Bmatrix} $$

$$ \frac{dW(x)}{dx} {=} \underline{B}(x) \begin{Bmatrix} W_i \\ W_{i+1} \end{Bmatrix} $$

Recall the element degree of freedom conventions:



$$ \begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix} {=} \begin{Bmatrix} d_1^{(i)} \\ d_2^{(i)} \end{Bmatrix} {=} \underline{d}^{(i)} $$

$$ \begin{Bmatrix} w_i \\ w_{i+1} \end{Bmatrix} {=} \begin{Bmatrix} w_1^{(i)} \\ w_2^{(i)} \end{Bmatrix} {=} \underline{w}^{(i)} $$

$$ \beta {=} \int_{x^i}^{x^{i+1}} \underbrace{(\underline{B} \underline{w}^{(i)})}_{1x1} \underbrace{ (EA)}_{1x1} \underbrace{(\underline{B}\underline{d}^{(i)})}_{1x1}, dx (scalar) $$

$$ {=} \underline{w}^i \cdot (\underline{k}^{(i)} \underline{d}^{(i)}) $$

$$ \beta {=} \int_{x^i}^{x^{i+1}} (EA) \underbrace{(\underline{B} \underline{w}^{(i)})}_{1x1} \cdot \underbrace{(\underline{B}\underline{d}^{(i)})}_{1x1} dx $$

$$ (\underline{B}\underline{w}^{(i)})^T(\underline{B}\underline{d}^{(i)}) $$

$$ (\underline{B}\underline{w}^{(i)})^T {=} \underline{w}^{(i)T} \underline{B}^T {=} \underline{w}^{(i)} \cdot \underline{B}^T $$

\beta {=} \underline{w}^{(i)} \cdot (\int \underline{B}^T(EA)\underline{B} dx )d^{(i)}

$$ \underline{k}^{(i)}_{2x2} {=} \int_{x_i}^{x_{i+1}} \underbrace{\underline{B}^T(x)}_{2x1} \underbrace{(EA)}_{\underbrace{(x)}_{1x1}} \underbrace{\underline{B}(x)}_{1x2} dx $$

$$ B(x) {=} \begin{vmatrix} HW & \frac{1}{L^{(i)}} \end{vmatrix} $$

$$ L^{(i)} {=} x_{i+1} - x_i $$ (Length of element i)

HW 6: Consider EA = const.

$$ k^{(i)} {=} \frac{EA}{L^{(i)}} \begin{vmatrix} 1 & -1 \\ -1 & 1 \end{vmatrix} $$

Transfer of variable coordinates from x to $$ \tilde{x} $$

$$ \tilde{x} :{=} x - x_i $$

$$ d\tilde{x} {=} dx $$

$$ k^{(i)} {=} \int^{\tilde{x} {=} L^{(i)}}_{\tilde{x} {=} 0} \underline{B}^T(\tilde{x})(EA)(\tilde{x})(\underline{B}(\tilde{x}) d\tilde{x} $$

HW6: Find expression for $$ \underline{k}^{(i)} $$ using above:



$$ A(\tilde{x}) {=} N_1^{(i)}(\tilde{x})A_1 + N_2^{(i)}(\tilde{x})A_2 $$

$$ E(\tilde{x}) {=} N_1^{(i)}(\tilde{x})E_1 + N_2^{(i)}(\tilde{x})E_2 $$

$$ \underline{k}^{(i)} {=} ? $$

Meeting 40: Friday, 5 December 2008. EML4500
Note: Eq. (2) (39-2), dimension analysis:

$$ [u] {=} L $$

(31-3): $$ [N_1] {=} [N_4] {=} 1 $$

$$ Eq. (2) (39-2) \underbrace{[N_1]}_{1} \underbrace{[\tilde{d}_1]}_{L} + \underbrace{[N_4]}_{1} \underbrace{[\tilde{d}_4]}_{L} $$

$$ [V] {=} L $$

$$ \underbrace{[N_2]}_{1} \underbrace{[\tilde{d}_2]}_{L} {=} L $$ (displ. transv.)

$$ \underbrace{[N_2]}_{1} \underbrace{[\tilde{d}_2]}_{L} {=} L $$ (rot.)

$$ \underbrace{[N_2]}_{1} \underbrace{[\tilde{d}_2]}_{L} and \underbrace{[N_2]}_{1} \underbrace{[\tilde{d}_2]}_{L} {=} L $$ are to be done for homework

Derivation of the beam shape functions

$$ N_2, N_3, N_5, N_6 $$ (38-4)

p. (38-3) plots

Recall: Governing PDE for beams

p. (37-4), Eq. (1), without $$ f_t $$ (distributed transverse load) and without inertia force $$ m\ddot{v} $$ (static case):

$$ \frac{\delta^2}{\delta x^2} \left \{(EI) \frac{\delta^2v}{\delta x^2}\right \} {=} 0 $$

Further, consider constant EI:

$$ \frac{\delta^4}{\delta x^4} v {=} 0 $$ Integrate 4 times to get 4 constants.

$$ \Rightarrow V(x) {=} C_0 + C_1x^1 + C_2x^2 + C_3x^3 $$

To obtain $$ N_2(x) $$ ($$ \tilde{x} $$ {=} x for simplicity)

V(0) = 1,  V(L) = 0

V'(0)   =   V'(L) = 0

Use above boundary conditions to solve for $$ L_0, ... , L_3 $$

V(0) = 1 = C0

(1) $$ V(L) {=} 1 + C_1(L) + C_2(L)^2 + C_3(L)^3 {=} 0 $$

$$ V^{\prime}(x) {=} C_1 + 2C_2x + 3C_3x^2 $$

$$ V^{\prime}(0) {=} C_1 {=} 0 $$

(2) $$ V^{\prime}(L) {=} 2C_2(L) + 3C_3(L)^2 {=} 0 $$

(2) $$ \rightarrow C_3 {=} - \frac{2}{3} \frac{C_2}{L} $$

(1) $$ \rightarrow 0 {=} 1 + C_2L^2 + (\frac{-2}{3} \frac{C_2}{L})L^3 {=} 1 + C_2L^2[1 - \frac{2}{3}] $$

$$ C_2 {=} - \frac{3}{L^2} $$

$$ C_3 {=} - \frac{2}{3} \frac{1}{L}(- \frac{3}{L^2}) {=} \frac{2}{L^3} $$

Compare with expression for N2 on p. (38-4)

For N3: boundary conditions         $$ \tilde{d}_3 $$ rotation

V(0) = V(L) = 0

V'(0) = 1, V'(L) = 0

For N5: V(0) = 0, V(L) = 1    $$ \tilde{d}_5 $$ displacement

V'(0) = 0, V'(L) = 0

For N6: $$ d_6 $$ rotation

V(0) = V(L) = 0

V'(0) = 0, V'(L) = 1

See p. (39-1) for plots of N5, N6

Derive coefficient in $$ \underline{\tilde{k}} $$ (elem. stiffness matrix)

Coefficient for EA: Done

Coefficient for EI: To Be Done

\tilde{k}_{22} {=} \frac{12EI}{L^3} \int^L_0 \frac{d^2N_2}{dx^2}(EI)\frac{d^2N_2}{dx^2}dx