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

 Homework 6 =Problem 6.1 Using Quadratic Lagrange Element Basis Functions to Obtain Trial and Exact Solution of Torsion Bar=

Given: Data set G1DM1.0/D1b defined on Lecture 26-2
Strong form for the Torsion Bar

where $${a_2}\left( x \right){\text{ =  2 + 3x}}$$  $$\qquad\text{and}\quad $$   $$f\left( {x,t} \right) = 5x $$  so we get Essential boundary conditions, Natural boundary conditions,

Find

1. For nel=2 compute  $$\tilde{\mathbf{K}}=\sum_{e=1}^{2}\tilde{\mathbf{K^{e}}}$$   with   $$\tilde{\mathbf{K^{e}}}$$   by Lecture 30-5 display  $$\tilde{\mathbf{K^{e}}}, e=1,2 $$

2. Compute $$\mathbf{k^{e}}$$ ,  $$\mathbf{L^{e}}$$  for e=1,2

3. Compute $$ \tilde{K}^{e}=L^{eT}K^{e}L^{e}$$ for e=1,2 and compare with 1)

4. Plot all QLEBF for nel=3

5. Plot $$ u_{\tilde{n}}^{h}  \ vs\   u$$ and $$ u_{\tilde{n}}^{h}\left(0.5 \right)-u\left(0.5 \right)\ vs\ \tilde{n}$$

1) Compute $$\tilde{\mathbf{K}}$$  for nel=2
$$\tilde{\mathbf{K^{e}}}$$ can be written as

(Refer Lecture 30-5)

the above figure shows a torsion bar with 2 elements, $$\omega _{1}$$  and  $$\omega _{2}$$  ,since we are considering a QLEBF, number of elements per node is 3. The nodes are marked in the above figure as points  $$x_{1},x_{2},x_{3},x_{4} and x_{5}$$  since the length of the bar is 1, and we have taken equidistant element nodes, the nodes are located at (0,0.25,0.5,0.75,1) respectively.

Since there are 3 nodes per element there will be three basis functions per element corresponding to each of the nodes.

The basis functions for the 1st element can be written as

The basis functions for the 2nd element can be written as

From Eq(1.6) $$\tilde{\mathbf{K}}\mathbf{_{11}^{1}}=\int_{x_{1}}^{x_{3}}{b_{1}^{1}}'(x)a_{2}(x){b_{1}^{1}}'(x)dx$$

Substituting all the terms and integrating we get $$\tilde{\mathbf{K}}\mathbf{_{11}^{1}}=10.833$$

Similarly the other terms in the  $$\tilde{\mathbf{K}}\mathbf{_{ij}^{1}}$$  can be found out and using MATLAB (MATLAB code attached)

Similarly

Commment
Note the order of $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  is  $$\tilde{n} X \tilde{n}$$  where  $$\tilde{n} $$  is the total number of degrees of freedom,  $$\tilde{n}=n_{E}+n_{F} $$ , (  $$n_{E} $$  no. of prescribed dofs on ess. bc and  $$n_{F} $$  no. of free dofs) in this case for nel=2  $$n_{E}=2 $$  and  $$n_{F}=3 $$  , hence the order of  $$\tilde{\mathbf{K}}\mathbf{^{e}}=5 X 5$$

2) Compute $$\mathbf{k^{e}},\mathbf{L^{e}}$$  for e=1,2
In this case order of $$\mathbf{k^{e}}$$  is 3 X 3 and order of  $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  is 5 X 5 hence for the matrix multiplication to be possible the order of  $$\mathbf{L^{e}}$$  is 3 X 5

Hence for e=1,2

$$\mathbf{k^{e}}$$ can be calculated from Eq 1.6 similar to  $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  but the order of  $$\mathbf{k^{e}}$$  is 3 X 3 hence for element 1

Similarly for element 2

$$\mathbf{L^{1}}$$ for the 1st element is given by

for element 2 $$\mathbf{L^{2}}$$  is given by

3) Compute $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  for e=1,2
for element 1

For element 2

The $$\tilde{\mathbf{K}}\mathbf{^{1}}$$  and  $$\tilde{\mathbf{K}}\mathbf{^{2}}$$  matrix obtained from part 1) i.e from Eq. 1.13 and Eq. 1.14 are same as that obtained in part 3 i.e Eq. 1.22 and Eq. 1.23

5)Plot $$ u_{\tilde{n}}^{h} \ vs\   u$$ and $$ u_{\tilde{n}}^{h}\left(0.5 \right)-u\left(0.5 \right)\ vs\ \tilde{n}$$

=Problem 6.2 Using Quadratic Lagrange Element Basis Functions to Obtain Trial and Exact Solution with Data Set G1DM1.0/D1=

Given: Data set G1DM1.0/D1 defined on [[media:fe1.s11.mtg9.djvu|Mtg 9 (e)]]
Strong form

where $${a_2}\left( x \right){\text{ =  2 + 3x}}$$  $$\qquad\text{and}\quad $$   $$f\left( {x,t} \right) = 5x $$  so we get Essential boundary conditions, Natural boundary conditions,

Find

1. For nel=2 compute  $$\tilde{\mathbf{K}}=\sum_{e=1}^{2}\tilde{\mathbf{K^{e}}}$$   with   $$\tilde{\mathbf{K^{e}}}$$   by Lecture 30-5 display  $$\tilde{\mathbf{K^{e}}}, e=1,2 $$

2. Compute $$\mathbf{k^{e}}$$ ,  $$\mathbf{L^{e}}$$  for e=1,2

3. Compute $$ \tilde{K}^{e}=L^{eT}K^{e}L^{e}$$ for e=1,2 and compare with 1)

4. Plot all QLEBF for nel=3

5. Plot $$ u_{\tilde{n}}^{h}  \ vs\   u$$ and $$ u_{\tilde{n}}^{h}\left(0.5 \right)-u\left(0.5 \right)\ vs\ \tilde{n}$$

1) Compute $$\tilde{\mathbf{K}}$$  for nel=2
$$\tilde{\mathbf{K^{e}}}$$ can be written as

(Refer Lecture 30-5)

