Finite Element Method



Finite Element Method

Consider an ODE of the form:

[pic]

with the following Boundary Conditions:

[pic] and [pic]

1.) Substitute for U an unspecified trial function into governing equation, i.e.

[pic]

[pic]

where

n = Number of nodes in system (or, locally, in element).

Uj = Unknown coefficients (to be determined), which are not a function of space, in general.

Nj = Known (user-specified) Basis function, which are a function of space, but not a function of time, in general.

2.) Force the Residual to zero in the global sense, i.e. satisfy the governing equation in the "weak form".

a.) Multiply by a Weighting function:

[pic]

[pic]

where Wi = Known (user-specified) Weight function.

b.) Integrate (take inner product) over global domain and force to zero.

Definition of an Inner Product: [pic]

[pic] Implies that W is orthogonal to R

< ( ) > = [pic]

(dimensionally dependent)

[pic]

If you can solve this equation then you have solved the governing equations, at least in the weak form.

3.) Discretize global domain into subdomains (Elements) such that Σ (Elements) = contiguous, non-overlapping representation of the global domain.

[pic]

4.) Select a basis function for [pic] - a few helpful hints in this selection are

a.) select an orthogonal series basis function

b.) choose a computationally efficient series

Ex. Lagrange Polynomials applied locally on each element

[pic]

[pic]

Assume Linear Element (N = 2)

[pic]

Assume Quadratic Element (N = 3)

[pic]

[pic]

5.) Integration by parts:

Required for [pic] order PDEs using linear basis functions; Provides a direct means to apply Boundary Conditions

[pic]

[pic]

more generally,

[pic]

Adding in the other components of the ODE yields:

[pic]

6.) Select Weighting functions

Ex: Choose "Galerkin" [pic]

(However, for clarity retain the W symbol. It will help when assembling matrices later.)

Note that for each weighting function, Wi, one must assemble all contributions from each basis function Nj.

7.) Assemble Matrices – consider only 2 elements, and assemble all components associated with Row I (for Wi):

[pic]

One row of equations for each Wi

Note: the local definitions of Trial and Weighting functions

[pic] for all [pic]

= 0 for all [pic]

[pic] for row i only contributions from [pic]

Left Elem Right Elem

[pic] = [pic] + 0 = [pic]

[pic] = [pic] + [pic] = [pic]

[pic] = 0 + [pic] = [pic]

[pic] = [pic] + 0 = [pic]

[pic] = [pic] + [pic] = [pic]

[pic] = 0 + [pic] = [pic]

[pic] = [pic] + [pic] = [pic]

Assembly of row I

[pic]

If uniform spacing then h1 = h2 = h

[pic]

[pic] [pic] + [pic] Simpson’s Rule = g

Note the similarities to F.D.

In general, we are forming an expression

[pic]

[pic]

In a FE code: Loop over each element and assemble the local (elemental) contributions for the Lhs matrix [A] and the Rhs vector.

After this loop is completed, then (and only then) apply the boundary conditions to the set of equations.

FE codes contain a lot of bookkeeping. One common mapping array is the element connectivity (or Incidence) list. It has the form:

IN(K,L) where

L = Element Number, L = 1, NE ( and NE = Number of Elements)

K = Local Node Number, K = 1, 2, … to the number of nodes in an element.

IN(K,L) = global node number ( maps global node number to the local node number within element L)

Consider the following 1-D example:

The Element number can have significance if using a frontal matrix solver.

The Node numbering can have significance if using a banded matrix solver.

Node and Element numberings have less significance if using a sparse, iterative matrix solver.

The mapping for this example is:

IN(1,1) = 3 IN(2,1) = 1

IN(1,2) = 4 IN(2,2) = 6

IN(1,3) = 1 IN(2,3) = 4

IN(1,4) = 7 IN(2,4) = 2

IN(1,5) = 5 IN(2,5) = 7

IN(1,6) = 6 IN(2,6) = 5

The FE approximation to the governing equation is accumulated (summed) in 3 nested loops (L, I, and J)

Do L = 1, NE (Loop over all elements)

h = x(in(2,L)) – x(in(1,L)) (Calculate element length)

Do I = 1, 2 (Loop for all Weighting Functions)

dWdx = (1/h)

if (I = 1) dWdx = - dWdx

Do J = 1, 2 (Loop for all Trial Functions)

dNdx = (1/h)

if (J = 1) dNdx = - dNdx

term1 = - dNdx * dWdx * h

term2 = f * h/6

if (I = J) term2 = 2 * term2

AL(J, I) = term1 + term2 (Local 2x2 construction)

Irow = IN(I,L) (Global Mapping)

Jcol = IN(J,L)

AG(Jcol, Irow) += AL(J, I) (Global Matrix Construction)

End Do (End of J loop)

Rhs(Irow) += g * h/2 (Rhs Construction)

End Do (End of I loop)

End Do (Finished with all Elements)

Everything done at the element level!

Easy to automate

a) One Element (problem-specific)

b) Assembly (problem-independent)

8) Add in the Boundary Conditions

a.) Type 1 BC : Satisfy Exactly U(0) = 1

(Note this is at global node number 3 for our example.)

One less unknown in algebraic system.

Remove row_3 and then move column_3 of matrix to Rhs since it is known (not usually done, but possible)

[pic]

This process is somewhat cumbersome. It can change the bandwidth or structure of your system. It can cause other solution-solving anomalies if one is not careful.

Cleaner process:

Remove row 3 completely, i.e. place zeros in all entries.

Then place a 1 on diagonal of [A]

Put the solution (i.e. U(0)=1) into the Rhs

[pic]

b.) Type II B.C. Satisfy approximately in {BCs} vector.

Recall: [pic] and for node 2 (in our example) the flux at the boundary is 5. Also note: [pic]

At the boundary node Wi = 1 for I = 2 (global) and Wi = 0 for all other nodes.

Therefore to apply Type II Boundary Conditions in our example:

Rhs(2) = Rhs(2) – 5

The entire formulation is complete. Call a matrix solver and obtain the solutions for Uj.

-----------------------

[pic]

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download