text:
User:EML4500.f08.JAMAMA/FE

= 7.1  and  =

Given: 2D Strong Form and BCs

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

$$\Omega=\bar{w}=\Box$$ biunit square $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.1)
 * }

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

$$div(\underline{K} \cdot grad (u))+f=\rho c \frac{\partial u}{\partial t}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.2)
 * }

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

$$ \underline{K}= \underline{I}$$ (identity matrix) $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.3)
 * }

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

$$f=0 \quad$$,  $$\frac{\partial u}{\partial t}= 0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.4)
 * }

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

Essential boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * $$g=2 \quad$$ on  $$\partial\Omega \quad$$

$$
 *  $$ \displaystyle                                                                          (Eq. 7.1.5)
 * }

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

Natural boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * none

$$
 *  $$ \displaystyle\!
 * }

Find: Temperature inside and error at (0,0)

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

Solve for $$u^h$$ until error is less then $$10^{-6}$$ at center (x,y)=(0,0) for: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle\!
 * }

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

1) Static (steady state): $$f(x)=1 \quad$$ in $$\Omega=\Box$$
 * style="width:95%" |
 * style="width:95%" |

 * 1a) 2d LIBF

 * 1b) 2d LLEBF

 * 1b1) Uniform mesh

 * 1b2) Non-uniform mesh

2) Dynamic(transient):$$\rho c=3 \quad$$

 * Initial condition $$u(x,t=0)=xy \qquad \forall x \in \Omega$$

$$
 *  $$ \displaystyle\!
 * }

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

2a)
 * style="width:95%" |
 * style="width:95%" |

 * $$f(x,t)=0 \qquad \forall x \in \Omega,\ \forall t >0$$

 * $$g(x,t)=2 \mbox{ on }\Gamma_g=\partial \Omega$$

 * 2d LIBF

$$
 *
 $$ \displaystyle\!
 * }

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

2b)
 * style="width:95%" |
 * style="width:95%" |

 * $$f(x,t)=1 \qquad \forall x \in \Omega,\  \forall t >0$$

 * $$g(x,t)=2 \mbox{ on } \Gamma_g=\partial \Omega$$

 * 2b1) 2d LIBF

 * 2b2) 2d LLEBF

 * 2b2a) Uniform mesh

 * 2b2b) non-uniform mesh

$$
 *
 $$ \displaystyle\!
 * }

=Solution=

In general by Prof. Vu Quoc's

$$ \tilde{m}(w,u) + \tilde{k}(w,u) = \tilde{f}(w,u) \ $$

$$ \tilde{m}(w,u) = \int_\Omega w \rho c \frac{\partial u}{\partial t} d \Omega \ $$

$$ \tilde{k}(w,u) = \int_{\Omega^e} \triangledown{w} \kappa(x, y) \triangledown{u} d \Omega \ $$

$$ \tilde{f}(w) = - \int_{\Gamma_h}w h d \Gamma_h + \int_\Omega w f d \Omega \ $$

$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

, where $$ x = x^e \ $$, where $$x^e \ $$ represents local node.

For linear 2D case with boundary being square

Bases, $$ b_i \ $$ are defined to be $$ L_{[I,J]} (x, y) = N^e_{[I,J]} (x, y) \ $$

where $$ N^e_{[I,J]}(x, y) = N_I^e(x)N_J^e(y) \ $$ and $$ I, J \ $$ corresponds to numbering of nodes on element scale.

also noting that it's a rectangular shape, such that the distance from x: 1-to-2 = 1-to-3 and node 1 = node 4 in 'x' coordiantes, y: 1-to-4 = 1-to-3 and node 1 = node 2 in 'x' coordinates, we substitute $$ y^e_4 $$ for $$ y^e_3 \ $$ and $$ y^e_1 $$ for $$ y^e_2 \ $$ and the same for other coordinates

$$ N^e_{[1,1]} (x, y) = \frac{x^e-x^e_2}{x^e_1 - x^e_2} \frac{y^e-y^e_4}{y^e_1 - y^e_4} = \frac{1}{4}(x-\xi)(y - \eta) = b_1^e (x,y) \ $$

$$ N^e_{[2,1]} (x, y) = \frac{x^e-x^e_1}{x^e_2 - x^e_1} \frac{y^e-y^e_4}{y^e_1 - y^e_4} = \frac{1}{4}(x-\xi)(y - \eta)  = b_2^e(x,y)\ $$

$$ N^e_{[2,2]} (x, y) = \frac{x^e-x^e_1}{x^e_2 - x^e_1} \frac{y^e-y^e_1}{y^e_4 - y^e_1} = \frac{1}{4}(x-\xi)(y - \eta)  = b_3^e (x,y) \ $$

$$ N^e_{[1,2]} (x, y) = \frac{x^e-x^e_2}{x^e_1 - x^e_2} \frac{y^e-y^e_1}{y^e_4 - y^e_1} = \frac{1}{4}(x-\xi)(y - \eta)  = b_4^e (x,y) \ $$

$$ \mathbf{BN} =

\begin{bmatrix} \frac{\partial N^e_{[1,1]} (x, y) }{\xi} & \frac {\partial N^e_{[2,1]} (x, y) }{\xi} & \frac {\partial N^e_{[2,2]} (x, y) }{\xi} & \frac{\partial N^e_{[1,2]} (x, y) }{\xi}   \\

\frac{\partial N^e_{[1,1]} (x, y) }{\eta} & \frac {\partial N^e_{[2,1]} (x, y) }{\eta} & \frac {\partial N^e_{[2,2]} (x, y) }{\eta} & \frac{\partial N^e_{[1,2]} (x, y) }{\eta}

\end{bmatrix}

\ $$

$$ [\mathbf{x} \ \mathbf{y}] \ $$

Jacobian Matrix = BN *[mathbf{x} \ mathbf{y}]

taking the $$ x^e \ $$ and $$ y^e \ $$ values from the definition of parent element and substituting into Jacobian matrix we get:

$$ \mathbf{J}^e =

\begin{bmatrix} \frac{\partial x}{\partial \xi} &  \frac{\partial y}{\partial \xi}    \\ \\ \frac{\partial x}{\partial \eta} &  \frac{\partial y}{\partial \eta} \end{bmatrix} =

\begin{bmatrix} \frac{1}{2} (x_2^e - x_1^e) & 0    \\ \\ 0 &  \frac{1}{2} (y_4^e - y_1^e) \end{bmatrix}

$$

From given

$$ d \omega_x = det(\mathbf{J}^e) d \omega_\xi = d \omega_x = \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \ $$

or

$$ d \omega_\xi = \frac{ d \omega_x}{det(\mathbf{J}^e)} = d \omega_\xi = \frac{ d \omega_x}{\frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e)} $$

With this last information we have a completely defined stiffness element matrix $$ K^e \ $$. Also note that since we have, 4 degrees of freedom in 1 element node, we do expect the gram matrix $$ \Gamma \ $$ to be 4x4, therefore, $$ \kappa_{4x4} \ $$

$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

$$ K_{ij}^e = \int_{-1}^{+1}

\begin{bmatrix} \triangledown{b^e_1} \triangledown{b^e_1} & \triangledown{b^e_1} \triangledown{b^e_2} & \triangledown{b^e_1} \triangledown{b^e_3} & \triangledown{b^e_1} \triangledown{b^e_4}      \\ & \triangledown{b^e_2} \triangledown{b^e_2} & \triangledown{b^e_2} \triangledown{b^e_3} & \triangledown{b^e_2}\triangledown{b^e_4}      \\ & &  \triangledown{b^e_3} \triangledown{b^e_3} & \triangledown{b^e_3} \triangledown{b^e_4}      \\ symetry & &   & \triangledown{b^e_4} \triangledown{b^e_4} \end{bmatrix}

\begin{bmatrix} 1 & 0 & 0 & 0      \\ 0 & 1 &  0 & 0      \\ 0 & 0 &  1 & 0      \\ 0 & 0 &  0 & 1      \end{bmatrix}

\left ( d \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \right ) \ $$

note, to save some space in typing the formula (else it stretched out of the page), the independent variables for $$ b_i^e (\xi, \eta) \ $$ were avoided

for the initial square of length of '1', we divided into 4 smaller squares each 0.5 length, the element nodes stiffness matrix came out to be:

K(:,:,1) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24] K(:,:,2) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24]

this answer shows 2 elements stiffness matrices located along 'x' axes. Noting that because of the symmetry of each elements, the equivalent spacing of each node in elements, and because the material matrix $$ \kappa \ $$ is identity matrix, all stiffness element matrices will be the same. The only differ in element stiffness matrix will occur when we will change the number of stiffness.

,where

$$ b^e_1(\xi, \eta) = \frac{1}{A^e}(x_1^e \frac{1-\xi}{2} + [ \frac{1+\xi}{2}- 1 ] x^e_2)(y_1^e \frac{1-\eta}{2} + [ \frac{1+\eta}{2} - 1 ] ) y^e_4) \ $$

$$ b_2^e(\xi, \eta) =  -\frac{1}{A^e}( [\frac{1-\xi}{2}- 1]x^e_1 + x_2^e \frac{1+\xi}{2})(y_1^e \frac{1-\eta}{2} +  [\frac{1+\eta}{2} - 1] y^e_4)  \ $$

$$ b_3^e (\xi, \eta)  =  \frac{1}{A^e}( [\frac{1-\xi}{2}- 1]x^e_1 + x_2^e \frac{1+\xi}{2})([\frac{1-\eta}{2} -1]y^e_1 + y_4^e \frac{1+\eta}{2} ) \ $$