the above figure shows a torsion bar with 2 elements, $$\omega _{1}$$  and  $$\omega _{2}$$  ,since we are considering a QLEBF, number of elements per node is 3. The nodes are marked in the above figure as points  $$x_{1}, \quad x_{2},\quad x_{3},\quad x_{4}\quad  and \quad x_{5}$$  since the length of the bar is 1, and we have taken equidistant element nodes, the nodes are located at (0, 0.25, 0.5, 0.75, 1) respectively.

Since there are 3 nodes per element there will be three basis functions per element corresponding to each of the nodes.

The basis functions for the 1st element can be written as

The basis functions for the 2nd element can be written as

From Eq(2.6) $$\tilde{\mathbf{K}}\mathbf{_{11}^{1}}=\int_{x_{1}}^{x_{3}}{b_{1}^{1}}'(x)a_{2}(x){b_{1}^{1}}'(x)dx$$

Substituting all the terms and integrating we get $$\tilde{\mathbf{K}}\mathbf{_{11}^{1}}=10.833$$

Similarly the other terms in the  $$\tilde{\mathbf{K}}\mathbf{_{ij}^{1}}$$  can be found out and using MATLAB (MATLAB code attached)

Similarly

Commment
Note the order of $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  is  $$\tilde{n} X \tilde{n}$$  where  $$\tilde{n} $$  is the total number of degrees of freedom,  $$\tilde{n}=n_{E}+n_{F} $$ , (  $$n_{E} $$  no. of prescribed dofs on ess. bc and  $$n_{F} $$  no. of free dofs) in this case for nel=2  $$n_{E}=2 $$  and  $$n_{F}=3 $$  , hence the order of  $$\tilde{\mathbf{K}}\mathbf{^{e}}=5 X 5$$

2) Compute $$\mathbf{k^{e}},\mathbf{L^{e}}$$  for e=1,2
In this case order of $$\mathbf{k^{e}}$$  is 3 X 3 and order of  $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  is 5 X 5 hence for the matrix multiplication to be possible the order of  $$\mathbf{L^{e}}$$  is 3 X 5

Hence for e=1,2

$$\mathbf{k^{e}}$$ can be calculated from Eq 1.6 similar to  $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  but the order of  $$\mathbf{k^{e}}$$  is 3 X 3 hence for element 1

Similarly for element 2

$$\mathbf{L^{1}}$$ for the 1st element is given by

for element 2 $$\mathbf{L^{2}}$$  is given by

3) Compute $$\tilde{\mathbf{K}}\mathbf{^{e}}$$  for e=1,2
for element 1

For element 2

The $$\tilde{\mathbf{K}}\mathbf{^{1}}$$  and  $$\tilde{\mathbf{K}}\mathbf{^{2}}$$  matrix obtained from part 1) i.e from Eq. 2.13 and Eq. 2.14 are same as that obtained in part 3 i.e Eq. 2.22 and Eq. 2.23

5)Plot $$ u_{\tilde{n}}^{h} \ vs\   u$$ and $$ u_{\tilde{n}}^{h}\left(0.5 \right)-u\left(0.5 \right)\ vs\ \tilde{n}$$

=Problem 6.3 Reproduce all steps in tutorial by Beede =

Solution
The PDF of the explanation of all the examples of the tutorial by Beede is attached below.

PDF file of the Problem 6.3

=Problem 6.4= =Problem 6.4(a) Verify the Divergence Theorem=

Problem Statement
Given a vector field $$\;\;q_{x} = -y^2\;;\; q_{y} = -2xy \;\;$$ on the domain shown below. Verify the divergence theorem.

Solution
$$ \vec{q} = q_{x}\vec{i}+ q_{y}\vec{j}$$

$$\therefore\;\; \vec{q} = -y^2\vec{i} - 2xy\vec{j}$$

Divergence of $$\vec{q} = \vec{\triangledown }. \vec{q} = = \frac{\partial q_{x}}{\partial x}+\frac{\partial q_{y}}{\partial y} = = 0 - 2x = -2x$$ The divergence theorem states,

$$\int_{\Omega }\vec{\triangledown}.\vec{q} = \oint_{\Gamma }\vec{q}.\vec{n}$$

$$\oint_{\Gamma }\vec{q}.\vec{n} = \int_{AB}\vec{q}.\vec{n^{1}}\;d\Gamma_{AB}+\int_{BC}\vec{q}.\vec{n^{2}}\;d\Gamma_{BC}+\int_{CD}\vec{q}.\vec{n^{3}}\;d\Gamma_{CD}+\int_{DA}\vec{q}.\vec{n^{4}}\;d\Gamma_{DA}$$

$$\vec{q}.\vec{n^{1}} = 2xy\;;\vec{q}.\vec{n^{2}} = -y^2\;;\vec{q}.\vec{n^{3}} = -2xy\;;\vec{q}.\vec{n^{4}} = y^2$$

$$d\Gamma_{AB} = dx;\;d\Gamma_{BC} = dy;\;d\Gamma_{CD} = -dx;\;d\Gamma_{DA} = -dy$$

From equations (1) & (2), we have

$$\int_{\Omega }\vec{\triangledown}.\vec{q} = \oint_{\Gamma }\vec{q}.\vec{n}$$

Hence, Divergence theorem is verified

=Problem 6.4(b) Verify the Divergence Theorem=

Problem Statement
Given a vector field $$\;\;q_{x} = 3x^2y+y^3\;;\; q_{y} = 3x+y^3 \;\;$$ on the domain shown below. Verify the divergence theorem.The curved boundary of the domain is a parabola.

Solution
$$ \vec{q} = q_{x}\vec{i}+ q_{y}\vec{j}$$

$$\therefore\;\; \vec{q} = (3x^2y+y^3)\vec{i}+(3x+y^3) \vec{j}$$

Divergence of $$\vec{q} = \vec{\triangledown }. \vec{q} = \frac{\partial q_{x}}{\partial x}+\frac{\partial q_{y}}{\partial y} = 6xy+3y^2$$

The divergence theorem states,

$$\int_{\Omega }\vec{\triangledown}.\vec{q} = \oint_{\Gamma }\vec{q}.\vec{n}$$

The equation of the given parabola is $$x^2 = 4-y \;\;\Rightarrow \;\; y = 4-x^2$$

