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

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

\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{w\rho c\frac{\partial u}{\partial t}d\text{ }\!\!\Omega\!\!\text{ +}}\int\limits_{\text{ }\!\!\Omega\!\!\text{ }}{\nabla w\left( {{K}_{ij}}\nabla u \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{ }}

$$ (Eq )<8>
 * 
 * }

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_{\text{ }\!\!\Omega\!\!\text{ }}{\nabla w\left( {{K}_{ij}}\nabla u \right)d\text{ }\!\!\Omega\!\!\text{ }}+\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>

Let $$\displaystyle u=\left[ N \right]{{u}_{e}}\quad and \quad w=\left[ N \right]{{w}_{e}},$$

where $$\displaystyle [N] $$ is the vector of interpolation functions (Lagrange polynomial) and $$\displaystyle {u}_{e}\quad and \quad {w}_{e} $$ are the nodal values of $$\displaystyle u $$ and $$\displaystyle w $$ respectively.

Now by taking the derivative of basis functions with respect to $$\displaystyle {x}_{1}$$ and $$\displaystyle {x}_{2}$$, 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" |

\left[ \begin{align} & \frac{dN}{d{{x}_{1}}} \\ & \frac{dN}{d{{x}_{2}}} \\ \end{align} \right]=\left[ B \right]

$$ (Eq )<10>
 * 
 * }

Substituting the Lagrangian forms in Equation 9, we obtain:-


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

$$\displaystyle \int\limits_{\Omega }{{{w}_{e}}^{T}\rho c{{[N]}^{T}}[N]\frac{\partial u}{\partial t}d\Omega }+\int\limits_{\Omega }{{{w}_{e}}^{T}{{[B]}^{T}}[B]{{u}_{e}}d\Omega }+\int\limits_{{{w}_{e}}^{T}H{{[N]}^{T}}[N]{{u}_{e}}d{{\Gamma }_{h}}=-}\int\limits_{{{w}_{e}}^{T}h{{[N]}^{T}}d{{\Gamma }_{h}}}+\int\limits_{{{w}_{e}}^{T}H{{[N]}^{T}}{{u}_{\infty }}d{{\Gamma }_{h}}}+\int\limits_{\Omega }{w_{e}^{T}{{[N]}^{T}}f(x,t)d\Omega } $$ (Eq )<11>
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * 
 * }

From the above equation the conductivity matrix, the heat source vector and the capacitance matrix can be described as:- Capacitance matrix


 * {| 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)= \int\limits_{\Omega }{\rho c{{[N]}^{T}}[N]\frac{\partial u}{\partial t}d\Omega } $$
 * }

Conductivity matrix


 * {| 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{k}(w,u)= \int\limits_{\Omega }{{{[B]}^{T}}[B]d\Omega }+\int\limits_{H{{[N]}^{T}}[N]d{{\Gamma }_{h}}} $$
 * }

Heat source matrix


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

plotting
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" |

{{K}_{FF}}=\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 )<15>
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * style="width:92%; padding:10px; border:2px solid #000000" |
 * 
 * }