$$  b_4^e (\xi, \eta) =  -\frac{1}{A^e}(x_1^e \frac{1-\xi}{2} + [\frac{1+\xi}{2}- 1]x^e_2)( [\frac{1-\eta}{2}-1]y^e_1 + y_4^e \frac{1+\eta}{2})  \ $$

force
From prof. Vu Quoc's notes for

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega wf d \omega \ $$

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega w1 d \omega \ $$

converting to local coordinates

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) =

- \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

F(:,:,1) = [ 1/16,  1/16] [ -1/16,  1/16] [ -1/16, -1/16] [  1/16, -1/16] F(:,:,2) = [ 1/16,  1/16] [ -1/16,  1/16] [ -1/16, -1/16] [  1/16, -1/16]

Global Matrix creation
Then we constructed the Global stiffness matrix by the process of superimposing the local stiffness coefficients. Matlab output below shows the process of superimposing the elements to global stiffness matrix.

The matlab code for global matrix construction was written with respect to this element naming system, and then it was extended to encompass other numbers of elements

Original element matrices

K(:,:,1) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24]

K(:,:,2) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24] Superpositions to global matrix

KGlobal = 1/24         -1/96           0             -1/96          -1/48           0              0              0              0             -1/96           1/12          -1/96          -1/48          -1/48          -1/48           0              0              0              0             -1/96           1/24           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/12          -1/48           0             -1/96          -1/48           0             -1/48          -1/48          -1/48          -1/48           1/6           -1/48          -1/48          -1/48          -1/48           0             -1/48          -1/96           0             -1/48           1/12           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/24          -1/96           0              0              0              0             -1/48          -1/48          -1/48          -1/96           1/12          -1/96           0              0              0              0             -1/48          -1/96           0             -1/96           1/24

FGlobal = 1/8           0             -1/8            1/4            0             -1/4            1/8            0             -1/8

Boundary Conditions
Noting that essential BCs were specified to be on all sides of the 2x2 square. It would mean (viewing from the above global node figure) all nodes around node 5 will be tagged as essential BC.

Since we can either eliminate the nodes form Global Stiffness matrix to construct another one with Temperature ('u' function) as unknown or just '0' the node values in Global Stiffness matrix and then do inverse and multiply by Force matrix to find other non-constrained dofs.

We wrote a matlab code that automatically constrains the associate nodes and puts their contribution to force function.

Here how code applies BC's

KGlobal = 1/24         -1/96           0             -1/96          -1/48           0              0              0              0             -1/96           1/12          -1/96          -1/48          -1/48          -1/48           0              0              0              0             -1/96           1/24           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/12          -1/48           0             -1/96          -1/48           0             -1/48          -1/48          -1/48          -1/48           1/6           -1/48          -1/48          -1/48          -1/48           0             -1/48          -1/96           0             -1/48           1/12           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/24          -1/96           0              0              0              0             -1/48          -1/48          -1/48          -1/96           1/12          -1/96           0              0              0              0             -1/48          -1/96           0             -1/96           1/24

Kchange = 0             0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              1/6            0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0

FGlobal = -1/16          0              1/16          -1/8            0              1/8           -1/16           0              1/16

Fchange = -7/48         -7/48           0             -1/6            1/3            1/12          -7/48          -7/48           0

Then the matrices were reduced to the ones of the importance ( Global K matrix after the BCs where applied was isolated so no 0s would be present and attributing F function value(s) where picked by code as well)

KchangeF = 1/6

FchangeF = 1/3

Finally Temperature at unknown nodes was found by using this formula

$$ u = K^{-1}F \ $$

u = 2

Also output for 4x4 element free node values came out to be:

dcalc = 2             2              2              2

as you may see, only 1 node, in this example was present, so only 1 numerical value was generated.

Discussion
Note 1 important aspect of this output. That is, with forcing function changed to 1, the output of free degrees of freedom is still '2' as if no forcing function was present to begin with. I guess it might be explainable that because of constrained boundary conditions at all 4 ends, in order to keep the temperature constant at all 4 ends while heat was generated, the generated heat had to be instantaneously dissipated though the boundaries of the area, so that boundary temperature would stay unchanged.

As in matrix/numeral explanation, it should be noted, that because essential boundary condition is applied over all 4 sides, it has integration limits (working boundary range) as all element boundary. Therefore, with intorduction with heat source, force instead of being modified, was only scaled (remember force equation):

for f = 0 (HW 6.6)

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

for f= 1 (HW 7.1)

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) =

- \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

you may see that force term is a scale only term. Moreover, notice the force global term

For 4 element at f = 0: FGlobal = 1/8           0             -1/8            1/4        0            -1/4            1/8            0             -1/8

Fchange = 5/24          7/48          -1/16           7/24         -1/3           -5/24           5/24           7/48          -1/16

For 4 element at f = 1:

FGlobal =

-1/16          0              1/16          -1/8            0             1/8           -1/16           0              1/16

Fchange =

-7/48         -7/48           0             -1/6           -1/3            1/12          -7/48          -7/48           0

For 9 element at f = 1

FGlobal =

-1/36          0              0              1/36          -1/18           0        0            1/18          -1/18           0        0            1/18          -1/36           0              0              1/36

Fchange =

-7/108        -7/108         -7/108          0            -11/108         5/54           5/54           1/108         -5/54         5/54           5/54           1/54          -7/108         -7/108         -7/108          0

, where boxed answer shows free nodes.

Note that force matrix DOES differ, but note that because of the similarity of force functions (scalar multiples of each other, namely, caused because forcing function is a constant and essential boundary condition acts on the same area as forcing function), the global force matrix, since it contained '0's on initial force function, will contain '0's on any other force function, as long as forcing function is scalar and essential B.C is the same as whole element B.C. Therefore, with same used stiffness matrix 'K' the adjusted force matrix with initial values of '0' on free nodes, always came out to be the same for specific number of elements with differing scalar value of forcing function.

As a result, for free node displacements, we observed the constant value of '2'.

--Eml5526.s11.team5.JA 20:51, 6 April 2011 (UTC)

Transient Term
$$ \tilde{m}(w,u) = \int_\Omega w \rho c \frac{\partial u}{\partial t} d \Omega = \int_\Omega w xy \frac{\partial u}{\partial t} d \Omega =

\int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( x_1^e \frac{1-\xi}{2} + x_2^e \frac{1+\xi}{2} \right ) \left( y_1^e \frac{1-\eta}{2} + y_4^e \frac{1+\eta}{2} \right ) \frac{\partial u}{\partial t} \left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

Overall Equation
$$ \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( x_1^e \frac{1-\xi}{2} + x_2^e \frac{1+\xi}{2} \right ) \left( y_1^e \frac{1-\eta}{2} + y_4^e \frac{1+\eta}{2} \right ) \frac{\partial u}{\partial t} \left( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1} \triangledown{b^e_1} & \triangledown{b^e_1} \triangledown{b^e_2} & \triangledown{b^e_1} \triangledown{b^e_3} & \triangledown{b^e_1} \triangledown{b^e_4}      \\ & \triangledown{b^e_2} \triangledown{b^e_2} & \triangledown{b^e_2} \triangledown{b^e_3} & \triangledown{b^e_2}\triangledown{b^e_4}      \\ & &  \triangledown{b^e_3} \triangledown{b^e_3} & \triangledown{b^e_3} \triangledown{b^e_4}      \\ symetry & &   & \triangledown{b^e_4} \triangledown{b^e_4} \end{bmatrix}

\begin{bmatrix} 1 & 0 & 0 & 0      \\ 0 & 1 &  0 & 0      \\ 0 & 0 &  1 & 0      \\ 0 & 0 &  0 & 1      \end{bmatrix}

\left ( d \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \right ) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

Static using 2d LLEBF non-uniform mesh

Here the answer is also 2 at every node, which agrees with the uniform mesh LLEBF. There is a slight difference between nodes which can be attributed to matlab roundoff error.

Dynamic (case 2b) using 2d LLEBF uniform mesh

Dynamic (case 2b) using 2d LLEBF non-uniform mesh

Eml5526.s11.team5.savery

Eml5526.s11.team5.smith 20:37, 21 April 2011 (UTC)

Eml5526.s11.team5.smith 21:11, 21 April 2011 (UTC)

= 7.2 Finding the static and dynamic solution of a bi-unit square using 2D LIBF =

P.D.E

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

$$T div (grad u) + f = \rho \ \frac{\partial^2u}{\partial t^2} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.1)
 * }

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

$$ \quad T = 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.2)
 * }

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

$$ f (\mathbf{x}) = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.3)
 * }

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

$$ \quad \rho = 3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.4)
 * }

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

$$ \quad Basis \ Function \ : \ 2D \ Lagrange \ Interpolation \ Basis \ Function $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Essential Boundary Condition

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

$$ \quad g = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.5)
 * }

Static(steady state)

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

1. Find the static solution $$ \quad U^h $$ to $$ \quad 10^{-6} $$ accuracy at the center. $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Dynamic(transient)

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

1. Fundamental Frequency and Fundamental Period $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

2. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1] $$ for symmetric and Non-Symmetric $$ \quad u^0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

3. Produce a movie of the vibrating membrane for symmetric and Non-Symmetric $$ \quad u^0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Solution
Substituting the given values in the P.D.E, we get,
 * {| style="width:100%" border="0"

$$\quad 4 \ div (grad u) \ + \ 1 \ = \ 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$ 4 \left ( \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2}\right ) + 1 = 0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$ \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} = \frac{-1}{4} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.6)
 * }

2D LIBF

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

$$ N_I (x,y) = L_{i,m}(x).L_{j,n}(y) \quad Where \ I = i+(j-1)m $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.7)
 * }