$$\int_{\Omega }\vec{\triangledown}.\vec{q} = \int_{-2}^{2}\left ( \int_{0}^{4-x^2}(6xy+3y^2)\;dy\right )dx = \int_{-2}^{2}\left (  \int_{0}^{4-x^2}(6x(4-x^2)+3(4-x^2)^2)\;dy\right )dx = \int_{-2}^{2}\left ( -(x^2-4)^2(x^2-3x-4)\right )dx$$

$$\oint_{\Gamma }\vec{q}.\vec{n} = \int_{AB}\vec{q}.\vec{n^{1}}\;d\Gamma_{AB}+\int_{BCA}\vec{q}.\vec{n^{2}}\;d\Gamma_{BCA}$$

$$\vec{q}.\vec{n^{1}} = -(3x+y^3)\;;\vec{q}.\vec{n^{2}} = 3x+y^3$$

$$d\Gamma_{AB} = dx;\;d\Gamma_{BCA} = dy$$

Equation of line AB is $$y = 0$$

$$\therefore\;\;\int_{AB}\vec{q}.\vec{n^{1}}\;d\Gamma_{AB} = \int_{-2}^{2}-3xdx = 0$$

$$\int_{BCA}\vec{q}.\vec{n^{2}}\;d\Gamma_{BCA} = \int_{-2}^{2}3x+y^3$$

Equation of parabola is $$y = 4 - x^2$$

, using [| Wolframalpha]

From equations (1) & (2), we have

$$\int_{\Omega }\vec{\triangledown}.\vec{q} = \oint_{\Gamma }\vec{q}.\vec{n}$$

Hence, Divergence theorem is verified

=Problem 6.4(c) Use Divergence Theorem=

Problem Statement
Using the divergence theorem prove

$$\oint_{\Gamma}n\;\;\mathrm{d}\Gamma = 0$$

Solution
Let, $$\vec{q} = \vec{i}+\vec{j}\;\;\; and \;\; \vec{n} = \vec{n_{1}}+\vec{n_{2}}$$

Divergence Theorem states,

$$\int_{\Omega }\vec{\triangledown}.\vec{q}\mathrm{d}\Omega = \oint_{\Gamma }\vec{q}.\vec{n}\mathrm{d}\Gamma$$

$$\oint_{\Gamma }\vec{q}.\vec{n}\;\mathrm{d}\Gamma = \oint_{\Gamma }(n_{1}+n_{2})\;\mathrm{d}\Gamma = \oint_{\Gamma }n\;\mathrm{d}\Gamma$$

$$\therefore\;\;\oint_{\Gamma }n\;\mathrm{d}\Gamma = \int_{\Omega }\vec{\triangledown}.\vec{q}\mathrm{d}\Omega $$

$$\vec{\triangledown}.\vec{q} = \frac{\partial (1)}{\partial x}+\frac{\partial (1)}{\partial y} = 0$$

$$\therefore\; \int_{\Omega }\vec{\triangledown}.\vec{q}\mathrm{d}\Omega = \int_{\Omega }0\;\mathrm{d}\Omega = 0$$

$$\Rightarrow\;\;\oint_{\Gamma}n\;\;\mathrm{d}\Gamma = 0$$

=Problem 6.5 Using 2DLIBF to Obtain Trial and Exact Solution for heat problem with boundaries enclosed by a bi-unit square=

Given: Data set G2DM1.0/D1 defined on Lecture 35-4
Strong form for the heat problem is given by Lecture 33-2

where $$\boldsymbol{\kappa}\left( {x,t} \right) = \mathbf{I} $$ $$\qquad\text{and}\quad $$   $$f\left( {x,t} \right) = 0 $$   and $$\frac{\partial u}{\partial t}\left( {x,t} \right) = 0 $$

so we get

Essential boundary conditions, Natural boundary conditions,

Find
1. Use 2DLIBF with n=m=2,4,6,8... to solve G2DM1.0/D1 for $$u^{h}$$ until accuracy of $$10^{-6}$$ at centre (x,y)=(0,0)

1) Finding Exact solution
The solution for the PDE $$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0$$ will be of the form

$$u(x,y)=c_{1}xy+c_{2}x+c_{3}y+c_{4}$$ where,  $$c_{1},c_{2},c_{3} and c_{4}$$  are constants.

The essential BC is specified through out the boundary of the square hence u=2 for the points on the 4 vertices of the square. i.e for (-1,-1);(-1,1);(1,-1) and (1,1) we have u=2, substituting for these four points into the equation for u above and solving them simultaneously we get  $$c_{1}=0$$ ,  $$c_{2}=0$$  ,  $$c_{3}=0$$  and  $$c_{4}=2$$ hence the exact solution of the PDE is given by

Hence, the trial solution of the governing equation is

2) 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

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}_{4}}=2,{{d}_{6}}=2,{{d}_{7}}=2,{{d}_{8}}=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 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 $$
 * }

Substituting Eq 6.1 and 6.3 into the above Equation we get,

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

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

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

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

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

\tilde{f}(\omega)=0 $$
 * }

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

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

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

Using Matlab to construct above two 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" |

= \begin{bmatrix} 0.6222 & -0.2000  & -0.0333 &  -0.2000 &  -0.3556  &  0.1111 & -0.0333  &  0.1111 &  -0.0222 \\   -0.2000 &  1.9556  & -0.2000 &  -0.3556 &  -1.0667  & -0.3556 &   0.1111 &        0&    0.1111\\   -0.0333 &  -0.2000 &   0.6222&    0.1111&   -0.3556 &  -0.2000&   -0.0222&    0.1111&   -0.0333\\   -0.2000 &  -0.3556 &   0.1111&    1.9556&   -1.0667 &        0&   -0.2000&   -0.3556&    0.1111\\   -0.3556 &  -1.0667 &  -0.3556&   -1.0667&    5.6889 &  -1.0667&   -0.3556&   -1.0667&   -0.3556\\    0.1111 &  -0.3556 &  -0.2000&         0&   -1.0667 &   1.9556&    0.1111&   -0.3556&   -0.2000\\   -0.0333 &   0.1111 &  -0.0222&  -0.2000 &  -0.3556  &  0.1111 &   0.6222 &  -0.2000 &  -0.0333\\    0.1111 &        0 &   0.1111&   -0.3556&   -1.0667 &  -0.3556&   -0.2000&    1.9556&   -0.2000\\   -0.0222 &   0.1111 &  -0.0333&    0.1111&   -0.3556 &  -0.2000&   -0.0333&   -0.2000&    0.6222\\

