Contents

Chapter 1
A 1-D Example

1.1   Introduction

The finite element method (FEM) is a numerical method to solve differential equations. There are several basic steps to solving a problem by the FEM method. These steps are:

1. Formulation of the problem into a differential or integral equation.
2. Development of the weak formulation of the problem.
3. Finite element approximation of the domain and variables.
4. Assembling the stiffness matrix and load vector and solving the problem.
5. Post processing (analysis and presentation of the results).

1.2   Steady State Heat Conduction in a Bar

We will proceed to describe the steps in the finite element method by the aid of an example. Consider the problem of heat conduction through a bar.

1.2.1   Step 1: Formulation of the problem

The problem under consideration is that of the steady state flow of heat through a bar.

Figure 1.1 shows how heat flows through an element of this bar. Balancing the heat entering and exiting this element gives

q(x) +f(x) \triangle x = q(x+\triangle x ) ,

(1.1)

where q is the heat flux along the bar and f(x) is the heat flow through the lateral surfaces per unit length of the bar. Reorganizing and taking the limit as \triangle x goes to zero yields

dq dx = f(x).

(1.2)

Introducing Fourier's law of heat conduction ( q = - k[(d q)/dx]) into this equation results in

-k d2 q dx2 =f(x) ,

(1.3)

where q is the temperature and k is the coefficient of thermal conduction.

Since the problem is a second order ordinary differential equation one can prescribe two conditions. We will look at solving the above differential equation subject to any one of the following boundary conditions:

• q = q_a at x=a and q = q_b at x=b ,
• q = q_a at x=a and q=-k[(d q)/dx] = q_b at x=b ,
• q=-k[(d q)/dx] = q_a at x=a and q = q_b at x=b , or
• q=-k[(d q)/dx] = q_a at x=a and q=-k[(d q)/dx] = q_b at x=b .

1.2.2   Step 2: Formulation of the weak form

The solution to this problem is a temperature field which satisfies the differential equation

-k d2 q dx2 =f(x)

(1.4)

at every point in the bar. We will relax this requirement by replacing this equation by

- ó õ   k q   d2 q dx2 dx = ó õ   q   f(x) dx ,

(1.5)