Where,

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

$$ L_{i,m}(x) = \prod_{k=1\neq i}^{m}\frac{x-x_k}{x_i - x_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }
 * {| style="width:100%" border="0"

$$ L_{j,n}(x) = \prod_{k=1\neq j}^{m}\frac{y-y_k}{y_j - y_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Consider the Weak Form

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

$$  \tilde{m}(w,u) + \tilde{k}(w,u) = \tilde{f}(w) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.8)
 * }

where,

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

$$m(w,u) = \int\limits_\Omega {w\rho c\frac} d\Omega $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$ k(w,u) = \int\limits_\Omega {\nabla wK\nabla u} d\Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$ f(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Static case
For static case, $$ \frac{\partial u}{\partial t} = 0 $$ the above equation can be written as,

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

$$ \quad K d = F $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

where,

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

$$ {K_{ij}} = \int_ {(\nabla N_I^e)I(\nabla N_J^e)d{w^e}} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.9)
 * }

Using the above equations in MATLAB, we get,

K =  33.3257  -13.5386  -13.5386    2.8639 -13.5386  33.3257    2.8639  -13.5386       -13.5386    2.8639   33.3257  -13.5386         2.8639  -13.5386  -13.5386   33.3257

F = [0.5625  0.5625  0.5625  0.5625]^T

d = [0.0617  0.0617  0.0617  0.0617]^T

Dynamic case
For this case, It is given that the value of mass density $$ \quad \rho \ = \ 3$$

Free Vibration: Solve Eigen Value problem
The equation given to solve the Eigen value problem is,

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

$$ \mathbf{K} \mathbf{\phi} = \lambda \mathbf{M}\mathbf{\phi}  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.10)
 * }

Where,

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

$$\quad \mathbf{\phi} - Eigen \ vectors $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$ \quad \lambda - Eigen \ values $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Using the above equation and the expression for K, F matrix, we can calculate the fundamental frequency and fundamental time period.

Fundamental Frequency = 0.23725

Fundamental Time Period = 4.2149

Transient Analysis with symmetric Uo
Initial conditions

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

$$u(\mathbf{x},t=0) = u^h_s(x) \ \forall \  \mathbf{x} \ \epsilon  \ \Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.11)
 * }

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

$$ \dot{u} (\mathbf{x},t=0) = 0 \ \forall \  \mathbf{x} \ \epsilon  \ \Omega $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.12)
 * }

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

$$T = 4, \qquad f(\mathbf{x},t) = 0 \ in \  \Omega  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.13)
 * }

By using the above initial condition, we can produce the movie of the vibrating membrane.

Movie of the vibrating membrane
The Movie of dynamic membrane has been uploaded in the youtube and it can be seen from the link below

FE Analysis of a Dynamic Membrane - symmetric case

For T = 0,

Transient Analysis with Non-symmetric Uo
In this case, initial u value is given by,

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

$$u(\mathbf{x},t=0) = (x+1)(y+\frac{1}{2}) cos (\frac{\pi}{2}x) cos (\frac{\pi}{2}y) \ \forall \  \mathbf{x} \ \epsilon  \ \Omega$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle                                                                     (Eq. 7.2.14)
 * }

By using the above initial condition, we can produce the movie of the vibrating membrane.

Movie of the vibrating membrane
The Movie of dynamic membrane for the non-symmetric case, has been uploaded in the youtube and it can be seen from the link below

FE Analysis of a Dynamic membrane - Non Symmetric case

--Eml5526.s11.team5.srv 21:06, 21 April 2011 (UTC)

MATLAB code for producing the movie of the vibrating membrane
--Lokeshdahiya 21:09, 21 April 2011 (UTC)

= 7.3 =

Given/Find: Static Solution of the Circle

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

Let Omega a circle with a one-unit radius. Find the static solution such that: $$T=4$$ $$f(\underline x)=1$$ in Omega $$g=0$$ on $$\Gamma _g=\partial \Omega $$
 * style="width:95%" |
 * style="width:95%" |

all to 10^-6 accuracy at the center. Quadrilateral and Triangular elements should be used but for only 1/4 of the circle. Plot the deformed shape in 3D. $$
 *
 $$ \displaystyle\!
 * }

Solution: Static Solution of the Circle Using Quadrilateral and Triangular Elements to an Accuracy of 10^-6

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

The Partial Differential Equation for a vibrating membrane is the following via transient analysis: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

T(\textrm{div \ grad \ u})+\textrm{f}=\rho \frac{\partial ^2u}{\partial t^2} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

Where: $$ u(\underline x,t)=$$ Transverse Displacement $$T=\quad$$ Tension (Force/Length) = Constant $$\textrm{f}(\underline x,t)=$$ Distributed Transverse Force $$\rho (\underline x)=$$ Mass Density (Mass/Area) & $$\textrm{div \ grad\ u}=\frac{\partial }{\partial x_i}\frac{\partial u}{\partial x_j}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

The Weighted Residual Form is as follows: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\int_{\Omega } w\left [ T\frac{\partial ^2u}{\partial x_i\partial x_j}+\textrm{f}-\rho \frac{\partial ^2u}{\partial t^2} \right ]d\Omega =0 $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

Which can thereby be broken into three separate terms: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\int_{\Omega } wT\frac{\partial ^2u}{\partial x_i\partial x_j}d\Omega =\mathbb{K} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\int_{\Omega } w\textrm{f}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\int_{\Omega } w\rho \frac{\partial ^2u}{\partial t^2}\ d\Omega=\tilde{m}(w,u) $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

This final equation will be used later in the Weak Form. For now, the first term above (set equal to $$\mathbb{K}$$), is equivalent to: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\mathbb{K}=\int_{d\Omega } \frac{\partial }{\partial x_i}(wT\frac{\partial u}{\partial x_j})\ d\Omega-\int_{\Omega } \frac{\partial w}{\partial x_i}T\frac{\partial u}{\partial x_j}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

And by using divergence theorem, the first term of the equation above may be converted to yield: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\mathbb{K}=\int_{\partial \Omega } n_iwT\frac{\partial u}{\partial x_j}d(\partial \Omega)-\int_{\Omega } \frac{\partial w}{\partial x_i}T\frac{\partial u}{\partial x_j}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

Now since $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\partial \Omega =\Gamma _g\cup \Gamma _h\ \textrm{but }\ w=0\ \ \textrm{on}\ \ \Gamma _g $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

and the following are true: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\frac{\partial w}{\partial x_i}=\textrm{grad}\ w=\bigtriangledown w $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\frac{\partial u}{\partial x_j}=\textrm{grad}\ u=\bigtriangledown u $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

Then the K equation becomes: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\mathbb{K}=\int_{\Gamma _h}(n_i\ wT\frac{\partial u}{\partial x_j}\ \partial \Gamma _h)-\int_{\Omega }\bigtriangledown wT\bigtriangledown u\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

Otherwise written as: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\mathbb{K}=\int_{\Gamma _h}w\ \underline nT\cdot \textrm{grad \ u}\ \partial \Gamma _h-T\int_{\Omega }\bigtriangledown w\bigtriangledown u\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

where the last term is conveniently used to write ultimately the weak form of the weighted residual form which is the following: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

\tilde{\textrm{f}}(w)=\tilde{m}(w,u)+\tilde{k}(w,u)\ \ \ \forall w \ \ \ s.t.\ \ \ w=0\ \ \ on\ \ \ \Gamma _g $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\tilde{m}(w,u)=\int_{\Omega } w\rho \frac{\partial ^2u}{\partial t^2}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\tilde{k}(w,u)=T\int_{\Omega }\bigtriangledown w\bigtriangledown u \ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

\tilde{\textrm{f}}(w)=\int_{\Gamma _h}w \underline n \ T\ \textrm{grad u}\ \partial \Gamma _h+\int_{\Omega }w\textrm{f}\ d\Omega $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

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

From the Discrete Weighted Form, $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$\underline M_{FF}\underline d_F+\underline K_{FF}\underline d_F=\underline F_F-(\underline M_{FE}g+\underline K_{FE}g)$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

But because g=0 from the boundary conditions, as well as M due to the static requirement, the displacements may be derived: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$\underline K \underline d=\underline F$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$\therefore \underline d=\underline K^{-1}\underline F$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

Now to solve for the displacements: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$N_I(x,y)=L_{i,m}(x)\cdot L_{j,n}(y) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

Where $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$I=i+(j-1)m $$ and $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$L_{i,m}(x)=\prod_{K=1, K\neq i}^{m}\left ( \frac{x-x_k}{x_i-x_k} \right ) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

And setting the Shape Function to the following: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$N=\begin{bmatrix} \frac{1}{4}(1-\xi _1)(1-\xi _2)\\ \frac{1}{4}(1+\xi _1)(1-\xi _2)\\ \frac{1}{4}(1-\xi _1)(1+\xi _2)\\ \frac{1}{4}(1-\xi _1)(1+\xi _2) \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

The deformed shape coordinates may further be derived from these following two equations, where the summations are taken across all coordinate points taken from Abaqus output files: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$X_{1}^{e}=\varphi_{1}^{e}(\xi )=\sum_{I=1}^{nel}N_I(\xi )x_{1,I}^{e} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$(y=)X_{2}^{e}=\varphi_{2}^{e}(\xi )=\sum_{I=1}^{nel}N_I(\xi )x_{2,I}^{e} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

Using the following Jacobian matrix: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$J=\left [ \frac{\partial x_i}{\partial \xi _j} \right ]=\begin{bmatrix} \frac{\partial x_1}{\partial \xi _1} & \frac{\partial x_1}{\partial \xi _2}\\ \frac{\partial x_2}{\partial \xi _1} & \frac{\partial x_2}{\partial \xi _2}\\ \end{bmatrix} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