\end{bmatrix}

$$
 * }

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

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

=\begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}

$$ By the equation that $$ Kd=F  $$, we can solve the rest of $$d$$ 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" |

=\begin{bmatrix} 2\\ 2\\ 2\\ 2\\ 2\\ 2\\ 2\\ 2\\ 2 \end{bmatrix}

$$
 * }

Hence, the trial solution of the governing equation is
 * {| style="width:100%" border="0"

$$\displaystyle \begin{align} &U_{F}^{h}(x)= 0.5*x*y*(x - 1)*(y - 1) - 1.0*x*(1.0*x + 1.0)*(1.0*y + 1.0)*(y - 1)\\ & \quad - 1.0*y*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1) - 1.0*x*(1.0*y + 1.0)*(x - 1)*(y - 1)\\& \quad - 1.0*y*(1.0*x + 1.0)*(x - 1)*(y - 1) + 2.0*(1.0*x + 1.0)*(1.0*y + 1.0)*(x - 1)*(y - 1)\\& \quad + 0.5*x*y*(1.0*x + 1.0)*(1.0*y + 1.0) + 0.5*x*y*(1.0*x + 1.0)*(y - 1)\\& \quad + 0.5*x*y*(1.0*y + 1.0)*(x - 1) \\ \end{align}
 * style="width:92%; padding:10px; border:2px solid #ff0000" |
 * style="width:92%; padding:10px; border:2px solid #ff0000" |

$$
 * }

It is found that the exact and actual solution are exactly matching for n=m=3, the plot of exact and trial solution is as shown below

=Problem 6.6 Solve G2DM1.0/D1 using 2D LLEBF until accuracy 10^{-6} at center(0,0)=

Problem Statement
Refer to lecture notes [[media:fe1.s11.mtg35.djvu|35-4]], [[media:fe1.s11.mtg36.djvu|36-6]] and [[media:fe1.s11.mtg37.djvu|37-1]] for more detailed description.

Given
Strong form for the heat problem is given in Lecture 33-2:

where $$\boldsymbol{\kappa}\left( {x,t} \right) = \mathbf{I} $$ $$\qquad\text{and}\quad $$   $$f\left( {x,t} \right) = 0 $$   and $$\frac{\partial u}{\partial t}\left( {x,t} \right) = 0 $$

so we get

Essential boundary conditions, Natural boundary conditions,

2). Integrate element stiffness matrix in parent coordinate
Find appropriate u to integrate stiffness matrix exactly with Gauss-Legendre Quadrature. Given detailed argument.

3). Use distorted quadrilateral meshes, compare results with uniform meshes
Meshes can be generated from CGX, and analysis can be done using Fish and Belytschko's code and CCX.

Find Exact solution
The solution for the PDE $$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0$$ will be of the form

$$u(x,y)=c_{1}xy+c_{2}x+c_{3}y+c_{4}$$ where,  $$c_{1},c_{2},c_{3} and c_{4}$$  are constants.

Details of the solution of PDE can be find at EqWorld

The essential BC is specified through out the boundary of the square hence u=2 for the points on the 4 vertices of the square, i.e for (-1,-1);(-1,1);(1,-1) and (1,1) we have u=2.

Substituting for these four points into the equation for u above and solving them simultaneously we get $$c_{1}=0$$ ,

$$c_{2}=0$$ ,

$$c_{3}=0$$

and $$c_{4}=2$$

hence the exact solution of the PDE is given by

Results and plots
Below is the analysis results:

=Problem 6.7=

Problem Statement
Refer to lecture notes [[media:fe1.s11.mtg36.djvu|36-1,2]] for more detailed description.

Given

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

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

A=\left[ \begin{matrix} 1 & 1 & 1 \\   2 & -1 & 3  \\   3 & 2 & 6  \\ \end{matrix} \right], B=\left[ \begin{matrix} 1 & 3 & 5 \\   1 & -4 & 1  \\   2 & 5 & 8  \\ \end{matrix} \right]

$$
 * }

Objective
Prove the following two expressions using given A and B.

1). $${{\left( AB \right)}^{T}}={{B}^{T}}{{A}^{T}}$$

2). $${{\left( {{A}^{-1}} \right)}^{T}}={{\left( {{A}^{T}} \right)}^{-1}}={{A}^{-T}}$$

Solution
1). $${{\left( AB \right)}^{T}}={{B}^{T}}{{A}^{T}}$$
 * {| style="width:100%" border="0"

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

{{\left( AB \right)}^{T}}={{\left( \left[ \begin{matrix}  1 & 1 & 1  \\   2 & -1 & 3  \\   3 & 2 & 6  \\ \end{matrix} \right]\left[ \begin{matrix}   1 & 3 & 5  \\   1 & -4 & 1  \\   2 & 5 & 8  \\ \end{matrix} \right] \right)}^{T}}={{\left[ \begin{matrix} 4 & 4 & 14 \\   7 & 25 & 33  \\   17 & 31 & 65  \\ \end{matrix} \right]}^{T}}=\left[ \begin{matrix} 4 & 7 & 17 \\   4 & 25 & 31  \\   14 & 33 & 65  \\ \end{matrix} \right]

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

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

{{B}^{T}}{{A}^{T}}={{\left[ \begin{matrix} 1 & 3 & 5 \\   1 & -4 & 1  \\   2 & 5 & 8  \\ \end{matrix} \right]}^{T}}{{\left[ \begin{matrix} 1 & 1 & 1 \\   2 & -1 & 3  \\   3 & 2 & 6  \\ \end{matrix} \right]}^{T}}=\left[ \begin{matrix} 1 & 1 & 2 \\   3 & -4 & 5  \\   5 & 1 & 8  \\ \end{matrix} \right]\left[ \begin{matrix} 1 & 2 & 3 \\   1 & -1 & 2  \\   1 & 3 & 6  \\ \end{matrix} \right]=\left[ \begin{matrix} 4 & 7 & 17 \\   4 & 25 & 31  \\   14 & 33 & 65  \\ \end{matrix} \right]

$$ Hence,
 * }
 * {| 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( AB \right)}^{T}}={{B}^{T}}{{A}^{T}}

