User:Eml5526.s11.team3.akj/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.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" |

{{K}_{FF}}= \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" |

{{F}_{F\times 1}}=\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" |

{{d}_{F\times 1}}=\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