where  is any region of the bar and [`(q)] is a test function. If one requires that this equation hold for all possible test functions [`(q)], it will be possible to get the original differential equation from this integral equation. Any point in the bar can be isolated by selecting a test function which is only positive and is nonzero only in a small region around that point. The aim is not to require the integral equation to hold for all possible [`(q)] functions, but to select a set of test functions and require the integral equality to hold for this set of functions. Hence, we do not require the differential equation to be exactly satisfied at all points, but require the temperature field to only satisfy the integral equation for a select number of test functions (we have weakened the original equation).

We will now use integration by parts to get the weak form of the problem as

-[k q   d q dx ]x=ab + ó õ b a  k d q   dx d q dx dx = ó õ b a  q   f(x) dx ,

(1.6)

subject to the appropriate boundary conditions. The application of integration by parts also introduces a weakening of the problem since we are no longer required to have temperature fields which have a second derivative. It is now only required that the temperature field have a first derivative as can be seen from equation (1.6).

1.2.3   Step 3: Finite element approximation

We will use two node line elements to approximate both the temperature field and the test functions. The length of the bar will be partitioned into ne elements as shown in Figure 1.2. Each element has an element number and one node at each end. Nodes are given two numbers. The global node number is one which is unique for each node and which distinguishes it from all other nodes. For each element there is also a local numbering system. For a two node element the local nodes are distinguished in the element by designating a local node number of 1 to one node and 2 to the other node.

fig1-2.wmf

The temperature and the test function will be approximated over each element by

q = qe1N1(x)+qe2N2(x) = é ë N1(x) N2 (x) ù û é ê ë qe1 qe2 ù ú û ,

(1.7)

and

q   = q   e 1  N1(x)+ q   e 2  N2(x) = é ê ë q   e 1  q   e 2  ù ú û é ê ë N1(x) N2 (x) ù ú û ,

(1.8)

where q^e_i is the value of q at local node i, [`(q)]<sup>e</sup><sub>i</sub> is the value of [`(q)] at local node i, and N_i(x) are known functions of x. The N_i are known as shape functions. The shape functions are selected such that N₁=1 at local node 1, N₁=0 at local node 2, N₂=0 at local node 1, N₂=1 at local node 2. For a two node line element of length l_e and local coordinate s as shown in Figure 1.3, the shape functions will be

N1 = 1- s le ,

(1.9)

and

N2 = s le .

(1.10)

The shape functions are selected such that the sum of the shape functions at each point will add up to unity. It must be noted that the shape functions are different for different elements, even though this is not made explicit in the notation.

The two integrals in the weak formulation of the problem, given in equation (1.6), can each be written as a sum of integrals, each integral being over one element. This can be written as

ó õ b a  k d q   dx d q dx dx = ne å e=1  ó õ We  k d q   dx d q dx dx,

(1.11)

and

ó õ b a  q   f(x) dx = ne å e=1  ó õ We  q   f(x) dx ,

(1.12)

where W_e is the domain of element number e.

Over each element one can use the above approximations for the temperature and test function to get

ó õ We  k d q   dx d q dx dx = ó õ le 0  k d[ q   e 1  N1 + q   e 2  N2] ds d[qe1 N1 +q2 N2] ds ds ,

(1.13)

and

ó õ We  q   f(x) dx = ó õ le 0  [ q   e 1  N1 + q   e 2  N2] f(x) ds ,

(1.14)

since x=x₁+s, where x₁ is the location of the first node in the element. Since q^e₁, q^e₂, [`(q)]<sup>e</sup><sub>1</sub>, and [`(q)]^e₂ are constants, one can write

d q   dx = d[ q   e 1  N1 + q   e 2  N2] ds = é ê ë q   e 1  q   e 2  ù ú û é ê ê ê ê ê ê ê ê ê ê ë d N1 ds d N2 ds ù ú ú ú ú ú ú ú ú ú ú û ,

(1.15)

and

d q dx = d[qe1 N1 +qe2 N2] ds = é ê ê ê ë d N1 ds d N2 ds ù ú ú ú û é ê ë qe1 qe2 ù ú û .

(1.16)

Substitution of these relations into the integrals of equations (1.13) and (1.14) gives

ó õ We  k d q   dx d q dx dx = k é ê ë q   e 1  q   e 2  ù ú û ó õ le 0  é ê ê ê ê ê ê ê ê ê ê ë d N1 ds d N2 ds ù ú ú ú ú ú ú ú ú ú ú û é ê ê ê ë d N1 ds d N2 ds ù ú ú ú û ds é ê ë qe1 qe2 ù ú û ,

(1.17)

and

ó õ We  q   f(x) dx = é ê ë q   e 1  q   e 2  ù ú û ó õ le 0  é ê ë N1 N2 ù ú û f(x) ds .

(1.18)

These equations can be written as

ó õ We  k d q   dx d q dx dx = [ q   e   ]T [Ke][qe],

(1.19)

and

ó õ We  q   f(x) dx = [ q   e   ]T [fe] ,

(1.20)

where

[ q   e   ] = é ê ê ê ê ë q   e 1  q   e 2  ù ú ú ú ú û ,

(1.21)

[qe] = é ê ë qe1 qe2 ù ú û ,

(1.22)

[Ke] = é ê ë K11e K12e K21e K22e ù ú û = k ó õ le 0  é ê ê ê ê ê ê ê ê ë ( dN1 ds )2 dN1 ds dN2 ds dN2 ds dN1 ds ( dN2 ds )2 ù ú ú ú ú ú ú ú ú û ds,

(1.23)

[fe] = é ê ë f1e f2e ù ú û = ó õ le 0  é ê ë N1 N2 ù ú û f(x) ds .

(1.24)

The matrix [K^e] is known as the element stiffness matrix and the matrix [f^e] is known as the element load vector. For the shape functions in (1.9) and (1.10) we have [(dN₁)/ds] = -[1/(l_e)] and [(dN₂)/ds] = [1/(l_e)]. Therefore, the element stiffness matrix can be written as

[Ke] = k ó õ le 0  é ê ê ê ê ê ê ê ê ë 1 le2 - 1 le2 - 1 le2 1 le2 ù ú ú ú ú ú ú ú ú û ds = é ê ê ê ê ê ê ê ê ë k le - k le - k le k le ù ú ú ú ú ú ú ú ú û .

(1.25)

To evaluate the element load vector one needs to know f(x). We will later calculate this for a particular example.

The weak form of the problem can now be written as

[ q   q]x=ab + ne å e=1  [ q   e   ]T[Ke][qe] = ne å e=1  [ q   e   ]T[fe],

(1.26)

where Fourier's law is used to replace the derivative of temperature in the first term.

1.2.4   Step 4: Assembling the stiffness matrix and load vector

The next step is to organize the problem in the form of a set of algebraic equations. Consider the example of a problem with two elements and three nodes as shown in Figure 1.4. The two summations in equation (1.26) can be written as

2 å e=1  [ q   e   ]T[Ke][ q   e   ] = é ê ë q   1 1  q   1 2  ù ú û é ê ë K111 K112 K121 K122 ù ú û é ê ë q11 q12 ù ú û + é ê ë q   2 1  q   2 2  ù ú û é ê ë K211 K212 K221 K222 ù ú û é ê ë q21 q22 ù ú û

= é ê ë q   1  q   2  q   3  ù ú û é ê ê ê ë K111 K112 0 K121 K122+K211 K212 0 K221 K222 ù ú ú ú û é ê ê ê ë q1 q2 q3 ù ú ú ú û ,

(1.27)

and

2 å e=1  [ q   e   ][fe] = é ê ë q   1 1  q   1 2  ù ú û é ê ë f11 f12 ù ú û + é ê ë q   2 1  q   2 2  ù ú û é ê ë f21 f22 ù ú û

= é ê ë q   1  q   2  q   3  ù ú û é ê ê ê ë f11 f12+f21 f22 ù ú ú ú û

(1.28)

Therefore, for this problem the weak formulation given in equation (1.29) can be put in the form

[ q   3  q3- q   1  q1] +[ q   ]T [K][q] = [ q   ][f],

(1.29)

where [K] is known as the global stiffness matrix, [f] is known as the global load vector, and

[ q   ] = é ê ê ê ê ê ê ê ê ë q   1  q   2  q   3  ù ú ú ú ú ú ú ú ú û ,

(1.30)

[q] = é ê ê ê ë q1 q2 q3 ù ú ú ú û ,

(1.31)

[K] = é ê ê ê ë K11 K12 K13 K21 K22 K23 K31 K32 K33 ù ú ú ú û = é ê ê ê ë K111 K112 0 K121 K122+K211 K212 0 K221 K222 ù ú ú ú û ,

(1.32)

[f] = é ê ê ê ë f1 f2 f3 ù ú ú ú û = é ê ê ê ë f11 f12+f21 f22 ù ú ú ú û .

(1.33)

To illustrate the way in which one arrives at the final set of equations we will write the expanded form of equation (1.29) in the form

[ q   3  q3- q   1  q1] + é ê ë q   1  q   2  q   3  ù ú û é ê ê ê ë K11 K12 K13 K21 K22 K23 K31 K32 K33 ù ú ú ú û é ê ê ê ë q1 q2 q3 ù ú ú ú û = é ê ë q   1  q   2  q   3  ù ú û é ê ê ê ë f1 f2 f3 ù ú ú ú û .

(1.34)

The test functions [`(q)] are arbitrary functions which we will select in the manner shown in Figure 1.4. Three test functions are selected since we will have three unknowns in this problem ( q<sub>1</sub>, q<sub>2</sub>, q<sub>3</sub>; q<sub>1</sub>, q<sub>2</sub>, q<sub>3</sub>; q<sub>1</sub>, q<sub>2</sub>, q<sub>2</sub>; or q<sub>1</sub>, q<sub>2</sub>, q<sub>2</sub>). Each test function provides one algebraic equation for the unknowns. As can be seen from Figure 1.4, the first test function has [`(q)]₁ = 1, [`(q)]<sub>2</sub> = 0, and [`(q)]₃ = 0. The second test function has [`(q)]<sub>1</sub> = 0, [`(q)]₂ = 1, and [`(q)]<sub>3</sub> = 0. The third test function has [`(q)]₁ = 0, [`(q)]<sub>2</sub> = 0, and [`(q)]₃ = 1.

The following three equations are obtained from the substitution of the selected test functions into equation (1.34).

-q1 + K11q1 +K12 q2 +K13q3 = f1 ,

K21q1 +K22 q2 +K23q3 = f2 ,

q3+K31q1 +K32 q2 +K33q3 = f3 .

(1.35)

These equations can now be solved for the three unknowns. If q₁ and q₃ are specified, then this system can be written as

é ê ê ê ë K11 K12 K13 K21 K22 K23 K31 K32 K33 ù ú ú ú û é ê ê ê ë q1 q2 q3 ù ú ú ú û = é ê ê ê ë f1+q1 f2 f3- q3 ù ú ú ú û .

(1.36)

This would be the problem which one would solve if this was a well posed problem. It turns out that this particular set of boundary values are not proper. Since the temperature is not fixed at any point in the bar, the temperature profile can slide up and down. An examination of the stiffness matrix shows that its determinant is zero (i.e., the system is singular).

If q₁ and q₃ are specified, then the middle equation can be used to find q₂, and the first and last equation are used to find the heat flux at the two ends. The, mixed type boundary conditions can also be solved in a similar manner.

Before we proceed it is important to state that in practice the boundary conditions are not imposed as was presented above. Theoretically, if we have known values of temperature at the boundary we should eliminated these known values from the set of unknowns and introduce the unknown heat fluxes in the boundary term into the list of unknowns. To avoid this complex procedure, the penalty method is used for imposing the temperature boundary. This method modifies the stiffness matrix and load vector to fix the value of temperature to a given value at a boundary node. The method is based on the fact that for any boundary node i there will be an equation

Kii qi +     other    terms      = fi .

(1.37)

To impose the condition q_i = [^(q)] one can add a term pK_iiq_i to the left hand side of the equation and the term pK_ii[^(q)] to the right hand side of the equation to get

Kii qi + pKiiqi +    other    terms      = fi+pKii ^ q

(1.38)

where p is a large number. After dividing by pK_ii one gets

qi p + qi+     other    terms      pKii = fi pKii + ^ q

(1.39)

All term in this equation become small accept the two newly added ones. Therefore, the equation will become dominated by the equation q_i = [^(q)] and the other terms become unimportant. Since the method eliminates the importance of any other terms entering this equation, it will not be necessary to introduce the unknown heat flux into the list of unknowns.

In short, to impose a temperature boundary condition at global node i one can multiply K_ii in the stiffness matrix by (1+p) (i.e., add pK_ii to K_ii), and add pK_ii[^(q)] to f_i. To impose a heat flux boundary condition at global node j one adds q to f_j, where q is a heat flux directed into the bar.

1.3   An example

As an example consider the problem of the bar shown in Figure 1.5. The bar is of length l=1, is connected to a constant temperature thermal bath of temperature [^(q)]=2 from the left and heat is being drawn from it at a constant rate [^q]=1 from the right. The heat flow into the bar per unit length from the lateral surfaces is given by f(x)=x. The coefficient of thermal conduction k = 1.

Three two node line elements are selected as shown in Figure 1.5. The local stiffness matrices can be calculated using equation (1.25) and will be

[K1] = é ê ë 2 -2 -2 2 ù ú û ,

[K2] = é ê ë 4 -4 -4 4 ù ú û ,

[K3] = é ê ë 4 -4 -4 4 ù ú û .

The global stiffness matrix will, therefore, be

[K] = é ê ê ê ê ë 2 -2 0 0 -2 (2+4) -4 0 0 -4 (4+4) -4 0 0 -4 4 ù ú ú ú ú û .

The load vectors can be calculated using equation (1.24). Using equation (1.9) and (1.10), the shape functions for the first element are

N1 = 1-2x,         N2 = 2x .

Therefore,

f11 = ó õ [1/2] 0  (1-2x)x dx = 1 24 ,

f12 = ó õ [1/2] 0  2x2 dx = 1 12 .

The shape functions for the second element are

N1 = 1-4(x- 1 2 ),         N2 = 4(x- 1 2 ).

Therefore,

f21 = ó õ [3/4] [1/2]  [1-4(x- 1 2 ) ]x dx = 7 96 ,

f22 = ó õ [3/4] [1/2]  4(x- 1 2 ) x dx = 1 12 .

The shape functions for the third element are

N1 = 1-4(x- 3 4 ),         N2 = 4(x- 3 4 ).

Therefore,

f31 = ó õ 1 [3/4]  [1-4(x- 3 4 ) ]x dx = 5 48 ,

f32 = ó õ 1 [3/4]  4(x- 3 4 ) x dx = 11 96 .

The global load vector will become

[f] = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 1 24 1 12 + 7 96 1 12 + 5 48 11 96 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û .

The problem that we must solve is

é ê ê ê ê ë 2(1+p) -2 0 0 -2 6 -4 0 0 -4 8 -4 0 0 -4 4 ù ú ú ú ú û é ê ê ê ê ë q1 q2 q3 q3 ù ú ú ú ú û = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 1 24 +4p 5 32 3 16 11 96 -1 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û ,

where p is a large number and the stiffness matrix and load vectors have been modified to impose the boundary conditions.

For p=1 the solution will be

é ê ê ê ê ë q1 q2 q3 q3 ù ú ú ú ú û = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 7 4 71 48 167 128 13 12 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û = é ê ê ê ê ë 1.75 1.47916 1.30468 1.08333 ù ú ú ú ú û .

For p=10 the solution is

é ê ê ê ê ë q1 q2 q3 q3 ù ú ú ú ú û = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 79 40 409 240 979 640 157 120 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û = é ê ê ê ê ë 1.975 1.70416 1.52968 1.30833 ù ú ú ú ú û .

For p=100 the solution is

é ê ê ê ê ë q1 q2 q3 q3 ù ú ú ú ú û = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 799 400 259 150 4967 3200 1597 1200 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û = é ê ê ê ê ë 1.9975 1.72666 1.55218 1.33083 ù ú ú ú ú û .

The exact solution is

q = 2- x3 6 - x 2 ,

which yields

é ê ê ê ê ë q1 q2 q3 q3 ù ú ú ú ú û = é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë 2 83 48 199 128 4 3 ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú ú û = é ê ê ê ê ë 2 1.72916 1.55468 1.33333 ù ú ú ú ú û .

Comparison of the exact solution and the finite element method solution shows that for p=100 there is less than 1% error in the finite element solution.

Chapter 2
A 2-D Example

The following is an example of how to solve a two dimensional problem using the finite element method. The problem under consideration is the torsion of a noncircular member.

2.1   Step 1: The torsion problem

The differential equation for the torsion problem is given by

¶2 f ¶x2 + ¶2 f ¶y2 = -2Gq,

(2.1)

where f is a function of x and y and is called the potential function, G is the shear modulus, and q is the amount of twist per unit length of the bar.¹ To solve the torsion problem one must solve for the potential function f which satisfies this differential equation at every point in the cross section W and which also satisfies the boundary conditions

f = 0

(2.2)

on the boundary curve ¶W.

If one defines the gradient operator in the x-y plain as

Ñ = ^ i   ¶ ¶x + ^ j   ¶ ¶y ,

(2.3)

then one will have

Ñ°Ñ = ¶2 ¶x2 + ¶2 ¶y2 .

(2.4)

Therefore, one can write the problem as

Ñ°Ñf = -2 G q          in       W

(2.5)

and

f = 0           on      ¶W

(2.6)

In this problem the shear stress at any point of the cross section can be found from f using

txy = ¶f ¶z       and       txz = - ¶f ¶y .

(2.7)

The torque T required to crate a twist of q is given by

T = 2 ó õ W  fdA ,

(2.8)

where dA is the differential element of area.

2.2   Step 2: The weak form of the torsion problem

In a simple way one can describe the generation of the weak form of an equation as follows. If for every c one can show ca = cb, then one can conclude that a=b. The weak form of equation a=b in this case is ca=cb.

To formulate the weak form of the torsion problem we will multiply both sides of (2.5) by a arbitrary potential function [`(f)](x,y). This will give the relation

