User:Eml5526.s11.team3/Homework 7

 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.20 using the boundary conditions in 1.21, we solve the problem.



The Matlab code to generate the plots and LIBF is:

=Problem 7.2 Static and Dynamic Finite Element Modelling and Analysis for Vibrating Membrane=

Given: Problem defined on [[media:fe1.s11.mtg40.djvu|Mtg 40 (c)]]
Refer to Lecture notes [[media:fe1.s11.mtg40.djvu|Lecture 40-5]] [[media:fe1.s11.mtg41.djvu|Lecture 41-1,41-2]] for more detailed description.

Strong Form

Essential boundary conditions,

Find
1. Find the static solution $$ \quad U_{s}^h $$ to $$ \quad 10^{-6} $$ accuracy at x = 0.

2. Dynamic (Transient) Solution

2a. Free Vibration Analysis

2b. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1]$$  and produce a movie of the vibrating membrane for symmetric

2c. Plot $$ u^h(0,t) \ vs \ t \quad  for \ t \ \epsilon \ [0,2T_1]$$  and produce a movie of the vibrating membrane for non-symmetric

i. Weak Form
From the PDE the weak form can be written as:

Now from the theory of Integration by parts, we can write

Hence Eq. 2.1 can be now written as

Now we can write $$\displaystyle d\Omega $$ as $$\displaystyle d\Omega ={{\Gamma }_{g}}\cup {{\Gamma }_{h}}$$, where $$\displaystyle {{\Gamma }_{g}}$$ = essential boundary condition. $$\displaystyle {{\Gamma }_{h}}$$ = natural boundary condition.

So by applying the Gauss theorem on the first term we obtain,

We select $$\displaystyle w $$ such that $$\displaystyle w=0 $$ on $$\displaystyle {{\Gamma }_{g}}$$

This Equation is of the form:-

where,

ii. Finding approximate solution
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.

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

where,

Capacitance matrix

Conductivity matrix

Force matrix

1. Static solution
We have strong form of vibrating membrane equation (Eq. 2.1) where T = 4, $$f\left( x \right){\text{  =  1}}$$,  $$\quad $$   $$\quad \rho = 3 $$,  $$\quad $$ $$\qquad\frac{\partial^2u}{\partial t^2}=0$$

we can also write for 2D

$$\frac{{{\partial }^{2}}u}{\partial x_{i}^{2}}= \frac{\partial^2u}{\partial x_{1}^2} + \frac{\partial^2u}{\partial x_{2}^2}$$

If we rearrange Eq 2.1 for static case conditions, we get



2. Dynamic solution
$$ \quad \rho \ = \ 3$$

2a. Free Vibration Analysis
The equation is,

Where, $$ \quad \lambda - Eigenvalues $$

$$\quad \mathbf{\phi} - Eigenvectors $$

We calculate eigenvalues and then put them into the equation below to find frequencies.

The fundamental frequency is $$\omega _{1}=\lambda_{1}^{1/2}$$, and  $$f_{1}=\frac{\omega _{1}}{2\Pi }$$.

The period T is $$T_{1}=\frac{1}{f_{1}}$$

2b. Plotting displacement and the movie of the vibrating membrane for symmetric
Initial conditions:

Note: To see the movie click on the picture twice and wait for a second for movie to start.



2c. Plotting displacement and the movie of the vibrating membrane for non-symmetric
Initial conditions:

Note: To see the movie click on the picture twice and wait for a second for movie to start.



=Problem 7.3 Find static solution to $$10^{-6}\;\;$$ accuracy using Quadrilateral and Triangular elements=

Problem Statement
Let $$\;\;\Omega\;\;$$ be a circular domain. Find static solution such that T = 4 and $$\;\;f(x) = 1\;\;\;\;$$in $$\;\;\;\;\Omega\;\; and \;\;g = 0 \;\;on \;\;\Gamma_{g}= \partial\Omega\;\;$$. Use both quadrilateral and triangular elements. Plot the deformed shape in 3-D.

Solution
The nodes and the elements information were obtained from Calculix and are shown below, for both the cases, quadrilateral elements and triangular elements.



We ran the program for the case of quadrilateral elements. We were unable to get the code to run till the end as it took a lot of time. Hence we had to terminate the program.