The following equations may now be utilized: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$\underline B_I(\xi )=\bigtriangledown _xN_{I}^{e}(\xi )=(\underline J^e)^{-T}(\xi )\bigtriangledown_\xi N_{I}^{e}(\xi ) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$k_{I3}^{e}=\int_{\omega =el}N_I(\xi )\cdot T\cdot \underline B_J (\xi )J(\xi )d(el) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

$$F_i=\int_{\omega =el}N_I(\xi )f\ det(J(\xi ))d(el) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

Implementing MATLAB at this point, in conjunction with output data files from Abaqus which contain element coordinate points for three separate quadrilateral meshes and three separate triangular meshes (included), the final displacements may be derived, to an accuracy of 10^-6: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

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

The above MATLAB script for the 6 different meshes generate the following results; the first three show quadrilateral elements in order of increasing elements. The last three show triangular elements in order of increasing elements. $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

The figure to the left shows a unit-less displacement of: 0.0623258

The figure to the left shows a unit-less displacement of: 0.0624613

The figure to the left shows a unit-less displacement of: 0.0625011

The figure to the left shows a unit-less displacement of: 0.0624613

The figure to the left shows a unit-less displacement of: 0.0639197

The figure to the left shows a unit-less displacement of: 0.0619208

The figure to the left shows a unit-less displacement of: 0.0621967

Comparision With ABAQUS
The problem can be converted to an equivalent Heat problem with the following assumptions:

Conductance K=4,

Temperature T on the boundary = 0

Heat generation f=1

The equivalent model was analyzed in ABAQUS and the following result was obtained

The Temperature at the center turned out to be 6.25e-02 which is equal to what we got using MATLAB.

--Eml5526.s11.team5.vijay 21:17, 21 April 2011 (UTC)

--Eml5526.s11.team5.mcdaniel 21:23, 21 April 2011 (UTC)

= 7.4 : Newton's Law of Cooling http://en.wikiversity.org/w/index.php?title=File:Fe1.s11.mtg41.djvu&page=5 =

Given
Data Set G2DM2.0/D1

$$ \Omega = Domain \quad $$

$$ K=I, f(x,t)=0$$

$$ \frac {\partial u}{\partial t} = 0, ~g=2, ~h=3$$

$$ H=.5 ~u_\infty =.1$$

$$ \partial \Omega = \Gamma_g \bigcup \Gamma_h \bigcup \Gamma_H $$

and, The Partial Differential Equation is:

Find
1) The weak form for heat conduction in 2-D with boundary convection. 2) Develop semi-discrete equations (ODEs). Give detailed expression of all quantities. 3) Solve G2DM2.0/D1 using 2D LIBF till $$\displaystyle 10^{-6} $$ accuracy at the center. Plot solution in 3D with contour.

1.The Weak Form for heat conduction in 2D with boundary convection
From the PDE (Eqn 6.10.1) the weak form can be written as,

Using integration by parts,

Substituting into the integral we obtain,

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

 * style="width:95%" |$$
 * style="width:95%" |$$

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial }{\partial {{x}_{i}}}\left( w{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0 $$ $$
 *
 $$ \displaystyle (Eq. 7.10.3)
 * }

Applying the Gauss Theorem on the first term gives,

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}w{{K}_{ij}}\frac{\partial u}{\partial{{x}_{j}}}{{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0 $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle (Eq. 7.10.4)
 * }

The Heat Flux is given by,

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

q_{i}=-K_{ij}\frac{\partial u}{\partial x_{j}} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

Substituting the heat flux and rearranging the terms in Eq. 7.10.4 gives,

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{-w{{q}_{i}}{n}_{i}}d\text{ }\!\!\Gamma\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)}d\text{ }\!\!\Omega\!\!\text{ +}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{wf\left( x,t \right)}d\text{ }\!\!\Omega\!\!\text{ }-\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}}d\text{ }\!\!\Omega\!\!\text{ }=0 $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle (Eq. 7.10.5)
 * }

The first term can be rewritten as,

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

\int\limits_{\Gamma }{-w{{q}_{i}}{{n}_{i}}d\Gamma =}\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{g}}}-\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{h}}}-\int\limits_{-w{{q}_{i}}{{n}_{i}}d{{\Gamma }_{H}}} $$ $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle\!
 * }

Select $$ \quad w $$ such that $$\quad w=0 $$ on $$\quad {{\Gamma }_{g}}$$

Organizing, Eq. 7.10.5 then becomes the following weak form

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ +}}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\frac{\partial w}{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)d\text{ }\!\!\Omega\!\!\text{ }}+\int\limits_{wHu}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}=-\int\limits_{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}+\int\limits_{wH{{u}_{\infty }}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}+\int\limits_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }},$$ $$ $$ \quad \forall w $$ such that $$\displaystyle w=0 \quad on \quad {{\Gamma }_{g}} $$
 * style="width:95%" |$$
 * style="width:95%" |$$
 *
 $$ \displaystyle (Eq.7.10.6)
 * }

2) Develop Semidiscrete Equations (ODE's)
Let $$ u=[ N]{{u}_{e}}\quad $$ and $$ \quad w=[ N ]{{w}_{e}},$$ where $$\displaystyle [N] $$ is the Lagrange polynomial and $${u}_{e}~ and ~ {w}_{e} $$ are the nodal values of $$ u \quad $$ and $$w \quad $$ respectively. Now, take the derivative of basis functions with respect to x1 and x2

Substituting the Langrangian form into the weak form (Eq 7.10.6), yields the semi-dicrete form:

From Eq. 6.10.8, we can obtain the Heat Source vector, Conductivity matrix, and Capacitance Matrix.

The heat source vector is $$\tilde{f}(w)= -\int\limits_{h{{[N]}^{T}}d{{\Gamma }_{h}}}+\int\limits_{H{{[N]}^{T}}{{u}_{\infty }}d{{\Gamma }_{h}}}+\int\limits_{\Omega }{{{[N]}^{T}}f(x,t)d\Omega } $$

The Conducitivity matrix is $$\tilde{k}(w,u)= \int\limits_{\Omega }{{{[B]}^{T}}[B]d\Omega }+\int\limits_{H{{[N]}^{T}}[N]d{{\Gamma }_{h}}} $$

The Capacitance Matrix is $$\tilde{m}(w,u)= \int\limits_{\Omega }{\rho c{{[N]}^{T}}[N]\frac{\partial u}{\partial t}d\Omega } $$

3.Solve G2DM2.0/D1 Using 2D LIBF till $$10^{-6}$$ accuracy at center. Plot Solution in 3D w/Contour.
To solve the problem using LIBF, divide the x-axis and Y-axis into 2 parts 2D LIBF

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

$$ N_I (x,y) = L_{i,m}(x).L_{j,n}(y) \quad Where \ I = i+(j-1)m $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!                                                                     $$
 * }

Where,

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

$$ L_{i,m}(x) = \prod_{k=1\neq i}^{m}\frac{x-x_k}{x_i - x_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }
 * {| style="width:100%" border="0"

$$ L_{j,n}(x) = \prod_{k=1\neq j}^{m}\frac{y-y_k}{y_j - y_k} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle\!
 * }

Consider a case of 9 nodes as shown below

The trial solution is

The trial solution must satisfy the essential boundary condition. Throughout boundary $$ \Gamma_g, \quad $$ the temperature should be a constant 2. Thus $$ d_1, d_2, d_3, d_6, d_9 =2. \quad $$ --Eml5526.s11.team5.smith 21:22, 21 April 2011 (UTC)

=Contributing Members =

title: User:EML4500.f08.JAMAMA/FE


text:
Template:TopicsWithSubpages

Using Template:TopicsWithSubpages
This template was transcluded from wikieducator - Template:TopicsWithSubpages

This Navigational template may be useful for complex content that has a hierarchical structure to it such as:

Topic 1

section 1

section 2

Topic 2

section 1

section 2

When inserted into a page it might look something like this:

(For a simpler template suitable for Wikiversity projects see: Template:EasyNavBar.)

Variables
The template has the following variables:

state = can be either left blank, or collapsed or plain can be entered here. Blank (ie. nothing entered) is the default and will provide the user with a way to collapse the navigation template with a show/hide link; collapsed will default the navigation template into the collapsed or closed state; and plain removes the show/hide option all together.

title = This is the title of the resource. Can be left blank.

image = You can include an image or logo to display in the template. Uses normal wiki image syntax.

above = Text to go in the header cell of the template. In the above example the text: "Test of a Navigation Template with Topics and Sections" is here.

topic1 = this would be a link to your first Topic.

section1 = this may be series of links separated by of the sections for topic 1.

topic2 = to add new topics increment the number; ie. topic2, topic3.

section2 = sections/pages for the topic2 table cell.

below = text that could go at the bottom of the navigation box.

Automated Creation of TopicWithSubpage Navigation for projects
Instead of copying and pasting all the template syntax into every page of your content, it may be easier to use this form to create a unique navigation template for your all of your content pages. You can then insert this newly created navigation template into each page of your content using a greatly simplified syntax. This also makes it far easier to modify the navigation in a single place (ie, add or subtract pages, etc.) and have changes reflected in all pages rather than having to modify the syntax in every page.

Usage
Replace the MyProjectNav text below with the name of your project, something like Template:AccountingBasicsNav and click on the Create New Navigation Template button. You will be taken to a page with further instructions.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: Template:TopicsWithSubpages


text:
User:EML4500.f08.JAMAMA/FE/JA

