Properties of materials - Biomedical Engineering Desk Reference - page 116

The finite element equations

ðK þ qMÞ y ¼ F (3.1.3.6)

that is, a set of simultaneous equations for the nodal
parameters y .

There are four basic methods of formulating the equations
of FEA. These are: (i) the direct or displacement method,
(ii) the variational method, (iii) the weighted residual
method, and (iv) the energy balance method. Only the
more popular variational and weighted residual methods
will be described here. The integral equation 2 will be used
to illustrate the variational method, while the differential
equation system 3.1.3.1a-c will be used to illustrate the
weighted residual method.

Weighted residual approach

The weighted residual approach can be applied directly to
any system of differential equations such as 3.1.3.1a-c and
even to those problems for which a variational principle
may not exist. This approach is therefore more general.
Themethod assumes an approximation u h ( x , y ) for the real
solution u ( x , y ). Because u h is approximate, its substitution
into Eq. 3.1.3.1a will result in an error or residual R h :

The variational approach

The FEM is introduced in the following way. The region
is divided into a finite number of elements of size h
( Fig. 3.1.3-4 ). The h notation is to be interpreted as re-
ferring to the subdivided domain. Instead of seeking the
function y that minimizes I (y) in the continuous domain,
i.e., the exact solution, we instead seek an approximate
solution by looking for the function y h that minimizes
I (y h ) in the discrete domain. The following trial functions
are defined over the discretized domain:
y h ðx; yÞ¼ X

R h ¼V 2 u h þ qu h f (3.1.3.7)

The weighted residual approach requires that some
weighted average of the error due to nonsatisfaction of
the differential equation by the approximate solution u h
(Eq. 3.1.3.7 ) vanish over the domain of interest:
ðð

R h wdU ¼ ðð
U

ðV 2 u h þ qu h fÞwdU ¼ 0

U

n

(3.1.3.8)

y i N i ðx; yÞ (3.1.3.3)

i ¼ 1

where w ( x , y ) is a weighting function. A function u h that
satisfies Eq. 3.1.3.8 for all possible w selected from a cer-
tain class of functions must necessarily satisfy the original
differential equations 3.1.3.1a-c . It actually does so only
in an average or ''weak'' sense. Equation 3.1.3.8 is there-
fore known as a ''weak form'' of the original equation
3.1.3.1a. The second-order derivatives of the V 2 term are
usually reduced to first order derivatives by an integration
by parts (Harris and Stocker, 1998). The result is another
weak form:
ðð

where N i are global basis or shape functions and y i are
nodal parameters. The sum is over the total number of
nodes n in the mesh. Using Eq. 3.1.3.3 in Eq. 3.1.3.2 , the
functional becomes
Ið y h Þ¼ X
i;j

y i y j ðð
U

VN i VN j d U þ q X
i;j

y i y j

ðð y i ðð

U

VN i N j dU 2 X
i; j

fN i d U

U

n Vu h Vw þ qu h w fw o dU ¼ 0

(3.1.3.4)

(3.1.3.9)

U

which can be written in matrix notation as

which has the advantage that approximating functions can
now be chosen from a much larger space, a space where
the function only needs be once-differentiable. Again, we
divide the region into a finite number of elements and
assume that the approximate solution can be represented
by the sum of the product of unknown nodal values y j and
interpolating functions N j ( x , y), defined at each node j of
the mesh:

Ið y h Þ¼ y T K y þ q y T M y 2y T F (3.1.3.5)

where
K ¼ ðð
U

VN i VN j dU; M ¼ ðð
U

N i N j dU;

F ¼ ðð
U

fN i dU

u h ¼ X

n

y j N j (3.1.3.10)

y T represents the transpose of the vector y; K is known as
the stiffness matrix, M as the mass matrix, and F as
the local load vector. The function y h that minimizes
Eq. 3.1.3.5 should satisfy d I (y h ) ¼ 0. This gives

j ¼ 1

When Eq. 3.1.3.10 is substituted into Eq. 3.1.3.9 with
w ¼ N i

, Eq. 3.1.3.6 results as before, proving the