f   Ñ°Ñf = -2 G q f             in       W.

(2.9)

Since this equation must hold for every point in the cross section W, we can integrate this relation over any region  of the cross section to get

ó õ   f   Ñ°ÑfdA = -2 ó õ   G q f   dA .

(2.10)

This equation is regarded as the weak form of the differential equation (2.5). For this equation to be equal to equation (2.5), this relation must hold for every function [`(f)].

We will rewrite equation (2.10) in a more convenient form using first the relation

Ñ°( f   Ñf) = f   Ñ°Ñf+Ñ f   °Ñf.

(2.11)

Substitution of [`(f)] Ñ°Ñf from this equation into (2.10) yields

ó õ   Ñ°( f   Ñf) dA - ó õ   Ñ f   °ÑfdA = -2 ó õ   G q f   dA .

(2.12)

Next we use the divergence theorem which states

ó õ   Ñ°v dA = ó õ C  v °n ds ,

(2.13)

for every vector v and where n is the outward normal to C and ds is the differential element of length along C. Selecting v = [`(f)]Ñf one can convert the first integral over  into an integral over C. This results into the final weak form of the equation as

ó õ C  f   Ñf°n dA - ó õ   Ñ f   °ÑfdA = -2 ó õ   G q f   dA .

(2.14)

2.3   Step 3: Discretization of the domain

The third step is to discretize the domain. Here we will use the three node triangular element which is the simplest two dimensional element to discretize the cross section W. Figure 2.2 shows this element. This three node triangular element has one node at each of its corners.

Figure 2.2 also shows an example of how the domain of the cross section can be modeled using this element. In this example the domain is modeled by six triangular elements. The element number is denoted by a number contained in a circle. This model of the domain is called a finite element mesh. This particular mesh has seven nodes.

The nodes in each element have both a ``local node number'' and a ``global node number''. The global node number refers to the node numbers given when discretize the domain. Each element also has a set of local node numbers. In this case the local node numbers for each element are 1, 2, and 3. For example, Figure 2.3 shows both the local and global node numbers for element 2. As can be seen from the figure the local node 1 refers to the global node 2, the local node 2 refers to the global node 8 and the local node 3 refers to the global node 9.

There is no order in how one numbers the elements and gives the global node numbers. To avoid keeping track of sign changes, in this problem we will always select local node numbers as counter-clock-wise.

It is essential that we have a method for going from the local node numbers to the global node number. To do this we will introduce node connectivity which is a matrix denoted by IJK(I,N), where IJK is the global node number for a node with a local node number I in element N. As an example consider the local and global node numbers given in Figure 2.3. The connectivity for element 2 will be given by

IJK(1,2)=2 ,         IJK(2,2)=8 ,        IJK(3,2)=9 .

(2.15)

The function f is approximated on each element by

f = f1 N1(x,y) +f2 N2(x,y) +f3 N3(x,y) ,

(2.16)

where f₁, f₂, and f₃ are the values of f at local node number 1, 2, and 3, respectively, and N₁, N₂ and N₃ are known functions of x and y. N_i are known as shape functions. Let (x₁,y₁), (x₂,y₂) and (x₃,y₃) denote the coordinates of the three nodes of an element. The three conditions f(x₁,y₁) = f₁, f(x₂,y₂) = f₂ and f(x₃,y₃) = f₃ can be satisfied by nine restrictions on the three N_i. At the first node

N1(x1,y1) = 1,         N2(x1,y1) = 0,         N3(x1,y1)=0 .

(2.17)

At the second node

N1(x2,y2) = 0,         N2(x2,y2) = 1,         N3(x2,y2)=0 .

(2.18)

At the third node

N1(x3,y3) = 0,         N2(x3,y3) = 0,         N3(x3,y3)=1 .

(2.19)

Since there are three restrictions on each shape function, one can represent each shape function by a linear function of x an y and solve for the unknown constants. Let us select for each element e the shape functions

N1e = a1e +b1e x + c1e y (2.20) N2e = a2e +b2e x + c2e y (2.21) N3e = a3e +b3e x + c3e y (2.22)

where the value of the constants for each element can be evaluated from the above restrictions on N_i. To simplify the notation, we will leave out the superscript e from the notation.

Since in addition to the area integrals we have line integrals on the boundary, we need to discretize the boundary ¶W. For the boundary we will use two node line elements. In this case each element has one node at each end. Figure 2.4 shows the discretization of the boundary. As in the discretization of the area, we have both a local and a global numbering system for the nodes. The connectivity of these two numbering systems is given by a matrix IJ(I,N), where IJ is the global node number for local node I of element N.

The function f in the boundary will be approximated by a function of the form

f = f1 ^ N   1  (x,y) +f2 ^ N   2  (x,y) ,

(2.23)

where f₁ and f₂ are the values of f at the two nodes and [^N]₁ and [^N]₂ are the shape functions. Using a local coordinate system as shown in figure 2.5, one can write these shape functions as

^ N   1  = 1- s l ,       and        ^ N   2  = s l

(2.24)

It must be mentioned that the sum of the shape functions at each point of the element must equal unity. This can be used as a check to see if one has selected the right shape functions.

2.4   Step 4: Finite element approximation of the weak form

For a problem with NN number of nodes the goal is to set the problem in the form

[K]{f} = { f} ,

(2.25)

where [K] is an (NN,NN) matrix of constants known as the global stiffness matrix, { f} is the (NN,1) matrix of unknown values of f at the nodes, and { f} is known as the global load vector and is an (NN,1) matrix of known constants. Once the problem is set into this form one can use a linear equation solver to get the unknown f.

Selecting the region  in equation (2.14) as the region of one element W_e, one can get one equation for each element. Adding these equations over all the elements one gets

NBE å e=1  ó õ ¶We  f   Ñf°nds - NE å e=1  ó õ We  Ñ f   °ÑfdA = -2 NE å e=1  ó õ We  Gq f   dA ,

(2.26)

where NE is the number of elements and NBE is the number of boundary elements. Note that the boundary integrals on the internal boundary between any two adjacent elements cancels out of the summation leaving only the boundary integral over ¶W.

On each element we will use the approximations for f and [`(f)] as