Given: 2D Strong Form and BCs

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

$$\Omega=\bar{w}=\Box$$ biunit square $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.1)
 * }

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

$$div(\underline{K} \cdot grad (u))+f=\rho c \frac{\partial u}{\partial t}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.2)
 * }

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

$$ \underline{K}= \underline{I}$$ (identity matrix) $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.3)
 * }

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

$$f=0 \quad$$,  $$\frac{\partial u}{\partial t}= 0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 7.1.4)
 * }

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

Essential boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * $$g=2 \quad$$ on  $$\partial\Omega \quad$$

$$
 *  $$ \displaystyle                                                                          (Eq. 7.1.5)
 * }

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

Natural boundary condition:
 * style="width:95%" |
 * style="width:95%" |
 * none

$$
 *  $$ \displaystyle
 * }

Find: Temperature inside and error at (0,0)

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

Solve for $$u^h$$ until error is less then $$10^{-6}$$ at center (x,y)=(0,0) for: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

1) Static (steady state): $$f(x)=1 \quad$$ in $$\Omega=\Box$$
 * style="width:95%" |
 * style="width:95%" |

 * 1a) 2d LIBF

 * 1b) 2d LLEBF

 * 1b1) Uniform mesh

 * 1b2) Non-uniform mesh

2) Dynamic(transient):$$\rho c=3 \quad$$

 * Initial condition $$u(x,t=0)=xy \qquad \forall x \in \Omega$$

$$
 *  $$ \displaystyle
 * }

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

2a)
 * style="width:95%" |
 * style="width:95%" |

 * $$f(x,t)=0 \qquad \forall x \in \Omega,\ \forall t >0$$

 * $$g(x,t)=2 \mbox{ on }\Gamma_g=\partial \Omega$$

 * 2d LIBF

$$
 *
 $$ \displaystyle
 * }

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

2b)
 * style="width:95%" |
 * style="width:95%" |

 * $$f(x,t)=1 \qquad \forall x \in \Omega,\  \forall t >0$$

 * $$g(x,t)=2 \mbox{ on } \Gamma_g=\partial \Omega$$

 * 2b1) 2d LIBF

 * 2b2) 2d LLEBF

 * 2b2a) Uniform mesh

 * 2b2b) non-uniform mesh

$$
 *
 $$ \displaystyle
 * }

=Solution=

In general by Prof. Vu Quoc's

$$ \tilde{m}(w,u) + \tilde{k}(w,u) = \tilde{f}(w,u) \ $$

$$ \tilde{m}(w,u) = \int_\Omega w \rho c \frac{\partial u}{\partial t} d \Omega \ $$

$$ \tilde{k}(w,u) = \int_{\Omega^e} \triangledown{w} \kappa(x, y) \triangledown{u} d \Omega \ $$

$$ \tilde{f}(w) = - \int_{\Gamma_h}w h d \Gamma_h + \int_\Omega w f d \Omega \ $$

$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

, where $$ x = x^e \ $$, where $$x^e \ $$ represents local node.

For linear 2D case with boundary being square

Bases, $$ b_i \ $$ are defined to be $$ L_{[I,J]} (x, y) = N^e_{[I,J]} (x, y) \ $$

where $$ N^e_{[I,J]}(x, y) = N_I^e(x)N_J^e(y) \ $$ and $$ I, J \ $$ corresponds to numbering of nodes on element scale.

also noting that it's a rectangular shape, such that the distance from x: 1-to-2 = 1-to-3 and node 1 = node 4 in 'x' coordiantes, y: 1-to-4 = 1-to-3 and node 1 = node 2 in 'x' coordinates, we substitute $$ y^e_4 $$ for $$ y^e_3 \ $$ and $$ y^e_1 $$ for $$ y^e_2 \ $$ and the same for other coordinates

$$ N^e_{[1,1]} (x, y) = \frac{x^e-x^e_2}{x^e_1 - x^e_2} \frac{y^e-y^e_4}{y^e_1 - y^e_4} = \frac{1}{4}(x-\xi)(y - \eta) = b_1^e (x,y) \ $$

$$ N^e_{[2,1]} (x, y) = \frac{x^e-x^e_1}{x^e_2 - x^e_1} \frac{y^e-y^e_4}{y^e_1 - y^e_4} = \frac{1}{4}(x-\xi)(y - \eta)  = b_2^e(x,y)\ $$

$$ N^e_{[2,2]} (x, y) = \frac{x^e-x^e_1}{x^e_2 - x^e_1} \frac{y^e-y^e_1}{y^e_4 - y^e_1} = \frac{1}{4}(x-\xi)(y - \eta)  = b_3^e (x,y) \ $$

$$ N^e_{[1,2]} (x, y) = \frac{x^e-x^e_2}{x^e_1 - x^e_2} \frac{y^e-y^e_1}{y^e_4 - y^e_1} = \frac{1}{4}(x-\xi)(y - \eta)  = b_4^e (x,y) \ $$

$$ \mathbf{BN} =

\begin{bmatrix} \frac{\partial N^e_{[1,1]} (x, y) }{\xi} & \frac {\partial N^e_{[2,1]} (x, y) }{\xi} & \frac {\partial N^e_{[2,2]} (x, y) }{\xi} & \frac{\partial N^e_{[1,2]} (x, y) }{\xi}   \\

\frac{\partial N^e_{[1,1]} (x, y) }{\eta} & \frac {\partial N^e_{[2,1]} (x, y) }{\eta} & \frac {\partial N^e_{[2,2]} (x, y) }{\eta} & \frac{\partial N^e_{[1,2]} (x, y) }{\eta}

\end{bmatrix}

\ $$

$$ [\mathbf{x} \ \mathbf{y}] \ $$

Jacobian Matrix = BN *[mathbf{x} \ mathbf{y}]

taking the $$ x^e \ $$ and $$ y^e \ $$ values from the definition of parent element and substituting into Jacobian matrix we get:

$$ \mathbf{J}^e =

\begin{bmatrix} \frac{\partial x}{\partial \xi} &  \frac{\partial y}{\partial \xi}    \\ \\ \frac{\partial x}{\partial \eta} &  \frac{\partial y}{\partial \eta} \end{bmatrix} =

\begin{bmatrix} \frac{1}{2} (x_2^e - x_1^e) & 0    \\ \\ 0 &  \frac{1}{2} (y_4^e - y_1^e) \end{bmatrix}

$$

From given

$$ d \omega_x = det(\mathbf{J}^e) d \omega_\xi = d \omega_x = \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \ $$

or

$$ d \omega_\xi = \frac{ d \omega_x}{det(\mathbf{J}^e)} = d \omega_\xi = \frac{ d \omega_x}{\frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e)} $$

With this last information we have a completely defined stiffness element matrix $$ K^e \ $$. Also note that since we have, 4 degrees of freedom in 1 element node, we do expect the gram matrix $$ \Gamma \ $$ to be 4x4, therefore, $$ \kappa_{4x4} \ $$

$$ K_{ij}^e = \int_{\Omega^e} \triangledown{b^e_i}'(x, y) \kappa(x, y) \triangledown{b^e_j}'(x, y) d \omega_x \ $$

$$ K_{ij}^e = \int_{-1}^{+1}

\begin{bmatrix} \triangledown{b^e_1} \triangledown{b^e_1} & \triangledown{b^e_1} \triangledown{b^e_2} & \triangledown{b^e_1} \triangledown{b^e_3} & \triangledown{b^e_1} \triangledown{b^e_4}      \\ & \triangledown{b^e_2} \triangledown{b^e_2} & \triangledown{b^e_2} \triangledown{b^e_3} & \triangledown{b^e_2}\triangledown{b^e_4}      \\ & &  \triangledown{b^e_3} \triangledown{b^e_3} & \triangledown{b^e_3} \triangledown{b^e_4}      \\ symetry & &   & \triangledown{b^e_4} \triangledown{b^e_4} \end{bmatrix}

\begin{bmatrix} 1 & 0 & 0 & 0      \\ 0 & 1 &  0 & 0      \\ 0 & 0 &  1 & 0      \\ 0 & 0 &  0 & 1      \end{bmatrix}

\left ( d \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \right ) \ $$

note, to save some space in typing the formula (else it stretched out of the page), the independent variables for $$ b_i^e (\xi, \eta) \ $$ were avoided

for the initial square of length of '1', we divided into 4 smaller squares each 0.5 length, the element nodes stiffness matrix came out to be:

K(:,:,1) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24] K(:,:,2) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24]

this answer shows 2 elements stiffness matrices located along 'x' axes. Noting that because of the symmetry of each elements, the equivalent spacing of each node in elements, and because the material matrix $$ \kappa \ $$ is identity matrix, all stiffness element matrices will be the same. The only differ in element stiffness matrix will occur when we will change the number of stiffness.

,where

$$ b^e_1(\xi, \eta) = \frac{1}{A^e}(x_1^e \frac{1-\xi}{2} + [ \frac{1+\xi}{2}- 1 ] x^e_2)(y_1^e \frac{1-\eta}{2} + [ \frac{1+\eta}{2} - 1 ] ) y^e_4) \ $$

$$ b_2^e(\xi, \eta) =  -\frac{1}{A^e}( [\frac{1-\xi}{2}- 1]x^e_1 + x_2^e \frac{1+\xi}{2})(y_1^e \frac{1-\eta}{2} +  [\frac{1+\eta}{2} - 1] y^e_4)  \ $$

$$ b_3^e (\xi, \eta)  =  \frac{1}{A^e}( [\frac{1-\xi}{2}- 1]x^e_1 + x_2^e \frac{1+\xi}{2})([\frac{1-\eta}{2} -1]y^e_1 + y_4^e \frac{1+\eta}{2} ) \ $$

