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

=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:



Comparison with using distorted quadrilateral meshes
=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.

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.

Below is the analysis results: