# 7.1 2D heat problem and here

## Given: 2D Strong Form and BCs

 ${\displaystyle \Omega ={\bar {w}}=\Box }$ biunit square ${\displaystyle \displaystyle (Eq.7.1.1)}$

 ${\displaystyle div({\underline {K}}\cdot grad(u))+f=\rho c{\frac {\partial u}{\partial t}}}$ ${\displaystyle \displaystyle (Eq.7.1.2)}$

 ${\displaystyle {\underline {K}}={\underline {I}}}$ (identity matrix) ${\displaystyle \displaystyle (Eq.7.1.3)}$

 ${\displaystyle f=0\quad }$ , ${\displaystyle {\frac {\partial u}{\partial t}}=0}$ ${\displaystyle \displaystyle (Eq.7.1.4)}$

 Essential boundary condition: ${\displaystyle g=2\quad }$ on ${\displaystyle \partial \Omega \quad }$ ${\displaystyle \displaystyle (Eq.7.1.5)}$

 Natural boundary condition: none $\displaystyle \displaystyle$

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

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

 1) Static (steady state): ${\displaystyle f(x)=1\quad }$ in ${\displaystyle \Omega =\Box }$ 1a) 2d LIBF 1b) 2d LLEBF 1b1) Uniform mesh 1b2) Non-uniform mesh 2) Dynamic(transient):${\displaystyle \rho c=3\quad }$ Initial condition ${\displaystyle u(x,t=0)=xy\qquad \forall x\in \Omega }$ $\displaystyle \displaystyle$

 2a) ${\displaystyle f(x,t)=0\qquad \forall x\in \Omega ,\ \forall t>0}$ ${\displaystyle g(x,t)=2{\mbox{ on }}\Gamma _{g}=\partial \Omega }$ 2d LIBF $\displaystyle \displaystyle$

 2b) ${\displaystyle f(x,t)=1\qquad \forall x\in \Omega ,\ \forall t>0}$ ${\displaystyle g(x,t)=2{\mbox{ on }}\Gamma _{g}=\partial \Omega }$ 2b1) 2d LIBF 2b2) 2d LLEBF 2b2a) Uniform mesh 2b2b) non-uniform mesh $\displaystyle \displaystyle$

# Solution

In general by Prof. Vu Quoc's notes

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

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

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

${\displaystyle {\tilde {f}}(w)=-\int _{\Gamma _{h}}whd\Gamma _{h}+\int _{\Omega }wfd\Omega \ }$

${\displaystyle K_{ij}^{e}=\int _{\Omega ^{e}}\triangledown {b_{i}^{e}}'(x,y)\kappa (x,y)\triangledown {b_{j}^{e}}'(x,y)d\omega _{x}\ }$ [1]

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

For linear 2D case with boundary being square

Bases, ${\displaystyle b_{i}\ }$ are defined to be ${\displaystyle L_{[I,J]}(x,y)=N_{[I,J]}^{e}(x,y)\ }$