$$
 * }.

2). $${{\left( {{A}^{-1}} \right)}^{T}}={{\left( {{A}^{T}} \right)}^{-1}}={{A}^{-T}}$$
 * {| style="width:100%" border="0"

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

{{\left( {{A}^{-1}} \right)}^{T}}={{\left( {{\left[ \begin{matrix}  1 & 1 & 1  \\   2 & -1 & 3  \\   3 & 2 & 6  \\ \end{matrix} \right]}^{-1}} \right)}^{T}}=\frac{1}{8}{{\left[ \begin{matrix} 12 & 4 & -4 \\   3 & -3 & 1  \\   -7 & -1 & 3  \\ \end{matrix} \right]}^{T}}=\frac{1}{8}\left[ \begin{matrix} 12 & 3 & -7 \\   4 & -3 & -1  \\   -4 & 1 & 3  \\ \end{matrix} \right]

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

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

{{\left( {{A}^{T}} \right)}^{-1}}={{\left( {{\left[ \begin{matrix}  1 & 1 & 1  \\   2 & -1 & 3  \\   3 & 2 & 6  \\ \end{matrix} \right]}^{T}} \right)}^{-1}}={{\left( \left[ \begin{matrix}   1 & 2 & 3  \\   1 & -1 & 2  \\   1 & 3 & 6  \\ \end{matrix} \right] \right)}^{-1}}=\frac{1}{8}\left[ \begin{matrix} 12 & 3 & -7 \\   4 & -3 & -1  \\   -4 & 1 & 3  \\ \end{matrix} \right]

$$ So,
 * }
 * {| 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( {{A}^{-1}} \right)}^{T}}={{\left( {{A}^{T}} \right)}^{-1}}={{A}^{-T}}

$$
 * }

WolframAlpha code:

Multiplying two matrices:

{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}*{{1, 3, 5}, {1, -4, 1}, {2, 5, 8}}

Transposing:

transpose[{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}*{{1, 3, 5}, {1, -4, 1}, {2, 5, 8}}]

Inverse-Transpose:

inverse{transpose{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}}

Transpose-Inverse:

transpose[inverse{{1, 1, 1}, {2, -1, 3}, {3, 2, 6}}]

=Problem 6.8=

Problem Statement
Verify the table of {(wi,xi),i=1,2,..,5} against NIST handbook and also with table on textbook page 89

Solution
The table from Lecture 36-3 is given below:

Nodes and weights for the 5-point Gauss–Legendre formula from | NIST Handbook are given below:

The values of the Nodes and Weights in the given table in Lecture 36-3 are calculated and the values are tabulated:

The table of Nodes and Weights with the calculated values from the text book Introduction to Finite Element Analysis by Fish and Belytschko is as below:

From the above tables, it can be concluded that the values of Nodes and Weights from the table in Lecture 36-3, NIST handbiik and the text book agree with each other.

Problem Statement
Use Gauss quadrature to obtain exact values for the following integrals. Verify by analytical integration:

$$(a).\;\int_{0}^{4}(x^2+1)\mathrm{d}x$$

$$(b).\;\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

Solution
(a).

$$I =\int_{0}^{4}(x^2+1)\mathrm{d}x$$

As $$n_{gp} = 2\; (two-point \;\;quadrature),\; p = 2n_{gp} - 1 = 3$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \end{bmatrix} $$

Solving the above system of 4 equations and 4 unknowns, we get

$$W_{1} = W_{2} = 1\;\; and \; \xi_{1} = \frac{1}{\sqrt3}\;\;\xi_{2} = \frac{-1}{\sqrt3}$$

We have, $$a = 0\;\;;b = 4$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = 2+2\xi$$

The value of the function $$f(\xi) = (2+2\xi)^2+1$$

$$I = \frac{l}{2}\int_{-1}^{1}((2+2\xi)^2+1)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}((2+2\xi_{1})^2+1)\; + \;W_{2}((2+2\xi_{2})^2+1)\right] \;and \; l = b - a = 4 - 0 = 4$$

$$= 2\left[10.953+1.715\right] = 25.334$$

Analytical Method

$$I =\int_{0}^{4}(x^2+1)\mathrm{d}x$$

$$= \left[\frac{x^3}{3}+x\right]_{0}^{4} = \frac{4^3}{3}+4 = 25.333$$

The analytical method also yields exactly the same result as the Gauss Quadrature.

(b). $$I =\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

$$n_{gp} = 4\; (fo-point \;\;quadrature),\; p = 2n_{gp} - 1 = 7$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \xi_{1}^6&\xi_{2}^6 \\ \xi_{1}^7&\xi_{2}^7 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ W_{4}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \frac{2}{7}\\ 0\\ \end{bmatrix} $$

Solving the above system of 8 equations and 8 unknowns, we get

$$W_{1} = W_{3} = 0.348\;W_{2} = W_{4} = 0.652\; and \; \xi_{1} = 0.861\;\;\xi_{2} = 0.339\;\;\xi_{3} = -0.861\;\;\xi_{4} = -0.339$$

We have, $$a = -1\;\;;b = 1$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

The value of the function $$f(\xi) = (\xi^4+2\xi^2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(\xi^4+2\xi^2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(\xi_{1}^4+2\xi_{1}^2)\; + \;W_{2}(\xi_{2}^4+2\xi_{2}^2)\;+ (\xi_{3}^4+2\xi_{3}^2)\;+ (\xi_{4}^4+2\xi_{4}^2)\;\right] \;and \; l = b - a = 4 - 0 = 4$$

$$= 1.4144+0.3169 = 1.731$$

Analytical Method

$$I =\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

$$= \left[\frac{\xi^5}{5}+\frac{2\xi^3}{3}\right]_{-1}^{1} = 2\left[\frac{1}{5}+\frac{2}{3}\right] = 1.731$$

The analytical method also yields exactly the same result as the Gauss Quadrature.

Problem Statement
Use three-point Gauss quadrature to evaluate the following integrals. Compare to the analytical integral.

(a).$$ \int_{-1}^{1}\frac{\xi}{\xi^2+1}\;\mathrm{d}x $$

(b). $$\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi}$$

Solution
(a).

$$I =\int_{-1}^{1}\frac{\xi}{\xi^2+1}\;\mathrm{d}x $$

As $$n_{gp} = 3\; (three-point \;\;quadrature),\; p = 2n_{gp} - 1 = 5$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \end{bmatrix} $$

Solving the above system of 6 equations and 6 unknowns, we get

$$W_{1} = W_{3} = 0.557\;W_{2} = 0.889\; and \; \xi_{1} = 0.7745\;\;\xi_{2} = 0\;\;\xi_{3} = -0.7745$$

We have, $$a = -1\;\;;b = 1$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

$$\therefore\;\; \mathrm{dx} = \mathrm{d\xi}$$

The value of the function $$f(\xi) = \frac{\xi}{\xi^2+1}$$

$$I = \frac{l}{2}\int_{-1}^{1}(\frac{\xi}{\xi^2+1})\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(\frac{\xi_{1}}{\xi_{1}^2+1})\; + \;W_{2}(\frac{\xi_{2}}{\xi_{2}^2+1})\;+ W_{3}(\frac{\xi_{3}}{\xi_{3}^2+1})\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 0.556\left(\frac{0.7745}{0.7745^2+1}\right)+0.556\left(\frac{-0.7745}{0.7745^2+1}\right) = 0$$

Analytical Method

$$I =\int_{-1}^{1}(\frac{\xi}{\xi^2+1})\mathrm{d}\xi = 0,$$, using [Wolframalpha]

The analytical method also yields exactly the same result as the Gauss Quadrature.

(b). $$\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi}$$

(b).

$$I =\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi} $$

As $$n_{gp} = 3\; (three-point \;\;quadrature),\; p = 2n_{gp} - 1 = 5$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \end{bmatrix} $$

Solving the above system of 6 equations and 6 unknowns, we get

$$W_{1} = W_{3} = 0.557\;W_{2} = 0.889\; and \; \xi_{1} = 0.7745\;\;\xi_{2} = 0\;\;\xi_{3} = -0.7745$$

We have, $$a = -1\;\;;b = 1$$

$$\zeta = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

The value of the function $$f(\xi) = cos^2\pi\xi$$

$$I = \frac{l}{2}\int_{-1}^{1}cos^2\pi\xi\mathrm{d\xi}$$

$$I = \frac{l}{2}\left[W_{1}(cos^2\pi\xi_{1})\; + \;W_{2}(cos^2\pi\xi_{2})\;+W_{3}(cos^2\pi\xi_{3})\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$I = 0.556(cos^2\pi(0.7745))+0.889(cos^2\pi(0))+0.556(cos^2\pi(0.7745)) = 1.53$$

Analytical Method

$$I =\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi} = 1,$$, using [Wolframalpha]

The value from the analytical method is different from the value evaluated from Gauss Quadrature

Problem Statement
The integral $$\int_{-1}^{1}(3\xi^3+2)\;\mathrm{d\xi}$$ can be integrated exactly using two-point Gauss quadrature. How is the accuracy affected if a. one-point quadrature is employed; b. three-point quadrature is employed.

Solution
$$I =\int_{-1}^{1}(3\xi^3+2)\;\mathrm{d\xi}$$

As $$n_{gp} = 2\; (two-point \;\;quadrature),\; p = 2n_{gp} - 1 = 3$$

$$W_{1} = W_{2} = 1\;\; and \; \xi_{1} =0.5774\;\;\xi_{2} = -0.5774$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi_{1}^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(3\xi_{1}^3+2)\; + \;W_{2}(3\xi_{2}^3+2)\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 2+2 = 4$$

(a).One - Point Quadrature: $$W = 2\;\; and \; \xi =0$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W(3\xi^3+2)\right]\;and\;\; l = b - a = 1 - (-1) = 2$$

$$= 2(2) = 4$$

(b).Three - Point Quadrature:

$$W_{1} = W_{3} = 0.556\;\; and \; \xi_{1} =0.7745\;\;\xi_{2} = 0\;\;\xi_{3} =-0.7745$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi_{1}^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(3\xi_{1}^3+2)\; + \;W_{2}(3\xi_{2}^3+2)+\;\;W_{3}(3\xi_{3}^3+2)\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 2(0.556)(2)+(0.889)(2) = 4$$

'''The value of the Integral is same when evaluated with 0ne -,Two - and Three- point Quadratures. Hence all the three are accurate.'''

=Problem 6.8.2(a) Gauss Quadrature=

Problem Statement
Use Gauss quadrature to obtain exact values for the following integrals. Verify by analytical integration:

$$(a).\;\int_{0}^{4}(x^2+1)\mathrm{d}x$$

$$(b).\;\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

Solution
(a).

$$I =\int_{0}^{4}(x^2+1)\mathrm{d}x$$

As $$n_{gp} = 2\; (two-point \;\;quadrature),\; p = 2n_{gp} - 1 = 3$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \end{bmatrix} $$

Solving the above system of 4 equations and 4 unknowns, we get

$$W_{1} = W_{2} = 1\;\; and \; \xi_{1} = \frac{1}{\sqrt3}\;\;\xi_{2} = \frac{-1}{\sqrt3}$$

We have, $$a = 0\;\;;b = 4$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = 2+2\xi$$

The value of the function $$f(\xi) = (2+2\xi)^2+1$$

$$I = \frac{l}{2}\int_{-1}^{1}((2+2\xi)^2+1)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}((2+2\xi_{1})^2+1)\; + \;W_{2}((2+2\xi_{2})^2+1)\right] \;and \; l = b - a = 4 - 0 = 4$$

$$= 2\left[10.953+1.715\right] = 25.334$$

Analytical Method

$$I =\int_{0}^{4}(x^2+1)\mathrm{d}x$$

$$= \left[\frac{x^3}{3}+x\right]_{0}^{4} = \frac{4^3}{3}+4 = 25.333$$

The analytical method also yields exactly the same result as the Gauss Quadrature.

(b). $$I =\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

$$n_{gp} = 4\; (fo-point \;\;quadrature),\; p = 2n_{gp} - 1 = 7$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \xi_{1}^6&\xi_{2}^6 \\ \xi_{1}^7&\xi_{2}^7 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ W_{4}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \frac{2}{7}\\ 0\\ \end{bmatrix} $$

Solving the above system of 8 equations and 8 unknowns, we get

