User:Eml5526.s11.team01/Hwk6

 HOMEWORK 6

=Problem 6.1: Quadratic Lagrangian Element Basis Function (QLEBF)= Solving a ordinary differential equation using Quadratic Lagrangian Element Basis Function (QLEBF). From Mtg.31

Given
The partial differential equation is:

The boundary conditions are:

$$ u\left(0\right) = 4$$

$$\frac{du}{dx}\left(1\right) = \frac{12}{5} $$

Find

 * 1. For nel = 2, compute $$\tilde{K}=\sum_{1}^{2}\tilde{K}^e$$ with $$\tilde{K}^e$$ given by equation (1) p.30-5
 * 2. Compute $$\tilde{k}^e$$ and $$L^e$$ for e = 1,2


 * 3. Compute $$\tilde{K}^e=L^{eT}k^eL^e$$ for e = 1,2. Compare to 1)
 * 4. Plot all the QLEBF for nel = 3
 * 5. Plot the exact solution, the numerical solution, and the error as a function of the number of elements

For nel = 2, compute $$\tilde{K}=\sum_{1}^{2}\tilde{K}^e$$ with $$\tilde{K}^e$$ given by equation (1) p.30-5
To calculate the stiffness matrix of the element we need the weak form; but, the weak form was already obtained in home work 5 which can be reach at http://en.wikiversity.org/wiki/User:Eml5526.s11.team01/Hwk5 Then, here we only copy that equation.

Now we use QLEBF to interpolate the solution $$U^h$$ and the weight funtion w. The aproximated soution can be interpolated as:


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

$$U^h=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}u_1+\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}u_2+\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}u_3 $$
 * $$\displaystyle (6.1.2) $$
 * }
 * }

Where $$u_{1}, u_{2}, u_{3}$$ are the nodal values of the funtion at the nodes 1, 2, and 3 respectively. The last equation can be written as:


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

$$U^h=\left [ N \right ]\left [ u \right ]^{e} $$
 * $$\displaystyle (6.1.3) $$
 * }
 * }

where $$\left [ u \right ]^{e}=\begin{bmatrix} u_1\\ u_2\\

u_3\end{bmatrix} $$

In the same way we can interpolate the weight funtion w to obtain:


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

$$w=\left [ N \right ]\left [ w \right ]^{e} $$
 * $$\displaystyle (6.1.4) $$
 * }
 * }

where $$\left [ w \right ]^{e}$$ is the vector with the nodal values of w

The derivative of $$U^h$$ and w are:


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

$$\frac{dU^h}{dx}=\begin{bmatrix} \frac{2x-x_2-x_3}{\left (x_1-x_2 \right )\left (  x_1-x_3\right )} & \frac{2x-x_1-x_3}{\left (x_2-x_1  \right )\left (  x_2-x_3\right )}  & \frac{2x-x_1-x_2}{\left (x_3-x_1  \right )\left (  x_3-x_2\right )} \end{bmatrix}\left [ u \right ]^e $$ The latter equation can be written as:
 * $$\displaystyle (6.1.5) $$
 * }
 * }


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

$$\frac{dU^{h}}{dx}=\left[B\right]\left [ u \right ]^{e} $$
 * $$\displaystyle (6.1.6) $$
 * }
 * }

and in similar way w can be written as:


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

$$\frac{dw}{dx}=\left[B\right]\left [ w \right ]^{e} =\left [ B \right ]^T\left [ w \right ]^{eT}$$ Replacing equations (6.1.3), (6.1.4), (6.1.6), and (6.1.7) in the weak form equation (5.7.4) we obtaine:
 * $$\displaystyle (6.1.7) $$
 * }
 * }


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

$$\left [ w \right ]^{eT}\int_{0}^{l}\left [ B \right ]^{T}\left ( 2+3x \right )\left [ B \right ]dx\left [ u \right ]^{e}=5\left [ w\right ]^{eT}\int_{0}^{l}\left [ B \right ]^{T}x\times dx-30w^{1}$$
 * $$\displaystyle (6.1.8) $$
 * }
 * }

The last equation is the weak form for any element; that is why the integration was done from 0 to l where l is the lenght of the element. Also in the last equation the factor $$w^{1}$$ in the nodal value of w at x = 1; therefore that factor is present in the weak form for the last element. Cancelling the factor $$\left [ w \right ]^{eT}$$ in equation (6.1.8)we obtain an equation of the form:$$Kd=f$$; therefore, we can extract from that equation the stiffness matrix which is:


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

