User:Eml5526.s11.team3.akj/Homework 7

=Problem 7.1 Static and Dynamic Finite Element Modelling and Analysis for Heat Problem=

Given: Problem defined on [[media:fe1.s11.mtg40.djvu|Mtg 39 (a)]]
Strong Form

for the static case

Initial condition for dynamic case

Boundary conditions for Dynamic case

Find
1. Find the static solution $$ \quad U_{s}^h $$

2. Dynamic (Transient) Solution

2b. Plot $$ u^h(0,t) \ vs \ t $$ at center till convergence to steady state

ii. Finding approximate solution
2DLIBF are given by

$$L_{i,m}$$ is the Lagrange basis function along x and $$L_{i,m}$$ is the Lagrange basis function along y given by

The LIBF be defined by a function

In order to solve for the free degree of freedom $${{d}_{F}},_ – ^ – $$, we construct the Conductance matrix, Heat Source  matrix and Capacitance matrix as follow:
 * {| style="width:100%" border="0"

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

\tilde{\kappa}(\omega,u)=\int_{\Omega}\bigtriangledown \omega.\boldsymbol{\kappa }.\bigtriangledown ud\Omega

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

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

\tilde{f}(\omega)=-\int_{\Gamma _{h}}\omega h d\Gamma _{h}+\int_{\Omega}\omega f d\Omega $$
 * }


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

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

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

To solve for the static case from Eq. 1.2 since $$\frac{\partial u}{\partial t}\left( {x,t} \right) = 0, f\left( {x,t} \right) = 1 $$ hence
 * {| style="width:100%" border="0"

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

\tilde{m}(\omega,u)=0 $$
 * }


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

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

\tilde{f}(\omega)= \int_{\Omega}\omega d\Omega $$
 * }


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

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

\tilde{\kappa}(\omega,u)=\int_{\Omega}\bigtriangledown \omega.\bigtriangledown ud\Omega

$$
 * }

$$\displaystyle

\tilde{m}(\omega),\tilde{f}(\omega,u),\tilde{\kappa}(\omega,u)$$ calculate all the degrees of freedom of the system.

Approximated solution $$ u^{h} \left ( x \right ) $$ and $$ w^{h} \left ( x \right ) $$ :

where $$ \left ( c_{i} \right ) $$ and $$ \left ( d_{j} \right ) $$ are constants and $$ \left ( b_{i} \right ) $$ and $$ \left (b_{j}\right ) $$ are the LIBF. We will use these LIBF as a shape function satisfying the CBS.

where,

Capacitance matrix

Conductivity matrix

Force matrix

But we are interested in finding out the degree of freedom at the free ends.

So have to solve for the free degrees of freedom. For this, the MATLAB code is generated such that free ends are collected at one end.

The mas matrix $$\mathbf{\tilde{M}}$$ is given by

The $$\displaystyle \tilde{\kappa} $$ matrix is given by

The force $$\displaystyle  \tilde{F} $$ matrix is given by

Since we only solve for the the free degree of freedom, we only use the equation:-


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

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |

\mathbf{M_{FF}}\mathbf\dot+ \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{M_{FF}}\mathbf\dot{g}+\mathbf{K_{FE}}\mathbf{g})

$$ (Eq 1.17 )
 * 
 * }

In the existing problem, $$\displaystyle \mathbf\dot=0 $$ and $$\displaystyle \mathbf{g}$$ = constant. Thus $$\displaystyle \mathbf\dot{g} = 0.$$

Taking this considerations, Eq 1.17 reduces to :-
 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |

\mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-\mathbf{K_{FE}}\mathbf{g}

$$ (Eq 1.18 )
 * 
 * }

We can reaarange to find $$\displaystyle \mathbf{d_{F}} $$ as:-

Taking this considerations, Eq 1.17 reduces to :-
 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |

\mathbf{d_{F}}=\mathbf{K_{FF}^{-1}}(\mathbf{F_{F}}-\mathbf{K_{FE}}\mathbf{g})

$$ (Eq 1.19 )
 * 
 * }

The Matlab code to generate the plots and LIBF is:



Solution for dynamic case
In the dynamic case, $$\displaystyle \mathbf\dot \ne 0 $$ and $$\displaystyle \mathbf\dot{g} = 0.$$

So, Eq 1.17 can be re-written as:-


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

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |

\mathbf{M_{FF}}\mathbf\dot+ \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g})$$

$$\displaystyle \mathbf{K_{FF}}\mathbf{d_{F}}=\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g}+\mathbf{M_{FF}}\mathbf\dot)$$

$$\displaystyle \mathbf{d_{F}}=\mathbf{K_{FF}^{-1}}(\mathbf{F_{F}}-(\mathbf{K_{FE}}\mathbf{g}+\mathbf{M_{FF}}\mathbf\dot)) $$ (Eq 1.20 )
 * 
 * }

the initial conditions obtained is:


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

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

\mathbf{d(o)}=\mathbf{M_{FF}^{-1}}\mathbf{G}

$$ where
 * }


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

$$\displaystyle
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |

G={{\langle {{b}_{i}},{{u}_{0}}\rangle }_{M}}=\int\limits_{\Omega }\rho c{{u}_{o}}d\Omega $$

$$\displaystyle \mathbf{u}=\mathbf{xy}

$$ (Eq 1.21)
 * 
 * }

Using MATLAB code "ode45" and the Equation 1.12 using the boundary conditions in 1.13, we solve the problem.

The Matlab code to generate the plots and LIBF is: