Finite elements/Solution of heat equation

From Wikiversity
Jump to navigation Jump to search

Finite element solution for the Heat equation[edit | edit source]

Let us now try to create a finite element approximation for the variational initial boundary value problem for the heat equation . This problem can be stated as


Let be a finite dimensional approximation to and let be a finite dimensional approximation to . Let the weighting functions be such that on . Similarly, any other function also goes to zero on the boundary .

Assume that the trial solutions can be represented as

The additional quantity results in the satisfaction of the boundary condition on .

If we substitute these quantities into the variational initial boundary value problem, we get the Galerkin formulation. The steps in this process are shown below.

Approximate IBVP[edit | edit source]

The approximate form of the variational IBVP is

Replacing with we get

After expanding and rearranging terms, we get

In compact notation

Similarly, the initial condition can be written as

Substituting for and expanding out terms, we have

In compact form,

Therefore the approximate form of the variational BVP can be written as


Note that the derivatives with respect to time are still continuous in this approximate variational BVP.

Finite element approximation[edit | edit source]

We now discretize our domain into elements where and is the total number of elements. Nodes may exist anywhere within the element domains. But the most common locations are at the element vertices and along the interelement boundaries.

Let the set of global node numbers be

Let the set nodes at which a temperature (essential BC) is prescribed be

Then the set of nodes at which a temperature is to be determined is the complement

Let be the number of nodes in . Then the number of equations to be solved is also . See Figure 1 to get an idea about what these sets of nodes refer to.

Figure 1. FEM discretization for the heat conduction problem.

As we have seen before, a typical weighting function is assumed to have the form

Note that only if for every . On the other hand for every node in .

Similarly, a typical trial solution can be written as

where is the unknown temperature at node at time .

We also often represent the approximate form of the essential BCs in a similar way, i.e.,

Substitute equations (42), (43), and (44) into the first of equations (40). You will get

Assuming that does not change with time, we can take the constants and time-dependent parameters outside the integrals to get

Using the usual argument about the arbitrary nature of , we get

In compact notation,

where , and .

Let us rename parts of the above equation as follows:

Then we get

In matrix form,

As before, we can break up the global matrices into sums over elements. Thus,

where

and is the number of nodes in an element.

We can do the same for the initial conditions. However, a more common approach is to specify the values of directly at the nodes instead of going through an assembly process.

Computing M, K, f[edit | edit source]

The finite element system of equations at time :

Isoparametric Map[edit | edit source]

Isoparametric map.

Coordinate Transformation[edit | edit source]

Integrating Stiffness Matrix[edit | edit source]

Stiffness matrix terms:

Index notation:

Expanded out (in 2D) (assume ):

Transformation[edit | edit source]

Chain Rule:

Matrix notation:

Inverse Transformation[edit | edit source]

with

Returning to the integration of the stiffness matrix terms:

In parent element coordinates:

Gaussian Quadrature[edit | edit source]

Gaussian Quadrature Example[edit | edit source]

Find the constants Co, C1, and X so that the quadrature formula

This has the highest possible degree of precision.

Solution.

Since there are three unknowns, Co, C1 and X1, we will expect the formula to be exact for

Thus

Equation 2 and 3 will yield.

Hence

Now,

And

Thus the degree of the precision is 2

Robin boundary conditions[edit | edit source]

A similar approach can be used when w:Robin boundary conditions are needed to be applied. Let us assume that the Robin boundary conditions take the form

Recall from the previous section that the boundary term in the weak form of the heat equation is

where

.

Let us assume, for simplicity, that the material is isotropic. In that case becomes a scalar quantity and we can write

Now, we can write the Robin boundary condition as

Using this expression we can get of the flux term in to get

Then the boundary term takes the form

If we ignore the contributions to the trial function from the Dirichlet boundary conditions we can write

Plugging these expressions into the boundary term leads to

Now recall that spatially discretized equation for transient heat conduction can be expressed as (in simplified form; after getting rid of contributions due to Dirichlet boundary conditions and with a scalar )

Substituting the last term on the right hand side with the new expression for the boundary term, we have

Invoking the arbitrariness of and rearranging, we get

Define

Then the finite element system of equations is

or, in matrix form

Note that the matrix and the vector are different from the case where there are no Robin boundary conditions. The boundary terms in the discretized system of equations can be computed in the usual manner and will add terms to the diagonal of the matrix.