$$k^{e}=\int_{0}^{l}\left [ B \right ]^{T}\left ( 2+3x \right )\left [ B \right ]dx$$
 * $$\displaystyle (6.1.9) $$
 * }
 * }

We developed a Matlab code (shown below) to calculate the global stiffness matrix $$\tilde{K}$$ using equation (1) p.30-5. The obtained matrix was:

Compute $$\tilde{k}^e$$ and $$L^e$$ for e = 1,2
Using the weak form developed above we calculate the stiffness matrix of the elements 1 and 2. As each elemet has 3 degrees of freedom and the total number of degrees of freedom is 7, the Location matrix is 3X5. The stiffness matrix are:


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

$$k^{1}=\begin{bmatrix} 10.83 & -12.66 &1.83 \\ -12.66 & 29.33 & -16.66\\ 1.83 & -16.66 & 14.83 \end{bmatrix} $$
 * $$\displaystyle (6.1.11) $$
 * }
 * }


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

$$k^{2}=\begin{bmatrix} 17.83 & -20.66 &2.83 \\ -20.66 & 45.33 & -24.66\\ 2.83 & -24.66 & 21.83 \end{bmatrix} $$ And the location matrix are:
 * $$\displaystyle (6.1.12) $$
 * }
 * }


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

$$L^{2}=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 1 & 0 &0 & 0 & 0 & \\ 0 & 0 & 1 & 0 &0  & 0 & \end{bmatrix} $$
 * $$\displaystyle (6.1.13) $$
 * }
 * }


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

$$L^{2}=\begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & \\ 0 & 0 & 0 &0 & 1 & 0 & \\ 0 & 0 & 0 & 0 &0  & 1 & \end{bmatrix} $$
 * $$\displaystyle (6.1.14) $$
 * }
 * }

Now the global matrix are calculated by using equation (2) p.31-1 which is repeated here:


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

$$\tilde{k}^{e}=L^{eT}k^{e}L^{e} $$
 * $$\displaystyle (6.1.15) $$
 * }
 * }

Therefore, for the element 1 whe have:


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

$$\tilde{k}^{1}=\begin{bmatrix} 1 & 0 &0 \\ 0 & 1 &0 \\ 0 & 0 &1 \\ 0 &0 &0 \\ 0 & 0 &0 \\ 0 & 0 & 0\\  0& 0 & 0 \end{bmatrix}\begin{bmatrix} 10.83 & -12.66 & 1.83\\ -12.66 & 29.33 &-16.66 \\ 1.83 & -16.66 & 14.83 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0& 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{bmatrix}=\begin{bmatrix} 10.83 & -12.66 & 1.83 & 0 & 0\\ -12.66 & 29.33 & -16.66 & 0 &0 \\ 1.83& -16.66 & 14.84 & 0 & 0\\  0& 0 & 0 & 0 & 0\\  0& 0 & 0 &0  & 0 \end{bmatrix} $$
 * $$\displaystyle (6.1.16) $$
 * }
 * }

In similar way for the elment 2 we have:


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

$$\tilde{k}^{2}=\begin{bmatrix} 0 & 0 &0 \\ 0 & 0 &0 \\ 0 & 0 &0 \\ 0 &0 &0 \\ 1 & 0 &0 \\ 0 & 1 & 0\\  0& 0 & 1 \end{bmatrix}\begin{bmatrix} 17.83 & -20.66 & 2.83\\ -20.66 & 45.33 &-24.66 \\ 2.83 & -24.66 & 21.83 \end{bmatrix}\begin{bmatrix} 0 & 0 & 0& 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 &0 \\ 0& 0& 17.83 & -20.66 & 2.83\\  0& 0 & -20.66 & 45.33 & -24.66\\  0& 0 &2.830 &-24.66  & 21.83 \end{bmatrix}$$
 * $$\displaystyle (6.1.17) $$
 * }
 * }

Performing the sum we obtain:

We can observ that this equation is equal to equation (6.1.10) because it corresponds to the same system with the same weak form. The location matrix is only a methodology for asembly the global matrix and therefore it can not change the final result.

Plot all the QLEBF for nel = 3
Qubic Lagrange Element Basis Function for every element are given by:


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