=Problem 7.4=

Problem statement
Please refer to lecture note 37-2 and 38 for more details.

Solve the heat conduction problem in 2D with boundary convection.

1> Construct the weak form for heat conduction in 2D with boundary convection.

2> Develop semi-discrete equations (ODEs) similar to (3) 23-3 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.



Developing the Weak Form of Newton's law of Cooling
i> From Lecture 33-2, the PDE of the problem is given by:-


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

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

\frac{\partial }\left( {{K}_{ij}}\frac{\partial u} \right)+f\left( x,t \right)=\rho c\frac{\partial u}{\partial t}

$$ (Eq )<1>
 * 
 * }

From the PDE the weak form can be written as:-


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

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

\int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right)+f\left( x,t \right)-\rho c\frac{\partial u}{\partial t} \right)}d\text{ }\!\!\Omega\!\!\text{ =0}

$$ (Eq )<2>
 * 
 * }

Now from the theory of Integration by parts, we can write


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

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

\int{w\left( \frac{\partial }{\partial {{x}_{i}}}\left( {{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}} \right) \right)}d\text{ }\!\!\Omega\!\!\text{ =}\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}}}{{K}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}d}\text{ }\!\!\Omega\!\!\text{ }

$$
 * }

Hence equation 2 can be now written as:-


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

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

\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

$$ (Eq )<3>
 * 
 * }

Now we can write $$\displaystyle d\Omega $$ as $$\displaystyle d\Omega ={{\Gamma }_{g}}\cup {{\Gamma }_{h}}\cup {{\Gamma }_{H}}$$, where $$\displaystyle {{\Gamma }_{g}}$$ = essential boundary condition. $$\displaystyle {{\Gamma }_{h}}$$ = natural boundary condition. $$\displaystyle {{\Gamma }_{H}}$$ = Mixed boundary condition.

So by applying the Gauss theorem on the first term we obtain,


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

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

\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

$$ (Eq )<4>
 * 
 * }

Now the heat flux is given by:-


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

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

{{q}_{i}}=-{{\text{ }\!\!\Kappa\!\!\text{ }}_{ij}}\frac{\partial u}{\partial {{x}_{j}}}

$$ (Eq )<5>
 * 
 * }

So by rearranging Equation 5, we can write


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

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

\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

$$ (Eq )<6>
 * 
 * }

We can write the first term as-

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

We select $$\displaystyle w $$ such that $$\displaystyle w=0 $$ on $$\displaystyle {{\Gamma }_{g}}$$

On $$\displaystyle {{\Gamma }_{h}},\quad  q(x,t)n(x)=h(x,t) ,\forall x\in {{\Gamma }_{h}}$$

$$\displaystyle {{\Gamma }_{H}},\quad q(x,t)n(x)=H[u(x,t)-{u}_{\infty}] ,\forall x\in {{\Gamma }_{H}}$$

So we can write Equation 6 as :-


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

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

-\int\limits_{whd{{\text{ }\!\!\Gamma\!\!\text{ }}_{h}}}-\int\limits_{wH(u-{{u}_{\infty }})d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}}-\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_{\Omega }{wf\left( x,t \right)d\text{ }\!\!\Omega\!\!\text{ }-}\int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ }=0}

$$ (Eq )<7>
 * 
 * }

This can be re-arranged in the form:-


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

$$\displaystyle \int\limits_{\Omega }{w\rho c\frac{\partial u}{\partial t}d\Omega +\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega +}}\int\limits_{wHu}d{{\Gamma }_{H}}=-\int\limits_{wh}d{{\Gamma }_{h}}+\int\limits_{wH{{u}_{\infty }}}d{{\Gamma }_{H}}+\int\limits_{\Omega }{wf(x,t)}d\Omega $$ (Eq )<8>
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * 
 * }

This Equation is of the form:-


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

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

\tilde{m}(w,u)+\tilde{k}(w,u)=\tilde{f}(w) \quad \forall w $$ such that $$\displaystyle w=0 \quad on \quad {{\Gamma }_{g}}

$$
 * }

where,


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

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

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