where ${\displaystyle N_{[I,J]}^{e}(x,y)=N_{I}^{e}(x)N_{J}^{e}(y)\ }$ and ${\displaystyle 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 ${\displaystyle y_{4}^{e}}$ for ${\displaystyle y_{3}^{e}\ }$ and ${\displaystyle y_{1}^{e}}$ for ${\displaystyle y_{2}^{e}\ }$ and the same for other coordinates

${\displaystyle N_{[1,1]}^{e}(x,y)={\frac {x^{e}-x_{2}^{e}}{x_{1}^{e}-x_{2}^{e}}}{\frac {y^{e}-y_{4}^{e}}{y_{1}^{e}-y_{4}^{e}}}={\frac {1}{4}}(x-\xi )(y-\eta )=b_{1}^{e}(x,y)\ }$

${\displaystyle N_{[2,1]}^{e}(x,y)={\frac {x^{e}-x_{1}^{e}}{x_{2}^{e}-x_{1}^{e}}}{\frac {y^{e}-y_{4}^{e}}{y_{1}^{e}-y_{4}^{e}}}={\frac {1}{4}}(x-\xi )(y-\eta )=b_{2}^{e}(x,y)\ }$

${\displaystyle N_{[2,2]}^{e}(x,y)={\frac {x^{e}-x_{1}^{e}}{x_{2}^{e}-x_{1}^{e}}}{\frac {y^{e}-y_{1}^{e}}{y_{4}^{e}-y_{1}^{e}}}={\frac {1}{4}}(x-\xi )(y-\eta )=b_{3}^{e}(x,y)\ }$

${\displaystyle N_{[1,2]}^{e}(x,y)={\frac {x^{e}-x_{2}^{e}}{x_{1}^{e}-x_{2}^{e}}}{\frac {y^{e}-y_{1}^{e}}{y_{4}^{e}-y_{1}^{e}}}={\frac {1}{4}}(x-\xi )(y-\eta )=b_{4}^{e}(x,y)\ }$

${\displaystyle \mathbf {BN} ={\begin{bmatrix}{\frac {\partial N_{[1,1]}^{e}(x,y)}{\xi }}&{\frac {\partial N_{[2,1]}^{e}(x,y)}{\xi }}&{\frac {\partial N_{[2,2]}^{e}(x,y)}{\xi }}&{\frac {\partial N_{[1,2]}^{e}(x,y)}{\xi }}\\{\frac {\partial N_{[1,1]}^{e}(x,y)}{\eta }}&{\frac {\partial N_{[2,1]}^{e}(x,y)}{\eta }}&{\frac {\partial N_{[2,2]}^{e}(x,y)}{\eta }}&{\frac {\partial N_{[1,2]}^{e}(x,y)}{\eta }}\end{bmatrix}}\ }$

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

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

[/itex]

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

${\displaystyle \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 Prof. Vu Quoc's Notes

${\displaystyle 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

${\displaystyle 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 ${\displaystyle K^{e}\ }$. Also note that since we have, 4 degrees of freedom in 1 element node, we do expect the gram matrix ${\displaystyle \Gamma \ }$ to be 4x4, therefore, ${\displaystyle \kappa _{4x4}\ }$

${\displaystyle K_{ij}^{e}=\int _{\Omega ^{e}}\triangledown {b_{i}^{e}}'(x,y)\kappa (x,y)\triangledown {b_{j}^{e}}'(x,y)d\omega _{x}\ }$

${\displaystyle K_{ij}^{e}=\int _{-1}^{+1}{\begin{bmatrix}\triangledown {b_{1}^{e}}\triangledown {b_{1}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{2}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{4}^{e}}\\&\triangledown {b_{2}^{e}}\triangledown {b_{2}^{e}}&\triangledown {b_{2}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{2}^{e}}\triangledown {b_{4}^{e}}\\&&\triangledown {b_{3}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{3}^{e}}\triangledown {b_{4}^{e}}\\symetry&&&\triangledown {b_{4}^{e}}\triangledown {b_{4}^{e}}\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 ${\displaystyle 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 ${\displaystyle \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

${\displaystyle b_{1}^{e}(\xi ,\eta )={\frac {1}{A^{e}}}(x_{1}^{e}{\frac {1-\xi }{2}}+[{\frac {1+\xi }{2}}-1]x_{2}^{e})(y_{1}^{e}{\frac {1-\eta }{2}}+[{\frac {1+\eta }{2}}-1])y_{4}^{e})\ }$

${\displaystyle b_{2}^{e}(\xi ,\eta )=-{\frac {1}{A^{e}}}([{\frac {1-\xi }{2}}-1]x_{1}^{e}+x_{2}^{e}{\frac {1+\xi }{2}})(y_{1}^{e}{\frac {1-\eta }{2}}+[{\frac {1+\eta }{2}}-1]y_{4}^{e})\ }$

${\displaystyle b_{3}^{e}(\xi ,\eta )={\frac {1}{A^{e}}}([{\frac {1-\xi }{2}}-1]x_{1}^{e}+x_{2}^{e}{\frac {1+\xi }{2}})([{\frac {1-\eta }{2}}-1]y_{1}^{e}+y_{4}^{e}{\frac {1+\eta }{2}})\ }$

${\displaystyle b_{4}^{e}(\xi ,\eta )=-{\frac {1}{A^{e}}}(x_{1}^{e}{\frac {1-\xi }{2}}+[{\frac {1+\xi }{2}}-1]x_{2}^{e})([{\frac {1-\eta }{2}}-1]y_{1}^{e}+y_{4}^{e}{\frac {1+\eta }{2}})\ }$

### force

From prof. Vu Quoc's notes for elemental force expression

${\displaystyle f^{e}(w)=-\int _{\Gamma _{h}}whd{\Gamma _{h}}+\int _{\omega }wfd\omega \ }$

${\displaystyle f^{e}(w)=-\int _{\Gamma _{h}}whd{\Gamma _{h}}+\int _{\omega }w1d\omega \ }$

converting to local coordinates

${\displaystyle f^{e}(w)=-\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}2*\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)+\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)=-\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\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

${\displaystyle 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)

${\displaystyle f^{e}(w)=-\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}2*\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)\ }$

for f= 1 (HW 7.1)

${\displaystyle f^{e}(w)=-\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}2*\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)+\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)=-\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\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

${\displaystyle {\tilde {m}}(w,u)=\int _{\Omega }w\rho c{\frac {\partial u}{\partial t}}d\Omega =\int _{\Omega }wxy{\frac {\partial u}{\partial t}}d\Omega =\int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\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

${\displaystyle \int _{\omega _{\xi }}{\begin{bmatrix}\triangledown {b_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\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_{1}^{e}}\triangledown {b_{1}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{2}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{1}^{e}}\triangledown {b_{4}^{e}}\\&\triangledown {b_{2}^{e}}\triangledown {b_{2}^{e}}&\triangledown {b_{2}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{2}^{e}}\triangledown {b_{4}^{e}}\\&&\triangledown {b_{3}^{e}}\triangledown {b_{3}^{e}}&\triangledown {b_{3}^{e}}\triangledown {b_{4}^{e}}\\symetry&&&\triangledown {b_{4}^{e}}\triangledown {b_{4}^{e}}\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_{1}^{e}}\\\triangledown {b_{2}^{e}}\\\triangledown {b_{3}^{e}}\\\triangledown {b_{4}^{e}}\end{bmatrix}}2*\left({\frac {1}{4}}A_{e}\ d\omega _{\xi }\right)\ }$

### References

NOTE THAT FOR LLIB CASE, THERE WAS ONLY 1 ELEMENT PRESENT (THE GLOBAL MATRIX ELEMENT) AND SINCE ALL BCs ON ELEMENT WERE CONSTRAINED BY ESSENTIAL BOUNDARIES, DUE TO THE LACK OF A FREE ELEMENT, A SOLUTION COULD NOT BE FOUND REGARDLESS OF THE NUMBER OF BASIS FUCNTIONS.

## Static using 2d LLEBF uniform mesh

N=4
N=6
N=8
N=10

As you can see for this case the temperature is 2 at every node which is to be expected.

## Static using 2d LLEBF non-uniform mesh

N=4
N=6
N=8
N=10

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.

N=2
N=3
N=4
N=5

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

N=3
N=4
N=5

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

## Given

### P.D.E

 ${\displaystyle Tdiv(gradu)+f=\rho \ {\frac {\partial ^{2}u}{\partial t^{2}}}}$ ${\displaystyle \displaystyle (Eq.7.2.1)}$

 ${\displaystyle \quad T=4}$ ${\displaystyle \displaystyle (Eq.7.2.2)}$

 ${\displaystyle f(\mathbf {x} )=1}$ ${\displaystyle \displaystyle (Eq.7.2.3)}$

 ${\displaystyle \quad \rho =3}$ ${\displaystyle \displaystyle (Eq.7.2.4)}$

 ${\displaystyle \quad Basis\ Function\ :\ 2D\ Lagrange\ Interpolation\ Basis\ Function}$ $\displaystyle \displaystyle$

### Essential Boundary Condition

 ${\displaystyle \quad g=0}$ ${\displaystyle \displaystyle (Eq.7.2.5)}$

## Find

 1. Find the static solution ${\displaystyle \quad U^{h}}$ to ${\displaystyle \quad 10^{-6}}$ accuracy at the center. $\displaystyle \displaystyle$

### Dynamic(transient)

 1. Fundamental Frequency and Fundamental Period $\displaystyle \displaystyle$

 2. Plot ${\displaystyle u^{h}(0,t)\ vs\ t\quad for\ t\ \epsilon \ [0,2T_{1}]}$ for symmetric and Non-Symmetric ${\displaystyle \quad u^{0}}$ $\displaystyle \displaystyle$

 3. Produce a movie of the vibrating membrane for symmetric and Non-Symmetric ${\displaystyle \quad u^{0}}$ $\displaystyle \displaystyle$

## Solution

Substituting the given values in the P.D.E, we get,

 ${\displaystyle \quad 4\ div(gradu)\ +\ 1\ =\ 0}$ $\displaystyle \displaystyle$

 ${\displaystyle 4\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)+1=0}$ $\displaystyle \displaystyle$

 ${\displaystyle {\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}={\frac {-1}{4}}}$ ${\displaystyle \displaystyle (Eq.7.2.6)}$