$$L_{1}=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}; L_{2} =\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}; L_{3} = \frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}u_3 $$
 * <p style="text-align:right;">$$\displaystyle (6.1.19) $$
 * }
 * }

Where $$x_{1}, x_{2}, and x_{3}$$, are the coordinates of the nodes 1, 2, and 3 respectively. This funtions are repeated on all the elements. Figure 6.1.1 shows the QLEBF for 3 elements, where it can be observed that every basis funtion has the value of 1 at the respective node and 0 on the others, of course the basis function are valid only inside their respective element. In the Figure, Li,j corresponds to j the basis funtion of the element i



Plot the exact solution, the numerical solution, and the error as a function of the number of elements
The analitical solution was developed before in other homeworks; therefore,here we only copy it

Figures 6.1.2 shows the analitical and numerical solution obtained with different number of elements. Although good aproximation was achieved even whit low number of elements, the convergence with a tolerance of $$10^{-6}$$ was not achieved even with more than 100 elements. With a tolerance of $$10^{-4}$$ convergence was achieved with 54 elements. Figure 6.1.3 shows the error as a function of the number of elements. In that figure, it can be observed that initially the error decreases fast but from aproximately 25 elements the relation between the number of elements and the error in logaritmic sacale tend to be linear and to achieve a tolerance of $$10^{-6}$$ more than 200 elements could be needed; but, is probalby that with higher number of elements the error decrease even slower because with more element the round-off error increases.





Appendix
 Matlab Code: 

=Problem 2:Weak Form Using Quadratic Lagrangian elemental Basis Function With Uniform Discretization.=

Given
General 1-D Model 1 Data set 1. The Partial Differential equation for this data set is,

Find
1)For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$

2)Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$

3)Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1

4)Plot all QLEBF for $$n_{el}=3$$.

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

Solution
Lets Obtain the Weak form first of all as follow

Performing Integration by parts technique we get

Solving this we get,

As we do not know $$ \frac{du}{dx}\left ( 1\right )$$ we choose $$w\left ( 1 \right )=0$$ and moreover $$u^{h}\left(1 \right )=4$$. To use quadratic Lagrangian elemental basis function we divide the dominium in elements, on each element we apply the weak form by using quadratic Lagrangian Polynomial. For that we need 3 nodes per element. Each element with global nodes i,i+1,i+2 will have 3 local nodes 1,2,3. Lets call $$x_{1},x_{2},x_{3}$$the co-ordinates at node 1,2,3 respectively.

To use natural co-ordinate $$\epsilon$$ instead of Co-ordinate $$x$$ we use the equation 4 on page 30-3 for that we know the natural co-ordinates of the nodes are -1,0,1 & local nodes 1,2,3 respectively. Using Lagrange Polynomial to map between two Co-ordinate system we have:

substituting the above values in the $$\displaystyle (6.2.7)$$ we get,

In this way we have found a shape function, and the later equation can be written as follow,

In the same way we can interpolate the function $$u$$ & $$w$$, therefore,

The Derivative of function $$u$$ with respect to $$\epsilon$$ is given by, In the same way, If we place the Node 2 at the midpoint of the element we can write the relation between the two co-ordinates as Taking the derivative on both sides we get, The derivative of $$u$$ & $$w$$ with respect to $$x$$ are calculated by using chain rule which is as below, Replacing all the terms in weak form as a function of $$\epsilon$$ we obtain This can be further simplified as

From the later Equation we can get the Stiffness Matrix of one element as follows,

and the force vector of the element will be,

From the later Equation we can get the Stiffness Matrix of one element as follows,

There will be additional following factor,

From the later Equation we can get the Stiffness Matrix of one element as follows,

Now only count for degree of freedom at x=1.

1) For $$\displaystyle n_{el}=2$$ Compute $$\underline{\tilde{K}}$$ Using Equation 1 On page 30-5. Display $$\underline{\tilde{K}},e=1,2$$.
Now according to the Lagrange Interpolation we have following equation,

Also,

Replacing in the weak form above we get, The last term is taken only for the last element. Further the stiffness matrix is

Performing the integration for the code shown below for one element we obtain, Hence the global stiffness matrix will be,

2) Compute $$K^{e}$$, $$L^{e}$$ for $$e=1,2$$.
$$K^{e}$$ is already calculated above. For $$L^{e}$$ As the element has 3 degrees of freedom & the total number of degrees of freedom is 5 The L matrix will be of size $$3\times5$$. Therefore,

3) Compute $$\underline{\tilde{K}}^{e}=\underline{\tilde{L}}^{e}\underline{K}^{e}\underline{L}^{e}, for e=1,2$$ & Compare to 1.
Now, Solving the above matrix we get,

For element 2 we have,

Solving the above matrix we get,

Hence Combining equations $$\displaystyle (6.2.17)$$ & $$\displaystyle (6.2.18)$$ we get global stiffness matrix as follows,

Eventually Comparing $$\displaystyle (6.2.15)$$ with $$\displaystyle (6.2.19)$$ we found that the Stiffness matrices are identical.

4) Plot all QLEBF for $$n_{el}=3$$.
The Quadratic Lagrange Elemental basis Function is plotted using the following MATLAB code for $$\displaystyle n_{el}=3$$.

Appendix
 Matlab Code: 

The plot is as follow.



5) Plot $$\displaystyle \tilde{u_{n}}^{h}$$ vs $$\displaystyle u$$ and [$$\tilde{u_{n}}^{h}\left ( 0.5\right )-u\left (0.5\right )$$] vs $$\tilde{n}$$.
These Two plots are plotted as below Using the mentioned MATLAB code.

Appendix
 Matlab Code: 





===Comment: From the above plot we can observe that the solution converges at nel=16. The Convergence of $$10^{-6}$$ was difficult to achieve therefore the convergence criteria was modified to $$E=\left(u_{0.5}^{i}- u_{0.5}^{i-1}\right )$$.===

=Problem 3: Calculix=

Given
Given the tutorial found in this tutorial.

Find
Reproduce the actions of the author in Calculix

Solution
The tutorial was followed and reproduced sucessfully. The actions and results were summerized in the following

=Problem 4: FB pp.148-149 Problems 6.1,6.2,6.3=

Given
6.1)

Given vector $$ q_x = -y^2 ~,~ q_y = -2xy $$ on domain shown in Figure 6.2.

6.2)

Given a vector field $$ q_x = 3x^2 y +y^3 ~,~ q_y = 3x + y^3 $$ on the domain shown in Figure 6.9. The curved boundary of the domain is a parabola.

6.3) Given $$ \oint_{T} \mathbf{n} d\Gamma = 0 $$

Find
6.1)

Verify the divergence theorem.

6.2)

Verify the divergence theorem.

6.3)

Prove using divergence theorem

FB 6.1)


We have to show that:

Where $$ \Omega $$ is the domain, $$ \Gamma $$ is the boundary, n is the normal vector and the notation $$ q_{i,i} $$  means $$ \frac{ \partial q_i}{ \partial x_i} $$

The integration on the boundary is

FB 6.2)
To calculate the integral on the surface we see that part of the surface is the parabola and the other part is in the x axis; therefore the integral will be:

The first integral is on the parabola and the second part is on the x axis.

On the curved surface the, normal vector on a segment of surface ds is shown in the figure. From that figure we observe that:

Also, it is clear that on the x axis we have y=0. Replacing in the integral we obtain:

Taking the derivative on the equation of the parabola we obtain:

Replacing on the integral we have:

Changing the limits of integration in the first two terms and organizing we obtain:



The integral on the volume will be:

Performing first the integration with respect to y we have:

FB 6.3)
The theorem says that for any tensor $$ T_{i...j} $$ (of any order) inside a volume $$ \Omega $$ surrounded by a surface $$ \Gamma $$ with a normal vector n

In this case the tensor T is the constant 1; therefore, changing the integral of surface to an integral of volume we have

But as 1 is a constant, its derivatives with respect to any coordinate are zero and we have

=Problem 6.5: Solve G2DM10/D1 using 2DLIBF=

Given: 2DLIBF from p. 29-2
PDE

$$ K = I $$ (identity matrix)

$$ f = 0, \frac{\partial u}{\partial t} = 0 $$

$$ \frac{\partial}{\partial x_i} \left(k_{ij} \frac{\partial u}{\partial x_j} \right) = 0 $$

Essential B.C.

$$ g = 2 on  \partial \Omega $$

Natural B.C.

none

Find
Solve G2DM10/D1 for $$ u^h $$ until accuracy $$ 10^{-6} $$ at center $$ (x,y) = (0,0) $$