$$  b_4^e (\xi, \eta) =  -\frac{1}{A^e}(x_1^e \frac{1-\xi}{2} + [\frac{1+\xi}{2}- 1]x^e_2)( [\frac{1-\eta}{2}-1]y^e_1 + y_4^e \frac{1+\eta}{2})  \ $$

force
From prof. Vu Quoc's notes for

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega wf d \omega \ $$

$$ f^e(w) = - \int_{\Gamma_h} wh d {\Gamma_h} + \int_\omega w1 d \omega \ $$

converting to local coordinates

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) =

- \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

F(:,:,1) = [ 1/16,  1/16] [ -1/16,  1/16] [ -1/16, -1/16] [  1/16, -1/16] F(:,:,2) = [ 1/16,  1/16] [ -1/16,  1/16] [ -1/16, -1/16] [  1/16, -1/16]

Global Matrix creation
Then we constructed the Global stiffness matrix by the process of superimposing the local stiffness coefficients. Matlab output below shows the process of superimposing the elements to global stiffness matrix.

The matlab code for global matrix construction was written with respect to this element naming system, and then it was extended to encompass other numbers of elements

Original element matrices

K(:,:,1) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24]

K(:,:,2) = [ 1/24, -1/96, -1/48, -1/96] [ -1/96,  1/24, -1/96, -1/48] [ -1/48, -1/96,  1/24, -1/96] [ -1/96, -1/48, -1/96,  1/24] Superpositions to global matrix

KGlobal = 1/24         -1/96           0             -1/96          -1/48           0              0              0              0             -1/96           1/12          -1/96          -1/48          -1/48          -1/48           0              0              0              0             -1/96           1/24           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/12          -1/48           0             -1/96          -1/48           0             -1/48          -1/48          -1/48          -1/48           1/6           -1/48          -1/48          -1/48          -1/48           0             -1/48          -1/96           0             -1/48           1/12           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/24          -1/96           0              0              0              0             -1/48          -1/48          -1/48          -1/96           1/12          -1/96           0              0              0              0             -1/48          -1/96           0             -1/96           1/24

FGlobal = 1/8           0             -1/8            1/4            0             -1/4            1/8            0             -1/8

Boundary Conditions
Noting that essential BCs were specified to be on all sides of the 2x2 square. It would mean (viewing from the above global node figure) all nodes around node 5 will be tagged as essential BC.

Since we can either eliminate the nodes form Global Stiffness matrix to construct another one with Temperature ('u' function) as unknown or just '0' the node values in Global Stiffness matrix and then do inverse and multiply by Force matrix to find other non-constrained dofs.

We wrote a matlab code that automatically constrains the associate nodes and puts their contribution to force function.

Here how code applies BC's

KGlobal = 1/24         -1/96           0             -1/96          -1/48           0              0              0              0             -1/96           1/12          -1/96          -1/48          -1/48          -1/48           0              0              0              0             -1/96           1/24           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/12          -1/48           0             -1/96          -1/48           0             -1/48          -1/48          -1/48          -1/48           1/6           -1/48          -1/48          -1/48          -1/48           0             -1/48          -1/96           0             -1/48           1/12           0             -1/48          -1/96           0              0              0             -1/96          -1/48           0              1/24          -1/96           0              0              0              0             -1/48          -1/48          -1/48          -1/96           1/12          -1/96           0              0              0              0             -1/48          -1/96           0             -1/96           1/24

Kchange = 0             0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              1/6            0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0              0

FGlobal = -1/16          0              1/16          -1/8            0              1/8           -1/16           0              1/16

Fchange = -7/48         -7/48           0             -1/6            1/3            1/12          -7/48          -7/48           0

Then the matrices were reduced to the ones of the importance ( Global K matrix after the BCs where applied was isolated so no 0s would be present and attributing F function value(s) where picked by code as well)

KchangeF = 1/6

FchangeF = 1/3

Finally Temperature at unknown nodes was found by using this formula

$$ u = K^{-1}F \ $$

u = 2

Also output for 4x4 element free node values came out to be:

dcalc = 2             2              2              2

as you may see, only 1 node, in this example was present, so only 1 numerical value was generated.

Discussion
Note 1 important aspect of this output. That is, with forcing function changed to 1, the output of free degrees of freedom is still '2' as if no forcing function was present to begin with. I guess it might be explainable that because of constrained boundary conditions at all 4 ends, in order to keep the temperature constant at all 4 ends while heat was generated, the generated heat had to be instantaneously dissipated though the boundaries of the area, so that boundary temperature would stay unchanged.

As in matrix/numeral explanation, it should be noted, that because essential boundary condition is applied over all 4 sides, it has integration limits (working boundary range) as all element boundary. Therefore, with intorduction with heat source, force instead of being modified, was only scaled (remember force equation):

for f = 0 (HW 6.6)

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

for f= 1 (HW 7.1)

$$ f^e(w) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) =

- \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

you may see that force term is a scale only term. Moreover, notice the force global term

For 4 element at f = 0: FGlobal = 1/8           0             -1/8            1/4        0            -1/4            1/8            0             -1/8

Fchange = 5/24          7/48          -1/16           7/24         -1/3           -5/24           5/24           7/48          -1/16

For 4 element at f = 1:

FGlobal =

-1/16          0              1/16          -1/8            0             1/8           -1/16           0              1/16

Fchange =

-7/48         -7/48           0             -1/6           -1/3            1/12          -7/48          -7/48           0

For 9 element at f = 1

FGlobal =

-1/36          0              0              1/36          -1/18           0        0            1/18          -1/18           0        0            1/18          -1/36           0              0              1/36

Fchange =

-7/108        -7/108         -7/108          0            -11/108         5/54           5/54           1/108         -5/54         5/54           5/54           1/54          -7/108         -7/108         -7/108          0

, where boxed answer shows free nodes.

Note that force matrix DOES differ, but note that because of the similarity of force functions (scalar multiples of each other, namely, caused because forcing function is a constant and essential boundary condition acts on the same area as forcing function), the global force matrix, since it contained '0's on initial force function, will contain '0's on any other force function, as long as forcing function is scalar and essential B.C is the same as whole element B.C. Therefore, with same used stiffness matrix 'K' the adjusted force matrix with initial values of '0' on free nodes, always came out to be the same for specific number of elements with differing scalar value of forcing function.

As a result, for free node displacements, we observed the constant value of '2'.

--Eml5526.s11.team5.JA 20:51, 6 April 2011 (UTC)

Transient Term
$$ \tilde{m}(w,u) = \int_\Omega w \rho c \frac{\partial u}{\partial t} d \Omega = \int_\Omega w xy \frac{\partial u}{\partial t} d \Omega =

\int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( x_1^e \frac{1-\xi}{2} + x_2^e \frac{1+\xi}{2} \right ) \left( y_1^e \frac{1-\eta}{2} + y_4^e \frac{1+\eta}{2} \right ) \frac{\partial u}{\partial t} \left( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

Overall Equation
$$ \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

\left( x_1^e \frac{1-\xi}{2} + x_2^e \frac{1+\xi}{2} \right ) \left( y_1^e \frac{1-\eta}{2} + y_4^e \frac{1+\eta}{2} \right ) \frac{\partial u}{\partial t} \left( \frac{1}{4} A_e \ d\omega_\xi \right ) + \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1} \triangledown{b^e_1} & \triangledown{b^e_1} \triangledown{b^e_2} & \triangledown{b^e_1} \triangledown{b^e_3} & \triangledown{b^e_1} \triangledown{b^e_4}      \\ & \triangledown{b^e_2} \triangledown{b^e_2} & \triangledown{b^e_2} \triangledown{b^e_3} & \triangledown{b^e_2}\triangledown{b^e_4}      \\ & &  \triangledown{b^e_3} \triangledown{b^e_3} & \triangledown{b^e_3} \triangledown{b^e_4}      \\ symetry & &   & \triangledown{b^e_4} \triangledown{b^e_4} \end{bmatrix}

\begin{bmatrix} 1 & 0 & 0 & 0      \\ 0 & 1 &  0 & 0      \\ 0 & 0 &  1 & 0      \\ 0 & 0 &  0 & 1      \end{bmatrix}

\left ( d \frac{1}{2} (x_2^e - x_1^e)*\frac{1}{2} (y_4^e - y_1^e) d \omega_\xi \right ) = - \int_{\omega_\xi}

\begin{bmatrix} \triangledown{b^e_1}      \\ \triangledown{b^e_2}    \\ \triangledown{b^e_3}     \\ \triangledown{b^e_4} \end{bmatrix}

2* \left ( \frac{1}{4} A_e \ d\omega_\xi \right ) \ $$

Static using 2d LLEBF non-uniform mesh

Here the answer is also 2 at every node, which agrees with the uniform mesh LLEBF. There is a slight difference between nodes which can be attributed to matlab roundoff error.

Dynamic (case 2b) using 2d LLEBF uniform mesh

Dynamic (case 2b) using 2d LLEBF non-uniform mesh

Eml5526.s11.team5.savery

Eml5526.s11.team5.smith 20:37, 21 April 2011 (UTC)

Eml5526.s11.team5.smith 21:11, 21 April 2011 (UTC)

title: User:EML4500.f08.JAMAMA/FE/JA


text:
User:EML4500.f08.JAMAMA/FE/vijay

Solution

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

A) Numbering the elements and nodes of the truss structure shown in Figure 3.3.1 is an elementary procedure. The arrangement of elements and nodes in a truss structure is arbitrary, therefore their numbering may be accomplished without any calculation, as seen in Figure 3.3.2 below. $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