$$W_{1} = W_{3} = 0.348\;W_{2} = W_{4} = 0.652\; and \; \xi_{1} = 0.861\;\;\xi_{2} = 0.339\;\;\xi_{3} = -0.861\;\;\xi_{4} = -0.339$$

We have, $$a = -1\;\;;b = 1$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

The value of the function $$f(\xi) = (\xi^4+2\xi^2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(\xi^4+2\xi^2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(\xi_{1}^4+2\xi_{1}^2)\; + \;W_{2}(\xi_{2}^4+2\xi_{2}^2)\;+ (\xi_{3}^4+2\xi_{3}^2)\;+ (\xi_{4}^4+2\xi_{4}^2)\;\right] \;and \; l = b - a = 4 - 0 = 4$$

$$= 1.4144+0.3169 = 1.731$$

Analytical Method

$$I =\int_{-1}^{1}(\xi^4+2\xi^2)\mathrm{d}\xi$$

$$= \left[\frac{\xi^5}{5}+\frac{2\xi^3}{3}\right]_{-1}^{1} = 2\left[\frac{1}{5}+\frac{2}{3}\right] = 1.731$$

The analytical method also yields exactly the same result as the Gauss Quadrature.

=Problem 6.8.2(b)Three-point Gauss quadrature=

Problem Statement
Use three-point Gauss quadrature to evaluate the following integrals. Compare to the analytical integral.

(a).$$ \int_{-1}^{1}\frac{\xi}{\xi^2+1}\;\mathrm{d}x $$

(b). $$\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi}$$

Solution
(a).

$$I =\int_{-1}^{1}\frac{\xi}{\xi^2+1}\;\mathrm{d}x $$

As $$n_{gp} = 3\; (three-point \;\;quadrature),\; p = 2n_{gp} - 1 = 5$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \end{bmatrix} $$

Solving the above system of 6 equations and 6 unknowns, we get

$$W_{1} = W_{3} = 0.557\;W_{2} = 0.889\; and \; \xi_{1} = 0.7745\;\;\xi_{2} = 0\;\;\xi_{3} = -0.7745$$

We have, $$a = -1\;\;;b = 1$$

$$x = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

$$\therefore\;\; \mathrm{dx} = \mathrm{d\xi}$$

The value of the function $$f(\xi) = \frac{\xi}{\xi^2+1}$$

$$I = \frac{l}{2}\int_{-1}^{1}(\frac{\xi}{\xi^2+1})\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(\frac{\xi_{1}}{\xi_{1}^2+1})\; + \;W_{2}(\frac{\xi_{2}}{\xi_{2}^2+1})\;+ W_{3}(\frac{\xi_{3}}{\xi_{3}^2+1})\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 0.556\left(\frac{0.7745}{0.7745^2+1}\right)+0.556\left(\frac{-0.7745}{0.7745^2+1}\right) = 0$$

Analytical Method

$$I =\int_{-1}^{1}(\frac{\xi}{\xi^2+1})\mathrm{d}\xi = 0,$$, using [Wolframalpha]

The analytical method also yields exactly the same result as the Gauss Quadrature.

(b). $$\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi}$$

(b).

$$I =\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi} $$

As $$n_{gp} = 3\; (three-point \;\;quadrature),\; p = 2n_{gp} - 1 = 5$$

$$ \begin{bmatrix} 1 &1 \\ \xi_{1}&\xi_{2} \\ \xi_{1}^2&\xi_{2}^2 \\ \xi_{1}^3&\xi_{2}^3 \\ \xi_{1}^4&\xi_{2}^4 \\ \xi_{1}^5&\xi_{2}^5 \\ \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2}\\ W_{3}\\ \end{bmatrix} = \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0\\ \frac{2}{5}\\ 0\\ \end{bmatrix} $$

Solving the above system of 6 equations and 6 unknowns, we get

$$W_{1} = W_{3} = 0.557\;W_{2} = 0.889\; and \; \xi_{1} = 0.7745\;\;\xi_{2} = 0\;\;\xi_{3} = -0.7745$$

We have, $$a = -1\;\;;b = 1$$

$$\zeta = \frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) = \xi$$

The value of the function $$f(\xi) = cos^2\pi\xi$$

$$I = \frac{l}{2}\int_{-1}^{1}cos^2\pi\xi\mathrm{d\xi}$$

$$I = \frac{l}{2}\left[W_{1}(cos^2\pi\xi_{1})\; + \;W_{2}(cos^2\pi\xi_{2})\;+W_{3}(cos^2\pi\xi_{3})\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$I = 0.556(cos^2\pi(0.7745))+0.889(cos^2\pi(0))+0.556(cos^2\pi(0.7745)) = 1.53$$

Analytical Method

$$I =\int_{-1}^{1}cos^2\pi\zeta\mathrm{d\xi} = 1,$$, using [Wolframalpha]

The value from the analytical method is different from the value evaluated from Gauss Quadrature

=Problem 6.8.2(c) Check for accuracy in Gauss Quadrature=

Problem Statement
The integral $$\int_{-1}^{1}(3\xi^3+2)\;\mathrm{d\xi}$$ can be integrated exactly using two-point Gauss quadrature. How is the accuracy affected if a. one-point quadrature is employed; b. three-point quadrature is employed.

Solution
$$I =\int_{-1}^{1}(3\xi^3+2)\;\mathrm{d\xi}$$

As $$n_{gp} = 2\; (two-point \;\;quadrature),\; p = 2n_{gp} - 1 = 3$$

$$W_{1} = W_{2} = 1\;\; and \; \xi_{1} =0.5774\;\;\xi_{2} = -0.5774$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi_{1}^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(3\xi_{1}^3+2)\; + \;W_{2}(3\xi_{2}^3+2)\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 2+2 = 4$$

(a).One - Point Quadrature: $$W = 2\;\; and \; \xi =0$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W(3\xi^3+2)\right]\;and\;\; l = b - a = 1 - (-1) = 2$$

$$= 2(2) = 4$$

(b).Three - Point Quadrature:

$$W_{1} = W_{3} = 0.556\;\; and \; \xi_{1} =0.7745\;\;\xi_{2} = 0\;\;\xi_{3} =-0.7745$$

$$f(\xi) = (3\xi^3+2)$$

