User:Eml5526.s11.team6.kurth/HW6

=Problem 6.2: Solving Data Set G1DM1.0/D1 Using Quadratic Lagrangian Element Basis Functions and Uniform Discretization=

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [1]$$, is

The Natural Boundary condition as applied on the boundary $$\Gamma_h = [0]$$, is

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 2, 4, 6, 8, ..$$

Find
1.) For $$n_{el} = 2$$, compute $$\tilde{\underline{K}}=\sum_{e=1}^2\tilde{\underline{K}}^e$$, with $$\tilde{\underline{K}}^e$$ by [[media: Fe1.s11.mtg30.djvu| (1)p.30-5]], display ,$$\tilde{\underline{K}}^e,\;e=1,2.$$

2.) Comp. $$\underline{k}^e,\underline{L}^e$$, for $$\displaystyle e = 1,2.$$

3.) Comp. $$\tilde{\underline{K}}^e={\underline{L}^e}^T\underline{k}^e\underline{L}^e$$, for $$e = 1,2.$$ Compare to 1).

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

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

Solution
When using quadratic Lagrange element basis functions for this case we will employ 1-D quadratic elements, i.e. elements consisting of 3 nodes each. Figure 1 below shows the element structure for the case when nel=2.



Element Basis Functions and Their Derivatives
The quadratic Lagrange element basis functions are

Note that these functions take these values only within their corresponding element and are zero everywhere else.

Since the discretization used is uniform these basis functions may be simplified knowing the locations of the nodes. The simplified basis functions are

Then it follows that the derivatives of these functions are

Constructing the Global Stiffness Matrix
From equation [[media:Fe1.s11.mtg30.djvu|(1)p.30-5]]

Note: This matrix denotes the contribution due to element 1 to the global stiffness matrix. Since the basis functions $$b_4(x)$$ and $$b_5(x)$$ are zero within this element, all the terms in the fourth and fifth rows and columns are equivalently zero.

Similarly for the contribution from element 2:

Note: This matrix denotes the contribution due to element 2 to the global stiffness matrix. Since the basis functions $$b_1(x)$$ and $$b_2(x)$$ are zero within this element, all the terms in the first and second rows and columns are equivalently zero.

Now, constructing the global stiffness matrix using equation [[media:Fe1.s11.mtg30.djvu|(2)p.30-4]]

Computing the Element Stiffness Matrices
For the purpose of comparison we will now compute the global stiffness matrix by computing the element stiffness matrices and then using scatter matrices to scatter the coefficients into the proper place in the global stiffness matrix.

From equation [[media:Fe1.s11.mtg31.djvu|(1)p.31-3]]

A scatter matrix is a matrix that relates the local element degrees of freedom to the global degrees of freedom as seen in equation [[media:Fe1.s11.mtg31.djvu|(3)p.31-1]]

Due to the similarity of the structure of the element stiffness matrix and the element stiffness contribution matrix previously calculated, it can be seen that the element stiffness matrices are equivalent to the non-zero sub-matrices of the stiffness contribution matrices. i.e.

And by inspection the scatter matrices are found to be

Alternate Computation of the Global Stiffness Matrix
The element contribution to the global stiffness matrix is then computed via equation [[media:Fe1.s11.mtg31.djvu|(2)p.31-1]]

Now, as before,

Note:    $$\displaystyle (Eq. 6.2.21) = (Eq. 6.2.11)$$

QLEBF Used for nel=3
The following MATLAB code was used to generate the QLEBF used for nel=3, shown in figure 2 below



Numerical Solution
This problem was then solved using the QLEBF and increasing the number of elements, $$ nel = 2,4,6...$$ until the error converged below 10-6. The following MATLAB code was used to develop figures 3 and 4 below, showing a comparison between the numerical and exact solution and the error convergence, respectively.





Solved by William Kurth.

=Problem 6.9: Reproduce F&B Examples 8.3, 8.4 and Solve F&B Problem 8.6=

Objectives
We were asked to use a finite element analysis software or routine to:

1.) Reproduce example 8.3, FB, p. 200

2.) Reproduce example 8.4, FB, p. 205

3.) Use the same codes to solve problem 8.6, FB, p. 212

Solution
The MATLAB finite element codes developed by F&B, corresponding to chapter 8, were used to reproduce and solve the following problems, these codes are available at | the FB companion site and are also shown below.

The basic procedure for using these codes is as follows:


 * define the problem parameters, such as material properties, boundary conditions, mesh and element structure in an appropriate input file (i.e. input_file_64ele.m)
 * specify desired plots created in this input file
 * define the problem geometry, nodal coordinates, and element connectivity matrix either in this input file or in the mesh2d.m subroutine
 * specify in the preprocessor.m file which input file to load
 * if any new variables are introduced, add them to the include_flags.m file to establish them as global variables
 * run the heat2d.m file to solve the problem and produce the plots specified in the input file

The function assembly.m assembles passed element stiffness and force matrices into the current global stiffness and force matrices. Called iteratively for each element.

BmatHeat2D.m computes the basis function derivative matrix.

gauss.m gets Gauss points and weight factors.

get_flux.m calculates the flux vectors.

heat2d.m calls the necessary functions and .m files to carry out the computations.

heat2Delem.m calculates the element stiffness and force matrices as well as the elemental source term.

include_flags.m initializes global variables

input_file_1ele.m is the input file configured to compute example 8.3 with 1 element.

input_file_16ele.m is the input file containing parameters to solve example 8.3 with 16 elements.

input_file_64ele.m is the input file containing parameters to solve example 8.3 with 64 elements.

mesh2d.m constructs the nodes and mesh for example 8.3.

NmatHeat2D.m constructs the matrix of basis functions.

plotmesh.m is used to plot the mesh/element structure.

postprocess.m displays the calculated results.

preprocessor.m specifies which input file to load.

setup_ID_LM.m sorts the nodes which are constrained on the essential boundary condition from those which are not.

solvedr.m partitions the K and F matrices and solves for the degrees of freedom.

src_and_flux.m computes the nodal flux vectors.

Example 8.3
Consider a plate with geometry shown in figure 1 below. The conductivity is isotropic with $$D = k\bigl[\begin{smallmatrix} 1&0\\0&1\end{smallmatrix}\bigr]$$ with k = 5 W°C-1. The temperature is prescribed at T=0 at the bottom and right edges. Heat fluxes $$ \bar{q} = 0$$ and $$ \bar{q} = 20 Wm^{-1}$$ on the right and top boundaries, respectively. A constant heat source of $$\displaystyle s = 6 Wm^{-2}$$ is applied over the plate. Solve the problem using 16 quadrilateral finite elements.



The accompaniment codes come already configured to solve this problem, one must only simply specify in the preprocessor.m file that the input_file_16ele.m is the input file to use (comment out the others and save the file). Running heat2d.m displays these results:





Example 8.4
Consider a plate with geometry shown in figure 4 below, with prescribed temperature at x = ± b and prescribed flux at y = ± b.



The essential boundary conditions are given by

The natural boundary conditions on $$\Gamma_q$$ are $$(\bar{q}=-k\mathbf{n}^T\nabla T)$$

Lastly the source term is defined as

The MATLAB code used for example 8.3 was adapted to solve this problem by rewriting the input file. The results are shown below.









Problem 8.6
Consider a chimney constructed of two isotropic materials: dense concrete (k = 2.0 W°C-1 ) and bricks (k = 0.9 W°C-1 ). The temperature of the hot gasses on the inside surface of the chimney is 140 °C, whereas the outside is exposed to the surrounding air, which is at T = 10 °C. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney crosssectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials. Analyze the problem with 2 x 2, 4 x 4 and 8 x 8 quadrilateral elements for 1/8 of the problem domain. A 2 x 2 finite element mesh is shown in Figure 8.23. 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.



The previous MATLAB code was once again adapted to solve the given problem. Below are the results for the cases of nel = 4, 16, and 64.













Solved by William Kurth.