B) To assemble the global and stiffness force matrix, the truss in Figure 3.3.2 may be further broken down into it's separate constituent elements. Figure 3.3.3 below shows each local element and global coordinate system.  $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

Element 1 has been numbered with global nodes 1 and 4. Shown in Figure 3.3.3 it is positioned at an angle $$\phi^{(1)}=90^{o}\quad$$, with respect to the x-axis. With $$cos(90^{o})=0, \ \ sin(90^{o})=1, \ \ l^{(1)}=1m, \ \ k^{(1)}=\frac{A^{(1)}E^{(1)}}{l^{(1)}}=\frac{AE}{1m}\quad$$, the element 1 local stiffness matrix $$\textbf{K}^{(1)}\quad$$ may be solved: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

$$\textbf{K}^{(1)}=\frac{AE}{1m} \begin{matrix}
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & -1\\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}

\begin{matrix} \\ (1)\\ \\ (4)\\ \\ \end{matrix}

\\

\begin{matrix} &(1) & & (4)  & \end{matrix}

\end{matrix}$$ $$
 *  $$ \displaystyle
 * }

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

The local stiffness matrices for all elements may be derived from the following equation, courtesy of "A First Course in Finite Elements," by Jacob Fish and Ted Belytschko, 2009, pg. 30: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

$$\textbf{K}^{(e)}=k^{e}\begin{bmatrix} cos^{2}\phi^{e} &cos\phi^{e}sin\phi^{e} &-cos^{2}\phi^{e}  &-cos\phi^{e}sin\phi^{e} \\ cos\phi^{e}sin\phi^{e} &sin^{2}\phi^{e} &-cos\phi^{e}sin\phi^{e}  &-sin^{2}\phi^{e} \\ -cos^{2}\phi^{e}&-cos\phi^{e}sin\phi^{e} &cos^{2}\phi^{e}  &cos\phi^{e}sin\phi^{e} \\ -cos\phi^{e}sin\phi^{e} &-sin^{2}\phi^{e} &cos\phi^{e}sin\phi^{e}  &sin^{2}\phi^{e} \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

In a similar fashion, element 2 has been numbered with global nodes 4 and 3. Shown in Figure 3.3.3 it is positioned at an angle $$\phi^{(2)}=0^{o}\quad$$, with respect to the x-axis. With $$cos(0^{o})=1, \ \ sin(0^{o})=0, \ \ l^{(2)}=1m, \ \ k^{(2)}=\frac{A^{(2)}E^{(2)}}{l^{(2)}}=\frac{AE}{1m}\quad$$, the element 2 local stiffness matrix $$\textbf{K}^{(2)}\quad$$ may be solved: $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

$$\textbf{K}^{(2)}=\frac{AE}{1m} \begin{matrix}
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} 1 & 0 & -1 & 0\\ 0 & 0 & 0 & 0\\ -1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}

\begin{matrix} \\ (3)\\ \\ (4)\\ \\ \end{matrix}

\\

\begin{matrix} &(3) & & (4)  & \end{matrix}

\end{matrix}$$ $$
 *  $$ \displaystyle
 * }

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

Also, element 3 has been numbered with global nodes 2 and 4. Shown in Figure 3.3.3 it is positioned at an angle $$\phi^{(3)}=135^{o}\quad$$, with respect to the x-axis. With $$cos(135^{o})=\frac{-1}{\sqrt 2}, \ \ sin(135^{o})=\frac{1}{\sqrt 2}, \ \ l^{(3)}=\sqrt 2m, \ \ k^{(3)}=\frac{A^{(3)}E^{(3)}}{l^{(3)}}=\frac{AE}{\sqrt 2m}\quad$$, the element 3 local stiffness matrix $$\textbf{K}^{(3)}\quad$$ may be solved: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{K}^{(3)}=\frac{AE}{\sqrt 2m} \begin{matrix}
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} & -\frac{1}{2}\\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \end{bmatrix}

\begin{matrix} \\ (2)\\ \\ (4)\\ \\ \end{matrix}

\\

\begin{matrix} &(2) & & (4)  & \end{matrix}

\end{matrix}$$ $$
 *
 $$ \displaystyle
 * }

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

Lastly, element 4 has been numbered with global nodes 2 and 3. Shown in Figure 3.3.3 it is positioned at an angle $$\phi^{(4)}=90^{o}\quad$$, with respect to the x-axis. With $$cos(90^{o})=0, \ \ sin(90^{o})=1, \ \ l^{(4)}=1m, \ \ k^{(4)}=\frac{A^{(4)}E^{(4)}}{l^{(4)}}=\frac{AE}{1m}\quad$$, the element 4 local stiffness matrix $$\textbf{K}^{(4)}\quad$$ may be solved: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{K}^{(4)}=\frac{AE}{1m} \begin{matrix}
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 0 & -1\\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}

\begin{matrix} \\ (2)\\ \\ (3)\\ \\ \end{matrix}

\\

\begin{matrix} &(2) & & (3)  & \end{matrix}

\end{matrix}$$ $$
 *
 $$ \displaystyle
 * }

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

Direct Assembly of the local stiffness matrices yields the global stiffness matrix and the corresponding global displacement, force and reaction matrices: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

$$ K=\frac{AE}{1m} \begin{matrix}

\begin{bmatrix} 0 &0 &0 &0 &0 &0 &0 &0\\ 0 &1 &0 &0 &0 &0 &0 &-1 \\ 0 &0 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} \\ 0 &0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2}+1 & 0 & -1 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \\ 0 & 0 & 0 & 0 & 1 & 0 & -1 & 0\\ 0 & 0 & 0 & -1 & 0 & 1 & 0 & 0\\ 0 &0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2} & -1 & 0 & 1+\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} \\ 0 &-1 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} & 0 & 0 & -\frac{1}{2\sqrt{2}} & 1+\frac{1}{2\sqrt{2}} \\ \end{bmatrix}

\begin{matrix} \\ (1)\\ \\ (2)\\ \\ (3)\\ \\ (4)\\ \\ \end{matrix}

\\

\begin{matrix} (1)& & &(2)& & &(3)& & &(4) \end{matrix}

\end{matrix} $$

$$ d=\begin{bmatrix} 0\\ 0\\ 0\\ 0\\ u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix},

f=\begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 10\\ 0\\ 0\\ 0 \end{bmatrix},

r=\begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ 0\\ 0\\ 0 \end{bmatrix} $$

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

C) The global system may now be partitioned after four rows and four columns, and the nodal displacements may then be solved: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

 * style="width:95%" |
 * style="width:95%" |

$$
 *
 $$ \displaystyle
 * }

$$ K_{F}= \begin{bmatrix} 1& 0 & -1 & 0\\ 0 & 1 & 0&0\\ -1 & 0 & 1+\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ 0 & 0 & -\frac{1}{2\sqrt{2}} & 1+\frac{1}{2\sqrt{2}}\\ \end{bmatrix},

K_{EF}= \begin{bmatrix} 0& 0 & 0 & 0\\ 0& 0 & 0 & -1\\ 0 & 0 & -\frac{1}{2\sqrt2} & \frac{1}{2\sqrt2}\\ 0 & -1 & \frac{1}{2\sqrt2} & -\frac{1}{2\sqrt2} \end{bmatrix},

d_{F}= \begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix},

\bar d_{E}= \begin{bmatrix} \bar u_{1x}=0\\ \bar u_{1y}=0\\ \bar u_{2x}=0\\ \bar u_{2y}=0 \end{bmatrix},

f_{F}= \begin{bmatrix} 10\\ 0\\ 0\\ 0 \end{bmatrix},

r_{E}= \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y} \end{bmatrix}, $$

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

The unknown displacement matrix is found from the solution of the reduced system of equations below: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\frac{AE}{1m}\begin{bmatrix} 1 &0 &-1  &0 \\ 0 &1  &0  &0 \\ -1 &0  &1+\frac{1}{2\sqrt{2}}  &-\frac{1}{2\sqrt{2}} \\ 0&0 &-\frac{1}{2\sqrt{2}}  & 1+\frac{1}{2\sqrt{2}} \end{bmatrix}\begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\begin{bmatrix} 10\\ 0\\ 0\\ 0 \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$
 * style="width:95%" |
 * style="width:95%" |

\therefore \begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\begin{bmatrix} 10(1+\sqrt(2))\\ 0\\ 5(1+2\sqrt(2))\\ 5 \end{bmatrix}1*10^{-9}m$$ $$
 *
 $$ \displaystyle
 * }

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

D) Now in computing the reactions and the stresses: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$r_{E}= \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y} \end{bmatrix}=K_E\bar d_E+K_{EF}d_F=\begin{bmatrix} 0 &0 &0  &0 \\ 0 &0  &0  &-1 \\ 0 &0  &-\frac{1}{2\sqrt2}  &\frac{1}{2\sqrt2} \\ 0 &-1 &\frac{1}{2\sqrt2}  &-\frac{1}{2\sqrt2} \end{bmatrix}\cdot \frac{1m}{AE}\begin{bmatrix} 40\\ 0\\ 30\\ 10 \end{bmatrix}=\frac{1m}{AE}\begin{bmatrix} 0\\ -10\\ \frac{-10}{\sqrt2}\\ \frac{10}{\sqrt2} \end{bmatrix}=\begin{bmatrix} 0\\ -5\\ -3.53553\\ 3.535534 \end{bmatrix}1*10^{-9}N$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