2D LIBF

 ${\displaystyle N_{I}(x,y)=L_{i,m}(x).L_{j,n}(y)\quad Where\ I=i+(j-1)m}$ ${\displaystyle \displaystyle (Eq.7.2.7)}$

Where,

 ${\displaystyle L_{i,m}(x)=\prod _{k=1\neq i}^{m}{\frac {x-x_{k}}{x_{i}-x_{k}}}}$ $\displaystyle \displaystyle$
 ${\displaystyle L_{j,n}(x)=\prod _{k=1\neq j}^{m}{\frac {y-y_{k}}{y_{j}-y_{k}}}}$ $\displaystyle \displaystyle$

Consider the Weak Form

 ${\displaystyle {\tilde {m}}(w,u)+{\tilde {k}}(w,u)={\tilde {f}}(w)}$ ${\displaystyle \displaystyle (Eq.7.2.8)}$

where,

 ${\displaystyle m(w,u)=\int \limits _{\Omega }{w\rho c{\frac {\partial u}{\partial t}}}d\Omega }$ $\displaystyle \displaystyle$

 ${\displaystyle k(w,u)=\int \limits _{\Omega }{\nabla wK\nabla u}d\Omega }$ $\displaystyle \displaystyle$

 ${\displaystyle f(w,u)=-\int _{\Gamma h}{whd{\Gamma _{h}}}+\int \limits _{\Omega }{wf}d\Omega }$ $\displaystyle \displaystyle$