$$\displaystyle\tilde{k}(w,u)=\int\limits_{\Omega }{\nabla w.\kappa .\nabla ud\Omega }+\int\limits_{wHu}d{{\text{ }\!\!\Gamma\!\!\text{ }}_{H}}$$

$$\displaystyle\tilde{f}(w)=-\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{ }}

$$ (Eq )<9>
 * 
 * }

Developing Semi-Discrete equations
ii>

To form the semi-discrete equations, we take $$\displaystyle w \quad and \quad u$$ as:-


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

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

w^h=\sum\limits_{i} \quad and \quad u^h=\sum\limits_{j}

$$ (Eq )<10>
 * 
 * }

Accordingly the mass, capacitance and heat source equations can be written as:-

Capacitance matrix


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

$$\displaystyle \overset{\tilde{\ }}{\mathop{m}}\,(w^h,u^h)=\int\limits_{\Omega }{w^h\rho c\frac{\partial u^h}{\partial t}d\Omega } $$
 * style="width:95%" |
 * style="width:95%" |

$$\displaystyle=\int\limits_{\Omega }{\sum{{{c}_{i}}{{b}_{i}}\rho c\sum{{{d}_{j}}{{b}_{j}}d\Omega }}}$$

$$\displaystyle=\sum\limits_{i}{\sum\limits_{j}{{{c}_{i}}\left( \int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega \right)}}d_{j}^{'} $$
 * }


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

$$\displaystyle \overset{\tilde{\ }}{\mathop{M}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{{{b}_{i}}\rho c{{b}_{j}}}d\Omega $$ (Eq )<11>
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * <p style="text-align:right">
 * }

Conductivity matrix


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

$$\displaystyle k(w^h,u^h)=\int\limits_{\Omega }{\nabla w^h.\kappa .\nabla u^hd\Omega +\int\limits_{w^hHu^hd{{\Gamma }_{H}}}} $$
 * style="width:95%" |
 * style="width:95%" |

$$\displaystyle=\int\limits_{\Omega }{\frac{\partial \sum}{\partial {{x}_{i}}}}.\kappa .\frac{\partial \sum}{\partial {{x}_{j}}}d\Omega +\int\limits_{\sumH\sumd{{\Gamma }_{\Eta }}} $$

$$\displaystyle=\sum\limits_{i}{\sum\limits_{j}\int\limits_{\Omega }{\left( \frac{\partial {{b}_{i}}}{\partial {{x}_{i}}}.\kappa \frac{\partial {{b}_{j}}}{\partial {{x}_{j}}}d\Omega \right)}}{{d}_{j}}+\sum\limits_{i}{\sum\limits_{j}\int\limits_{\left( {{b}_{i}}H{{b}_{j}}d{{\Gamma }_{\Eta }} \right)}}{{d}_{j}}

$$
 * }


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

$$\displaystyle \overset{\tilde{\ }}{\mathop{K}}\,({{b}_{i}},{{b}_{j}})=\int\limits_{\Omega }{b_{i}^{'}\kappa b_{j}^{'}d\Omega +\int\limits_{{{b}_{i}}H{{b}_{j}}d}}{{\Gamma }_{\Eta }} $$ (Eq )<12>
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * <p style="text-align:right">
 * }

Force matrix


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