Lastly, in computing each truss constituent element stress, the following equation is implemented (ref: Fish and Belytschko P30, Eq (2.47)): $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \sigma^e = E^e \frac{u_{2x}^{'e}-u_{1x}^{'e}}{l^{e}} = \frac{E^e}{l^e} \begin{bmatrix} -1 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} u_{1x}^{'e}\\ u_{1y}^{'e}\\ u_{2x}^{'e}\\ u_{2y}^{'e} \end{bmatrix} =\frac{E^e}{l^e} \begin{bmatrix} -\cos \phi^e & -\sin \phi^e & \cos \phi^e & \sin \phi^e \end{bmatrix} \mathbf d^e $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

For element 1, $$\phi^{(1)}=90^{o}\quad$$, with respect to the x-axis. ($$cos(90^{o})=0, \ \ sin(90^{o})=1$$), and $$\textbf{d}^{(1)}\quad$$ may be equated to: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{d}^{1}=\begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\frac{1m}{AE}\begin{bmatrix} 0\\ 0\\ 30\\ 10 \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(1)}=[0\ -1\  0\  1]\begin{bmatrix} 0\\ 0\\ 30\\ 10 \end{bmatrix}\frac{1}{A}=\frac{10}{A}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(1)}=500.00 \ N$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

For element 2, $$\phi^{(2)}=0^{o}\quad$$, with respect to the x-axis. ($$cos(0^{o})=1, \ \ sin(0^{o})=0$$), and $$\textbf{d}^{(2)}\quad$$ may be equated to: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{d}^{2}=\begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\frac{1m}{AE}\begin{bmatrix} 40\\ 0\\ 30\\ 10 \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(2)}=[-1\ 0\  1\  0]\begin{bmatrix} 40\\ 0\\ 30\\ 10 \end{bmatrix}\frac{1}{A}=\frac{-10}{A}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(2)}=500.00 \ N$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

For element 3, $$\phi^{(3)}=135^{o}\quad$$, with respect to the x-axis. ($$cos(135^{o})=-\frac{1}{\sqrt 2}, \ \ sin(135^{o})=\frac{1}{\sqrt 2}$$), and $$\textbf{d}^{(3)}\quad$$ may be equated to: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{d}^{3}=\begin{bmatrix} u_{2x}\\ u_{2y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\frac{1m}{AE}\begin{bmatrix} 0\\ 0\\ 30\\ 10 \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(3)}=\frac{1}{\sqrt2}[\frac{1}{\sqrt2}\ -\frac{1}{\sqrt2}\  -\frac{1}{\sqrt2}\  \frac{1}{\sqrt2}]\begin{bmatrix} 0\\ 0\\ 30\\ 10 \end{bmatrix}\frac{1}{A}=\frac{-20}{A}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(3)}=-1000.00 \ N$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

For element 4, $$\phi^{(4)}=90^{o}\quad$$, with respect to the x-axis. ($$cos(90^{o})=0, \ \ sin(90^{o})=1$$), and $$\textbf{d}^{(4)}\quad$$ may be equated to: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\textbf{d}^{4}=\begin{bmatrix} u_{2x}\\ u_{2y}\\ u_{3x}\\ u_{3y} \end{bmatrix}=\frac{1m}{AE}\begin{bmatrix} 0\\ 0\\ 40\\ 0 \end{bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(4)}=[0\ -1\  0\  1]\begin{bmatrix} 0\\ 0\\ 40\\ 0 \end{bmatrix}\frac{1}{A}=0$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\therefore \sigma^{(4)}=0.00 \ N$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

title: User:EML4500.f08.JAMAMA/FE/vijay


text:
User:EML4500.f08.JAMAMA/FE/Dahiya

Dahiya

title: User:EML4500.f08.JAMAMA/FE/Dahiya


text:
Template:TWSbase

 This template was transcluded from wikieducator - Template:TWSbase

title: Template:TWSbase


text:
User:EML4500.f08.JAMAMA/FE/srv

srv

title: User:EML4500.f08.JAMAMA/FE/srv


text:
User:EML4500.f08.JAMAMA/FE/savery

Problem

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

Consider $$ F=\left \{ 1,x,x^2,x^3,x^4 \right \} $$ on the interval $$ \Omega =[0,1] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

Find

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

1. Construct the $$ \Gamma(F)$$ matrix and observe its propeties.
 * style="width:95%" |
 * style="width:95%" |
 * }

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

 * 2. Find the determinant of $$ \Gamma(F)$$.
 * }

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

3. Conclude if F is an orthogonal family., i.e., $$ \Gamma_{ij} = \delta_{ij} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

1)

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

The Gram Matrix is defined as:
 * style="width:95%" |
 * style="width:95%" |

 * $$\Gamma_{ij}==\int_\Omega b_i(x)b_j(x)\,dx$$

$$
 *  $$ \displaystyle (eq. 1.1)
 * }

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

For the case where $$ F=\left \{ 1,x,x^2,x^3,x^4 \right \} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }
 * {| style="width:100%" border="0"

$$ b_1(x)=1$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

$$ b_2(x)=x $$ $$
 * style="width:90%" |
 * style="width:90%" |
 *  $$ \displaystyle
 * }

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

$$  b_3(x)=x^2  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

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

$$  b_4(x)=x^3  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ b_5(x)=x^4  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

Now we need to calculate each element in the $$ \Gamma $$ Matrix using eqn. 1.1: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

Since $$ \int_\Omega b_i(x)b_j(x)\,dx = \int_\Omega b_j(x)b_i(x)\,dx$$,  $$\Gamma_{ij}$$ will be equivalent to $$\Gamma_{ji}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{1,1} = \int_{0}^{1} 1*1*dx = \left [ x\right ]_{0}^{1} = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{1,2} = \Gamma_{2,1} = \int_{0}^{1} 1*x*dx = \left [ x^2/2\right ]_{0}^{1} = 1/2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{1,3} = \Gamma_{3,1} = \int_{0}^{1} 1*x^2*dx = \left [ x^3/3 \right ]_{0}^{1} = 1/3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{1,4} = \Gamma_{4,1} = \int_{0}^{1} 1*x^3*dx = \left [ x^4/4 \right ]_{0}^{1} = 1/4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{1,5} = \Gamma_{5,1} = \int_{0}^{1} 1*x^4*dx = \left [ x^5/5 \right ]_{0}^{1} = 1/5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{2,2} = \int_{0}^{1} x*x*dx = \left [ x^3/3\right ]_{0}^{1} = 1/3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{2,3} = \Gamma_{3,2} = \int_{0}^{1} x*x^2*dx = \left [ x^4/4\right ]_{0}^{1} = 1/4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{3,3} = \int_{0}^{1} x^2*x^2*dx = \left [ x^5/5 \right ]_{0}^{1} = 1/5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{3,4} = \Gamma_{4,3} = \int_{0}^{1} x^2*x^3*dx = \left [ x^6/6\right ]_{0}^{1} = 1/6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{3,5} = \Gamma_{5,3} = \int_{0}^{1} x^2*x^4*dx = \left [ x^7/7\right ]_{0}^{1} = 1/7 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{4,4} = \int_{0}^{1} x^3*x^3*dx = \left [ x^7/7\right ]_{0}^{1} = 1/7 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{4,5} = \Gamma_{5,4} = \int_{0}^{1} x^3*x^4*dx = \left [ x^8/8\right ]_{0}^{1} = 1/8 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \Gamma_{5,5} = \Gamma_{5,5} = \int_{0}^{1} x^4*x^4*dx = \left [ x^9/9\right ]_{0}^{1} = 1/9 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

Assembling the matrix we get: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

\Gamma_ \left (F \right) = \begin{bmatrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5}\\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8}\\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{bmatrix} $$

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

The matrix is observed to only be symmetric $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

2)

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

Using the matlab function det(X) the determinant of $$\Gamma$$ can be found to be: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

 * $$ \det \left ( \Gamma \right ) = 0 $$

3)

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

To conclude if $$\Gamma$$ is orthogonal we compare it to a matrix formed by the Kronecker Delta. $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

The Kronecker Delta is defined by: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$\delta_{ij} = \begin{cases} 1, & \mbox{ for i=j} \\ 0, & \mbox{ for i≠j} \end{cases} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

Thus the general form of a the matrix is: $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

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

$$ \begin{bmatrix} 1    & \cdots & 0      \\ \vdots & \ddots & \vdots \\ 0     & \cdots & 1 \end{bmatrix}
 * style="width:95%" |
 * style="width:95%" |

$$ $$
 *
 $$ \displaystyle
 * }

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

Thus the $$ \Gamma $$ Matrix is not orthogonal $$
 * style="width:95%" |
 * style="width:95%" |
 *
 $$ \displaystyle
 * }

--Eml5526.s11.team5.savery 19:49, 2 February 2011 (UTC) = Problem 4.2, Exercise = 3.5, Fish and Belytschko =

Problem
Obtain the Weak Form for the equation of Heat Conduction with the Boundary Condition $$ T(0)= \ 100$$ and $$ \overline{q}= \ hT $$. The condition on the right is a convection condition.

Solution
{| style="width:100%" border="0" Here $$ \quad \tilde{m}(x) =\quad m(x).c$$. This is all due to change in temperature i.e. due to Convection. $$ \mathbf{Boundary \quad Conditions:}$$
 * style="width:95%" |
 * style="width:95%" |

on above Eq. 4.2.7,We get:

{| style="width:100%" border="0" As we know that $$ w \quad $$ has to satisfy Homogeneous Essential Boundary Condition i.e. $$ w(0)= 0 \quad $$. So, Eq.4.2.8 becomes
 * style="width:95%" |
 * style="width:95%" |

title: User:EML4500.f08.JAMAMA/FE/savery