$$I = \frac{l}{2}\int_{-1}^{1}(3\xi_{1}^3+2)\; \mathrm{d}\xi$$

$$I = \frac{l}{2}\left[W_{1}(3\xi_{1}^3+2)\; + \;W_{2}(3\xi_{2}^3+2)+\;\;W_{3}(3\xi_{3}^3+2)\right] \;and \; l = b - a = 1 - (-1) = 2$$

$$= 2(0.556)(2)+(0.889)(2) = 4$$

'''The value of the Integral is same when evaluated with 0ne -,Two - and Three- point Quadratures. Hence all the three are accurate.'''

=Problem 6.8 Verify the given table=

Problem Statement
Verify the table of {(wi,xi),i=1,2,..,5} against NIST handbook and also with table on textbook page 89

Solution
The table from Lecture 36-3 is given below:

Nodes and weights for the 5-point Gauss–Legendre formula from | NIST Handbook are given below:

The values of the Nodes and Weights in the given table in Lecture 36-3 are calculated and the values are tabulated:

The table of Nodes and Weights with the calculated values from the text book Introduction to Finite Element Analysis by Fish and Belytschko is as below:

From the above tables, it can be concluded that the values of Nodes and Weights from the table in Lecture 36-3, NIST handbiik and the text book agree with each other.

=Problem 6.9 Reproduce Exercise8.3, 8.4 and Problem 8.6 in Textbook = Refer to P200, P205 and P212 of Fish & Belytschko's A First Course in Finite Elements, 2007, John Wiely & Sons, Ltd, for more detailed description.

Problem Statement
The conductivity is isotropic and K=5. The temperature T=0 is prescribed along bottom edge and left edge.

The heat flux q=0 and q=20 are prescribed on the right edge and top edge, respectively.

A constant heat source s=6 is applied over the plate.

Analysis Results
Using provided Matlab codes, we can obtain following results for using different number of element.

Instruction on using the provided Matlab codes:
heat2d.m is the main program;

When running heat2d.m, it will call functions of include_flags, [K,f,d] = preprocessor, [ke, fe] = heat2Delem(e), [K,f] = assembly(K,f,e,ke,fe), src_and_flux(f), [d,f_E] = solvedr(K,f,d) and postprocessor(d) step by step.

In the include_flags function, we define global variables which may be used in following functions and subprograms.

In the [K,f,d] = preprocessor function, we define how many elements we gonna use in this program, and there are three types provided, i.e. 1element, 16elements and 64elements. While doing analysis, we can change element number here. In each element type input file, we define material property, essential BCs, natural BCs, mesh generation and also plot the meshed figure.

In the [ke, fe] = heat2Delem(e) function, we compute the elemental conductance matrix and nodal flux vector. In this function, the Gauss function for integration is also called.

In the [K,f] = assembly(K,f,e,ke,fe) function, we assemble the elemental conductance matrix and the heat source vector.

In the src_and_flux(f) function, we Compute and assemble nodal boundary flux vector and point sources.

In the [d,f_E] = solvedr(K,f,d) function, we partition and solve the system of equations.

In the postprocessor(d) function, we plot temperature and flux.

Matlab Codes

From below website, we can download all the required Matlab Codes. All Rights Reserved.

The specific codes for this case are in the package of Chapter 8.

Matlab Codes donwload

Problem Statement
Since the exercise haven’t specify the parameters a and b in the textbook, from the plot provided, we take a=0.25 and b=1 for following analysis. For heat equation with isotropic conductivity and K=1, the corresponding source term is given by
 * {| style="width:100%" border="0"

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

s=-{{\nabla }^{2}}T=2a\frac{1}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-4

$$ The essential boundary conditions on $${{\Gamma }_{T}}$$ are
 * }.
 * {| style="width:100%" border="0"

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

T(r=a)=0,T(x=\pm b,y)={{a}^{2}}+{{b}^{2}}+{{y}^{2}}-2a\sqrt{{{y}^{2}}+{{b}^{2}}}

$$ The natural boundary conditions on $${{\Gamma }_{q}}$$ are ($$\bar{q}=-k{{n}^{T}}\nabla T$$)
 * }.
 * {| style="width:100%" border="0"

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

\begin{align} & \bar{q}(x,y=b)=-\frac{\partial T}{\partial y}(x,y=b)=2b(\frac{2a}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-1) \\ & \bar{q}(x,y=-b)=-\frac{\partial T}{\partial y}(x,y=-b)=2b(1-\frac{2a}{\sqrt{{{x}^{2}}+{{y}^{2}}}}) \\ \end{align}

$$
 * }

Analysis Results
Based on provided Matlab codes for exercise 8.3, page 200, we modified the code to match the problem 8.4, page 205.

Matlab Code Package for Example Problem 8.4

Below is the analysis results:

Problem Statement
Consider a chimney constructed of two isotropic materials: dense concrete ($$k=2.0W{}^{\circ }{{C}^{-1}}$$) and bricks ($$k=0.9W{}^{\circ }{{C}^{-1}}$$).

The temperature of the hot gases on the inside surface of the chimney is $$140{}^{\circ }C$$, whereas the outside is exposed to the surrounding air, which is at $$T=10{}^{\circ }C$$. The dimensions of the chimney (in meters) are shown below.

For the analysis, exploit the symmetry and consider 1/8 of the chimney cross sectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials.

Analyze the problem with 2×2, 4×4 and 8×8 quadrilateral elements for 1/8 of the problem domain. A2×2 finite element mesh is shown in Figure below. Symmetry implies insulated boundary conditions on edges AD and BC. Note that element boundaries have to coincide with the interface between the concrete and bricks.

Analysis Results
Based on provided Matlab codes for exercise 8.3, page 200, we modified the code to match the problem 8.6, page 212.

Matlab Code Package for Problem 8.6

Below is the analysis results:

=Problem 6.10=

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
 * 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 the semi-discrete form which is as:-

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

$$ (Eq )<12>
 *
 * }

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

$$ (Eq )<13>
 *
 * }

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 }

$$ (Eq )<14>
 * <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 )<15>
 * 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 plot obtained from www.wolframaplha.com is shown as below. Contour Plot

=Contributing Members & Referenced Lecture=

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

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


text:
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" |
 *

 * }

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

 * }

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

 * }

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.

=Contributing Members & Referenced Lecture=

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