### Static case

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

 ${\displaystyle \quad Kd=F}$ $\displaystyle \displaystyle$

where,

 ${\displaystyle {K_{ij}}=\int _{w^{e}}{(\nabla N_{I}^{e})I(\nabla N_{J}^{e})d{w^{e}}}}$ ${\displaystyle \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 ${\displaystyle \quad \rho \ =\ 3}$

#### Free Vibration: Solve Eigen Value problem

The equation given to solve the Eigen value problem is,

 ${\displaystyle \mathbf {K} \mathbf {\phi } =\lambda \mathbf {M} \mathbf {\phi } }$ ${\displaystyle \displaystyle (Eq.7.2.10)}$

Where,

 ${\displaystyle \quad \mathbf {\phi } -Eigen\ vectors}$ $\displaystyle \displaystyle$

 ${\displaystyle \quad \lambda -Eigen\ values}$ $\displaystyle \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

 ${\displaystyle u(\mathbf {x} ,t=0)=u_{s}^{h}(x)\ \forall \ \mathbf {x} \ \epsilon \ \Omega }$ ${\displaystyle \displaystyle (Eq.7.2.11)}$

 ${\displaystyle {\dot {u}}(\mathbf {x} ,t=0)=0\ \forall \ \mathbf {x} \ \epsilon \ \Omega }$ ${\displaystyle \displaystyle (Eq.7.2.12)}$

 ${\displaystyle T=4,\qquad f(\mathbf {x} ,t)=0\ in\ \Omega }$ ${\displaystyle \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

For T = 0,

#### Transient Analysis with Non-symmetric Uo

In this case, initial u value is given by,

 ${\displaystyle 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 }$ ${\displaystyle \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

--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 Static Solution of a Circle Using Quadrilateral and Triangular Elements

## Given/Find: Static Solution of the Circle

 Let Omega a circle with a one-unit radius. Find the static solution such that: ${\displaystyle T=4}$ ${\displaystyle f({\underline {x}})=1}$ in Omega ${\displaystyle g=0}$ on ${\displaystyle \Gamma _{g}=\partial \Omega }$ 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 \displaystyle$

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

 The Partial Differential Equation for a vibrating membrane is the following via transient analysis: $\displaystyle \displaystyle$

 ${\displaystyle T({\textrm {div\ grad\ u}})+{\textrm {f}}=\rho {\frac {\partial ^{2}u}{\partial t^{2}}}}$ $\displaystyle \displaystyle$

 Where: ${\displaystyle u({\underline {x}},t)=}$ Transverse Displacement ${\displaystyle T=\quad }$ Tension (Force/Length) = Constant ${\displaystyle {\textrm {f}}({\underline {x}},t)=}$ Distributed Transverse Force ${\displaystyle \rho ({\underline {x}})=}$ Mass Density (Mass/Area) & ${\displaystyle {\textrm {div\ grad\ u}}={\frac {\partial }{\partial x_{i}}}{\frac {\partial u}{\partial x_{j}}}}$ $\displaystyle \displaystyle$

 The Weighted Residual Form is as follows: $\displaystyle \displaystyle$

 ${\displaystyle \int _{\Omega }w\left[T{\frac {\partial ^{2}u}{\partial x_{i}\partial x_{j}}}+{\textrm {f}}-\rho {\frac {\partial ^{2}u}{\partial t^{2}}}\right]d\Omega =0}$ $\displaystyle \displaystyle$

 Which can thereby be broken into three separate terms: $\displaystyle \displaystyle$

 ${\displaystyle \int _{\Omega }wT{\frac {\partial ^{2}u}{\partial x_{i}\partial x_{j}}}d\Omega =\mathbb {K} }$ $\displaystyle \displaystyle$

 ${\displaystyle \int _{\Omega }w{\textrm {f}}\ d\Omega }$ $\displaystyle \displaystyle$

 ${\displaystyle \int _{\Omega }w\rho {\frac {\partial ^{2}u}{\partial t^{2}}}\ d\Omega ={\tilde {m}}(w,u)}$ $\displaystyle \displaystyle$

 This final equation will be used later in the Weak Form. For now, the first term above (set equal to ${\displaystyle \mathbb {K} }$), is equivalent to: $\displaystyle \displaystyle$

 ${\displaystyle \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 }$ $\displaystyle \displaystyle$

 And by using divergence theorem, the first term of the equation above may be converted to yield: $\displaystyle \displaystyle$

 ${\displaystyle \mathbb {K} =\int _{\partial \Omega }n_{i}wT{\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 }$ $\displaystyle \displaystyle$

 Now since $\displaystyle \displaystyle$

 ${\displaystyle \partial \Omega =\Gamma _{g}\cup \Gamma _{h}\ {\textrm {but}}\ w=0\ \ {\textrm {on}}\ \ \Gamma _{g}}$ $\displaystyle \displaystyle$

 and the following are true: $\displaystyle \displaystyle$

 ${\displaystyle {\frac {\partial w}{\partial x_{i}}}={\textrm {grad}}\ w=\bigtriangledown w}$ $\displaystyle \displaystyle$

 ${\displaystyle {\frac {\partial u}{\partial x_{j}}}={\textrm {grad}}\ u=\bigtriangledown u}$ $\displaystyle \displaystyle$

 Then the K equation becomes: $\displaystyle \displaystyle$

 ${\displaystyle \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 }$ $\displaystyle \displaystyle$

 Otherwise written as: $\displaystyle \displaystyle$

 ${\displaystyle \mathbb {K} =\int _{\Gamma _{h}}w\ {\underline {n}}T\cdot {\textrm {grad\ u}}\ \partial \Gamma _{h}-T\int _{\Omega }\bigtriangledown w\bigtriangledown u\ d\Omega }$ $\displaystyle \displaystyle$

 where the last term is conveniently used to write ultimately the weak form of the weighted residual form which is the following: $\displaystyle \displaystyle$

 ${\displaystyle {\tilde {\textrm {f}}}(w)={\tilde {m}}(w,u)+{\tilde {k}}(w,u)\ \ \ \forall w\ \ \ s.t.\ \ \ w=0\ \ \ on\ \ \ \Gamma _{g}}$ $\displaystyle \displaystyle$

 ${\displaystyle {\tilde {m}}(w,u)=\int _{\Omega }w\rho {\frac {\partial ^{2}u}{\partial t^{2}}}\ d\Omega }$ $\displaystyle \displaystyle$

 ${\displaystyle {\tilde {k}}(w,u)=T\int _{\Omega }\bigtriangledown w\bigtriangledown u\ d\Omega }$ $\displaystyle \displaystyle$

 ${\displaystyle {\tilde {\textrm {f}}}(w)=\int _{\Gamma _{h}}w{\underline {n}}\ T\ {\textrm {gradu}}\ \partial \Gamma _{h}+\int _{\Omega }w{\textrm {f}}\ d\Omega }$ $\displaystyle \displaystyle$

 From the Discrete Weighted Form, $\displaystyle \displaystyle$

 ${\displaystyle {\underline {M}}_{FF}{\underline {d}}_{F}+{\underline {K}}_{FF}{\underline {d}}_{F}={\underline {F}}_{F}-({\underline {M}}_{FE}g+{\underline {K}}_{FE}g)}$ $\displaystyle \displaystyle$

 But because g=0 from the boundary conditions, as well as M due to the static requirement, the displacements may be derived: $\displaystyle \displaystyle$

 ${\displaystyle {\underline {K}}{\underline {d}}={\underline {F}}}$ $\displaystyle \displaystyle$

 ${\displaystyle \therefore {\underline {d}}={\underline {K}}^{-1}{\underline {F}}}$ $\displaystyle \displaystyle$

 Now to solve for the displacements: $\displaystyle \displaystyle$

 ${\displaystyle N_{I}(x,y)=L_{i,m}(x)\cdot L_{j,n}(y)}$ $\displaystyle \displaystyle$

 Where $\displaystyle \displaystyle$

 ${\displaystyle I=i+(j-1)m}$ and $\displaystyle \displaystyle$

 ${\displaystyle L_{i,m}(x)=\prod _{K=1,K\neq i}^{m}\left({\frac {x-x_{k}}{x_{i}-x_{k}}}\right)}$ $\displaystyle \displaystyle$

 And setting the Shape Function to the following: $\displaystyle \displaystyle$

 ${\displaystyle 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}}}$ $\displaystyle \displaystyle$

 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: $\displaystyle \displaystyle$

 ${\displaystyle X_{1}^{e}=\varphi _{1}^{e}(\xi )=\sum _{I=1}^{nel}N_{I}(\xi )x_{1,I}^{e}}$ $\displaystyle \displaystyle$

 ${\displaystyle (y=)X_{2}^{e}=\varphi _{2}^{e}(\xi )=\sum _{I=1}^{nel}N_{I}(\xi )x_{2,I}^{e}}$ $\displaystyle \displaystyle$

 Using the following Jacobian matrix: $\displaystyle \displaystyle$

 ${\displaystyle 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}}}$ $\displaystyle \displaystyle$

 The following equations may now be utilized: $\displaystyle \displaystyle$

 ${\displaystyle {\underline {B}}_{I}(\xi )=\bigtriangledown _{x}N_{I}^{e}(\xi )=({\underline {J}}^{e})^{-T}(\xi )\bigtriangledown _{\xi }N_{I}^{e}(\xi )}$ $\displaystyle \displaystyle$

 ${\displaystyle k_{I3}^{e}=\int _{\omega =el}N_{I}(\xi )\cdot T\cdot {\underline {B}}_{J}(\xi )J(\xi )d(el)}$ $\displaystyle \displaystyle$

 ${\displaystyle F_{i}=\int _{\omega =el}N_{I}(\xi )f\ det(J(\xi ))d(el)}$ $\displaystyle \displaystyle$

 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: $\displaystyle \displaystyle$

 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. $\displaystyle \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