Solution
Using the weak form to be

and

2DLIBF

where

$$ I = i + \left( j-1 \right) m $$

$$ n = m = 2, 4, 6, 8... $$

Thus

Therefore by substituting equation 6.5.3 into equation 6.5.1 and integrating, the stiffness matrix can be solved for. This is done by using Matlab.

The solution yields



=Problem 6=

Solution
=Problem 7=

Solution
= Problem 6.8 =

Verify the Table of {(wi, xi), i=1,2,3,4,5} against the NIST Handbook and that in the textbook
The Table of $$\left\{ \left( {{w}_{i}},{{x}_{i}} \right),\ i=1,2,3,4,5 \right\}$$ is shown as follows:

Complare the positions of points and the weight with respect to each of them in the Table above with NIST Handbook and Fish and Belytschko's book on page 89.

Solution
Calculation of the value of each term in the above table gives the following result:

The nodes and weights for the 5-point Gauss–Legendre formula in the NIST Handbook is the following:

The position of Gauss points and corresponding weights in Fish and Belytschko's book on page 89 can be seen as follows:

Conclusion:

It can be seen from the above several tables that the table of $$\left\{ \left( {{w}_{i}},{{x}_{i}} \right),\ i=1,2,3,4,5 \right\}$$ given in Lecture 36 gives us quite precise results of the position of Gauss points and corresponding weights compared to the other two authoritative sources.

Use Gauss quadrature to obtain exact values for the given integrals and verify by analytical integration.
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
(a)

Select the number of Gauss point, n:

$$2n-1\ge 2$$, select n=2.

Then,

By analytical integration:

So, the result performed by analytical integration is identical to that obtained by Gauss quadrature.

(b)

Select the number of Gauss point, n:

$$2n-1\ge 4$$, select n=3.

By analytical integration:

So, the result performed by analytical integration is identical to that obtained by Gauss quadrature.

(c)

Comments:

The results for (a) and (b) using MATLAB code are 25.3333 and 1.7333  respectively, which are the same with those calculated by analytical integration.

