User:Eml5526.s11.team4/HW7

= Problem 7.1 - Static and Dynamic 2D G2DM1.0/D1 using LIBF and LLEBF with uniform and non-uniform mesh =

Problem Statement
Part 1: Static (Steady State)- 1a.) Use 2-D LIBF (Lagrange Interpolation Basis Functions) with n = m = 2, 4, 6, 8, ... to solve G2DM1.0/D1 for $$ u^{h} $$ until accuracy reaches $$ 10^{-6} $$ at the center (x,y) = (0,0). 1b1.) Repeat this problem using 2D LLEBF (Linear Legrange Element Basis Functions) with a uniform mesh. 1b2.) Repeat this problem using 2D LLEBF (Linear Legrange Element Basis Functions) with a non-uniform mesh. Recall G2DM1.0/D1: $$\Omega = \overline{w}$$ = Biunit Square $$\mathbf{K} = \mathbf{I}$$ (Identity matrix) $$ f(x) = 1 $$ $$\frac{\partial u}{\partial t}=0$$ Part 2: Dynamic (Transient)- 2a.) Use 2-D LIBF (Lagrange Interpolation Basis Functions) with n = m = 2, 4, 6, 8, ... to solve G2DM1.0/D1 for $$ u^{h} $$ until accuracy reaches $$ 10^{-6} $$ at the center (x,y) = (0,0). Consider: $$ f(\mathbf{x},t)=0, \forall \; \mathbf{x} \; \varepsilon; \Omega, \forall \; t> 0 $$ $$ g(\mathbf{x},t)=2 $$ on $$\Gamma_{g}=\partial\Omega $$ 2b1.) Use 2-D LIBF (Lagrange Interpolation Basis Functions) with n = m = 2, 4, 6, 8, ... to solve G2DM1.0/D1 for $$ u^{h} $$ until accuracy reaches $$ 10^{-6} $$ at the center (x,y) = (0,0). Consider: $$ f(\mathbf{x},t)=1, \forall \; \mathbf{x} \; \varepsilon \; \Omega, \forall \; t> 0 $$ $$ g(\mathbf{x},t)=2 $$ on $$\Gamma_{g}=\partial\Omega $$ 2b2a.) Repeat problem 2b1 (shown above) using 2D LLEBF (Linear Legrange Element Basis Functions) with a uniform mesh. 2b2b.) Repeat problem 2b1 (shown above) using 2D LLEBF (Linear Legrange Element Basis Functions) with a non-uniform mesh.

Part a) The weak form
To get to the weak form we multiply by an arbitrary function, $$w$$, and integrate over the domain. This is the weighted residual form.

From the weighted residual form, we use partial integration and relax the requirement on the functions $$u$$ and $$w$$ to be once differentiable, instead of twice differentiable in the weighted residual form and strong form.

Recall partial integration

We use this on 7.4.5 where dv is the second order derivatives. Important, the boundary term is not a simple evaluation! It is a line integral along the boundary. In general the integration will be one dimension less than the dimensionality of the problem (A surface in a 3D problem, and line in 2D problems).

Now we expand the boundary terms

We now apply the boundary condition to our domain

Where the essential boundary terms have been removed by requiring $$w=0$$ at the essential boundary.

We note that on $$\Gamma_h$$, $$\underline{q}\cdot\underline{n} = h(x,t)$$ and that $$q_i=-\frac{\partial u}{\partial x_j}$$ so that we only have normal derivatives to the boundary, then

For the mixed boundary, $$\Gamma_H$$, $$ \underline{q}\cdot\underline{n} = H(u-u_\infty)$$

Finally, all boundary terms evaluated, we have the weak form where $$w=0$$ on $$\Gamma_g$$

Part b) Discrete weak form
We now discritize, where summation is implied by the indices

Into 7.4.12

Factoring out non-spatial terms, and removing $$c_i$$

Part c) Solve using Lagrange basis functions
Looking at 7.4.14, we note that all spatial terms may now be evaluated. The unknown can be factored out onto the right hand side of the evaluated matrix. Doing so results in 7.4.15

Now, with a basis function $$b_i(x)$$ we can invert the matrix on the left hand side and solve for $$d_j$$.