$$\displaystyle \overset{\tilde{\ }}{\mathop{f}}\,(w^h)=-\int\limits_{w^hhd{{\Gamma }_{h}}}+\int\limits_{w^hH{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{w^hf(x,t)d\Omega }$$
 * style="width:95%" |
 * style="width:95%" |

$$\displaystyle =-\int\limits_{\sum\limits_{i}hd{{\Gamma }_{h}}}+\int\limits_{\sum\limits_{i}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{\sum\limits_{i}f(x,t)d\Omega }$$

$$\displaystyle =-\sum\limits_{i}\int\limits_{{{b}_{i}}hd{{\Gamma }_{h}}}+\sum\limits_{i}\int\limits_{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\sum\limits_{i}\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega } $$
 * }


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

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

\overset{\tilde{\ }}{\mathop{F}}\,(b_i)=-\int\limits_{{{b}_{i}}hd{{\Gamma }_{h}}}+\int\limits_{{{b}_{i}}H{{u}_{\infty }}d{{\Gamma }_{H}}}+\int\limits_{\Omega }{{{b}_{i}}f(x,t)d\Omega }

$$ (Eq )<13>
 * <p style="text-align:right">
 * }

3D Contour plots
iii>

consider the case of n=m=3 the number of nodes is 9 The trial solution is
 * {| style="width:100%" border="0"

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

u_{p}^{h}(x)={{d}_{1}}*N_{1}+{{d}_{2}}N_{2}+{{d}_{3}}N_{3}+\cdots +{{d}_{9}}N_{9},_ – ^ – j=0,1,2,\cdots ,9

$$
 * }

Since the trial solution must satisfy the essential boundary condition,, we should have  $${{d}_{1}}=2,{{d}_{2}}=2,{{d}_{3}}=2,{{d}_{6}}=2,{{d}_{9}}=2$$., because it is specified that through out the boundary the temperature is 2 i.e for all the nodes on the boundary temperature is 2

In order to solve for the rest coefficient $${{d}_{j}},_ – ^ – j=1,2,\cdots $$, we construct the Conductance matrix, Heat Source matrix and Capacitance matrix as given by Equations 12,13 and 14.

Using Matlab to construct above matrices until satisfying the convergence criterion, we finally have
 * {| style="width:100%" border="0"

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

=\left( \begin{matrix} 151/200 &  -1/5 &  -1/30 & -133/1000 & -16/45 &    1/9 & -33/500 &    1/9 &  -1/45\\ -1/5 &  88/45 &   -1/5 & -16/45 & -16/15 & -16/45 & 1/9 & 0 & 1/9\\ -1/30 & -1/5 & 28/45 & 1/9 & -16/45 & -1/5 &  -1/45 & 1/9 & -1/30\\ -133/100 & -16/45 & 1/9 & 112/45 & -16/15 & 0 & -133/1000 & -16/45 & 1/9\\ -16/45 & -16/15 & -16/45 & -16/15  & 256/45 & -16/15 & -16/45 & -16/15 & -16/45\\ 1/9 & -16/45 & -1/5 & 0 & -16/15 & 88/45 & 1/9 & -16/45 & -1/5\\ -33/500 & 1/9 & -1/45 & -133/1000 & -16/45 & 1/9 & 1511/2000 & -1/5 & -1/30\\ 1/9 & 0 & 1/9 & -16/45 & -16/15 & -16/45 & -1/5 & 88/45 & -1/5\\ -1/45 & 1/9 & -1/30 & 1/9 & -16/45 & -1/5 & -1/30 & -1/5 & 28/45\\ \end{matrix} \right)

$$
 * }


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

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

{F}=\left( \begin{matrix} 166/1000\\ 0\\ 0\\ 33/500\\ 0\\ 0\\ -59/60\\ -4\\ -1\\ \end{matrix} \right)

$$
 * }


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

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

{d}=\left( \begin{matrix} 2.000\\ 2.000\\ 2.000\\ 0.335\\ 0.847\\ 2.000\\ -1.460\\ -1.331\\ 2.000\\ \end{matrix} \right)

$$
 * }

The trial solution using 2D LIBF obtained is as follows:-


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

$$\displaystyle \begin{align} &U_{F}^{h}=(x)0.50*x*y*(x - 1)*(y - 1) - 1*x*(1.0*x + 1.0)*(1.0*y + 1.0)*(y - 1)\\& \quad + 0.665*y*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1) - 0.1678*x*(1.0*y + 1.0)*(x - 1)*(y - 1)\\& \quad - 1.000*y*(1.0*x + 1.0)*(x - 1)*(y - 1) + 0.847*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1)*(y - 1) \\& \quad + 0.500*x*y*(1.0*x + 1.0)*(1.0*y + 1.0) + 0.500*x*y*(1.0*x + 1.0)*(y - 1) - 0.365*x*y*(1.0*y + 1.0)*(x - 1)\\ \end{align} $$ (Eq )<14>
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * <p style="text-align:right">
 * }

The Matlab code to generate the plots and LIBF is:

The 3D contour plot generated from MATLAB is dislpayed below.



=Contributing Members & Referenced Lecture=