f = f1 N1 + f2 N2 + f3 N3 , (2.27) f   = f   1  N1 + f   2  N2 + f   3  N3 , (2.28)

where it is understood that the numbers refer to local node numbers for element e, and that the N_i are the shape functions for element e. This equation can be written as

f = é ë N1 N2 N3 ù û é ê ê ê ë f1 f2 f3 ù ú ú ú û

(2.29)

and

f   = é ê ë f   1  f   2  f   3  ù ú û é ê ê ê ë N1 N2 N3 ù ú ú ú û

(2.30)

The derivatives of f and [`(f)] can be written as

¶f ¶x = é ê ê ê ë ¶N1 ¶x ¶N2 ¶x ¶N3 ¶x ù ú ú ú û é ê ê ê ë f1 f2 f3 ù ú ú ú û ,

(2.31)

¶f ¶y = é ê ê ê ë ¶N1 ¶y ¶N2 ¶y ¶N3 ¶y ù ú ú ú û é ê ê ê ë f1 f2 f3 ù ú ú ú û ,

(2.32)

¶ f   ¶x = é ê ë f   1  f   2  f   3  ù ú û é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë ¶N1 ¶x ¶N2 ¶x ¶N3 ¶x ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú û ,

(2.33)

¶ f   ¶y = é ê ë f   1  f   2  f   3  ù ú û é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë ¶N1 ¶y ¶N2 ¶y ¶N3 ¶y ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú û .

(2.34)

The first equation for the area integral over W_e can be written as

ó õ We  Ñ f   °ÑfdA = ó õ We  ¶ f   ¶x ¶f ¶x + ¶ f   ¶y ¶f ¶y dA .

(2.35)

The first term of this equation can be written as

ó õ We  ¶ f   ¶x ¶f ¶x dA = é ê ë f   1  f   2  f   3  ù ú û ó õ We  é ê ê ê ê ê ê ê ê ê ê ê ê ê ê ë ¶N1 ¶x ¶N2 ¶x ¶N3 ¶x ù ú ú ú ú ú ú ú ú ú ú ú ú ú ú û é ê ê ê ë ¶N1 ¶x ¶N2 ¶x ¶N3 ¶x ù ú ú ú û dA é ê ê ê ë f1 f2 f3 ù ú ú ú û .

(2.36)

A similar expression can be obtained for the second term of this integral. Therefore, this integral can be represented as