We use Lagrange basis functions, but need them in two dimensions. To do this we take the tensor product of two one dimensional Lagrange basis functions, as follows.

As an example, consider a first order basis function in two dimensions:

Note that this can be extended in order and in dimensions. To extend to three dimensions one can use $$b_i=b_l\otimes b_m \otimes b_o$$. To increase the order of interpolation increase the dimensions of each basis function before the tensor product.

A code which can do this in general is provided in MATLAB below.

Here we introduce a standard naming scheme that accounts for matrices that arise in the use of FEM codes. We first use it, then explain it. We write 7.4.15 as 7.4.18.

Where in B200, M211, and M222 the first letter denotes the type of matrix, the second the dimensionality, and the following integers the order of the derivatives in the expression. So for B200, its a boundary matrix of two dimensions with no derivatives. For M211 it's a mass matrix of two dimensions with a first derivative in the first direction times a first derivative in the first direction. For M222 it's a mass matrix of two dimensions with a first derivative in the second direction times a first derivative in the second direction.

A MATLAB code to calculate these has also been developed. This follows the execution of the previous code and works for any dimensionality and any order of interpolation.

For the case of 2-D and first order interpolation functions (example used above).

Now, substituting for the numeric values in 7.4.20

Now the essential boundary condition must be handled. This is done by assigning $$b_j$$ along the essential boundary and solving the remaining equations. One way to do this is to 'strike' the matix elements along essential values, and assign the right hand side to the essential value. That way $$b_i$$ will be forced to be the essential boundary condition value.

A method for doing this is presented in the MATLAB code below.

Results
Now, we increase the interpolation order and collect the value at the center point until we reach $$10^{-6}$$ accuracy (this is p-convergence instead of h-convergence). We will define our accuracy relative to our highest level solution.

Without using an analytic solution, we approximated the error by comparing with $$n=5$$. At $$n=4$$ we have reached $$10^{-4}$$ accuracy and expect the increase in precision to following the trend in the log plot below. With this assumption, we have reached $$10^{-6}$$ at $$n=5$$.

Plots
A convergence plot at the center and upper left (where $$\Gamma_h$$ and $$\Gamma_H$$ are applied)



The following are surface plots for the fitted interpolation polynomials of increasing order.



Complete Code Listing
The MATLAB code used, listed in full with comments. Note that the code 'caches' the compiled matrices as functions in the runtime directory. After the first run for a given order following runs are very fast.

Alternative Solution by Fangfang Zhu:
First derive the equilibrium equation for an arbitrary choose domain showns as:



As showns in the above figure,the equilibrium of heat can be written as:

Using the Taylor series expansion we can write:

So the equilibrium equation can be rewritten as:

Dividing dxdy both side yields:

Written in compact form:

Consider the constitutive equation for 2D in heat transfer. The 2D Fourier Law for 2D space can be written as: $$\displaystyle \vec{q}=-[k]\nabla T$$

Written in cartisan matrix form:

For isotropic material we have $$\displaystyle {{k}_{xy}}={{k}_{yx}}=0$$ and $$\displaystyle {{k}_{xx}}={{k}_{yy}}=k$$ So we can write:

or

Thus the strong form for this problem in index notation can be written as:

Following the routine of deriving weak form we can easily get the weak form for this problem:

Let

The weak form can be written as:

Descritization process: Set the approximated solution as a linear combination of shape functions $$\displaystyle {{N}_{i}}(x)$$ as:

And the approximated weighting function is alinear combination of shape functions $$\displaystyle {{N}_{i}}(x)$$ as:

Where $$\displaystyle \tilde{K}={{\int\limits_{D}{\left[ B \right]}}^{T}}\cdot \left[ K \right]\cdot \left[ B \right]$$

Similoarily:

Where $$\displaystyle \tilde{M}=\int\limits_{D}{{{N}^{T}}\cdot \rho c\cdot N}dD$$

where $$\displaystyle {{F}_{i}}=\left[ \int\limits_{-hTN_{i}^{T}ds}+\int\limits_{{{q}_{0}}N_{i}^{T}ds}+\int\limits_{h{{T}_{a}}{{N}_{i}}ds+\int\limits_{D} – f{{N}_{i}}}dD \right]$$

Thus we can write the discrete weak form:

The temprature distribution by Fangfang's code