Use three-point Gauss quadrature to evaluate the given integrals and compare to the analytical integration.
(c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB code.

Solution
(a)

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Integrate using WolframAlpha, we get the same result.

(b)

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Integrate using WolframAlpha, we get the result of I=1, which is different from what we get using three-point Gauss quadrature.

(c)

Comments:

The result for (a) using MATLAB code is 0, which is the same with that calculated by analytical integration.

And the result for (b) using MATLAB code is 1.53, which is different from that calculated by analytical integration.

Employ one-point and three-point Gauss quadrature and evaluate the accuracy of the given integral compared to the result by two-point quadrature.
The integral $$\int\limits_{-1}^{1}{\left( 3{{\xi }^{3}}+2 \right)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.

Check your calculations against MATLAB code.

Solution
a.

n=1, then we have the following according to the Table of {(wi, xi)}:

So,

b.

n=3, then we have the following according to the Table of {(wi, xi)}:

So,

Comments:

The results for (a) and (b) using MATLAB code are all equal to 4 , which are the same with those calculated by analytical integration.

Conclusion:

It can be seen from the above Gauss quadrature calculation and the verification using MATLAB code that the accuracy is not affected when the number of Gauss points is either increased or decreased in this problem.

Matlab Code
=Problem 9=

Solve the 2D heat conduction problem with 16 elements using the finite element code.
Consider the heat conduction problem depicted in Figure 6.9.1. The coordinates are given in meters. The conductivity is isotropic, with $$\mathbf{D}=k\left[ \begin{matrix} 1 & 0 \\   0 & 1  \\ \end{matrix} \right]$$, and $$k=5\text{W}{}^{\text{o}}{{\text{C}}^{\text{-1}}}$$. The temperature $$T=0$$ is prescribed along edges AB and AD. The heat fluxes $$\bar{q}=0$$ and $$\bar{q}=20\text{W}{{\text{m}}^{\text{-1}}}$$ are prescribed on edges BC and CD, respectively. A constant heat source $$s=6\text{W}{{\text{m}}^{\text{-2}}}$$ is applied over the plate. Consider this problem modeled with 16 quadrilateral finite elements as shown in Figure 6.9.2.

Solve this problem using the finite element code given in Section 12.5.



Solution


The MATLAB code are downloaded at FB Student Companion Site.

Show the temperature distribution for the coarsest (34-element) and the finest (502-elements) meshes.
In this example, we consider a manufactured solution of the form

defined over the domain of a square plate with a hole. For heat equation with isotropic conductivity and k=1, the corresponding source term is given by

The essential boundary conditions on $${{\Gamma }_{T}}$$ are

The natural boundary conditions on $${{\Gamma }_{q}}$$ are $$\left( \overline{q}=-k{{\mathbf{n}}^{T}}\nabla T \right)$$

Solution
Due to the hole in the midle of the domain, we find dificult to generate the mesh by hand, therefore we use a comercial program to solve the problem. First we start with a coarse mesh withwith 244 elements (Figure 6.9.5). The results are shown in Figure 6.5.6.



With 244 elements the results are show in the following figure.



Figure 6.9.7 shows the results across the plate at y=0.5.



Then the mesh was refine to 1116 elements (Figure 6.9.8), and again we show the results at y = 0.5 (Figure 6.9.9). When this results were compared with the analytical solution (Figure 6.9.11) an important difference was found. This diffence may be due to some error in our procedure, but also can be due to the size of the hole. As in the book the diemeter of the hole is not given we chose a diameter of 0.5 which apparently is lower than that used in the book.



The following figure shows the results with 1116 elements.



The analitical solution is shown below.



Determine the temperature and flux in the two materials of a chimney with 2x2, 4x4 and 8x8 quadrilateral elements for 1/8 of the problem domain.
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.



Solution
Figure 6.9.13 shows the initial mesh used to solve the problem. With this mesh, we obtained the temperature shown in Figure 6.9.14 and the heat flux shown in Figure 6.9.15. Then, we refine the mesh as shown in Figure 6.9.16 and obtained the temperature and heat flux shown in Figure 6.9.16 and 6.9.17 respectively. In the heat flux results we observed a high flux in the right border in the contact between the two materials. Finally, an extrafine mesh was used (Figure 6.9.19) to obtain the temperature shown in Figure 6.9.20 and the heat flux shown in Figure 6.9.21. It can be observed that with extrafine mesh the maximum values of heat flux obtained are higher than those obtained with coarse mesh and the regions of maximum heat flux are narrow. This behaviour was not observed on the results of temperature where the results are similar with the three mesh.



















=Problem 10: Newtons Law of Colling=

Given
The partial Differential Equation is

1.Develop Weak form similar to (3)-(6) page 34-3.
The weak form will be as below,

Lets calculate the derivative of the following factor, Therefore, Replacing in the integral we obtain, Applying the Gauss Theorem on the first term we obtain, We can recognize that the Heat Flux is Given by, Eventually We have,

Organizing We get, Simplifying we get, For the Case when Conductivity if the Identity Matrix, we have following equations, Also We can Write as,

2) Develop Semidiscrete Equations (ODE's) similar to (3) page 23-3. Give Detailed Description for all the quantities.
Performing the Discretization on the Dominium we have,

Where N is the vector of the Interpolation (Lagrange Polynomial) and $$\displaystyle u_{e},w_{e}$$ are the Nodal Values of u & w Respectively. The Derivative With respect to $$\displaystyle x_{1}$$, $$\displaystyle x_{2}$$ are Replacing on the Semi Discrete Weak form we obtain, From the Last Equation we can Obtain the Conductivity Matrix, the Heat source vector and the Capacitance matrix. The Conductivity Matrix is Given by

Further the Capacitance Matrix is Given by,

And the Heat Source is

3.Solve G2DM2.0/D1 Using 2D LIBF (similar to HW 6.6) till $$10^{-6}$$ accuracy at center. Plot Solution in 3D w/Contour.
To solve G2DM2.0/D1 Following Routine was developed. The following figure shows the result for 16 points.

To solve the problem using LIBF we will devide the x-axis and Y-axis in 2 parts initially. In that way the Approximated solution by LIBF is given by:

We Have developed a MATLAB code to plot this routine. The MATLAB code is as follows

Appendix
 Matlab Code: 

The Plot for the Routine Developed along with the intended 3-D contour is as below.



Comment: The Contour is plotted for $$\displaystyle n_{el}=4$$. The above MATLAB code is developed for larger values of $$n_{el}$$.
=Contributing Members=

= References =