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

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>
 * 
 * }

Non-discrete

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 } $$ $$\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:95%" |
 * style="width:95%" |
 * }
 * {| 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> Conductivity matrix
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * 
 * }


 * {| 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}}}} $$ $$\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:95%" |
 * style="width:95%" |

$$
 * }


 * {| 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> Force matrix
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * 
 * }


 * {| 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 }$$ $$\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:95%" |
 * style="width:95%" |
 * }


 * {| 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>
 * 
 * }

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" |
 * 
 * }

The Matlab code to generate the plots and LIBF is:

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