4

Principles for Weak Forms

Principles used for creating weak forms for MFree methods are outlined in a concise,
practical, and easy-to-understand manner in this chapter. We introduce the principles, and
explain the procedures for their use. The principles used in MFree methods are similar to
those used in the conventional finite element method (FEM). However, MFree methods
use these principles in much more innovative ways. Therefore, this chapter covers these
principles in a new systematic way for easy application to MFree methods. Readers with
experience in FEM may skip the chapter, but it may prove helpful to skim the chapter to
become familiar with the terminology conventions used in the book. Readers with less
experience with FEM may have difficulty fully understanding the details of all the pro-
cedures and equations presented here, but a global view of the general procedure and
principles will help in the following chapters. We recommend that such readers study this
chapter without delving too deeply into the equations. Subsequently, it may be a good
idea to review this chapter again after reading the following chapters to understand the
details of all the equations presented in this chapter.

4.1

Strong Forms vs. Weak Forms

The partial differential system equations developed in

Chapter 3

are

strong form

system

equations, such as Equations 3.20 through 3.22 for the mechanics of solids.

Weak form

, in

contrast to strong form, requires weaker consistency on the assumed function of dependent
variables (

u

,

v

,

w

in this case). Obtaining the

exact

solution for a strong form system

equation is ideal but, unfortunately, it is usually very difficult for practical engineering
problems that are very complex. The finite difference method (FDM), which uses finite
differential representation (Taylor series) of a function in a local domain, can be used to
solve system equations of strong form to obtain an

approximated

solution. However, FDM

requires a regular mesh of grids, and can usually work only for problems with regular
geometry and simple boundary conditions. One of the MFree methods for solving strong
form system equations to obtain an approximate solution is to use arbitrarily distributed
grids based on Taylor series expansions, and least square or moving least squares (MLS)
approximations. The formulation is very simple but less stable, and the accuracy of the
results often depends on the selection of the nodes for constructing the differential equa-
tions. There are also MFree methods for solving strong form system equations using an
integral representation of field variable functions, such as the smooth particle hydrody-
namics (SPH) methods. This formulation can deal well with dynamic problems of infinite
domain, such as problems in astrophysics (Monaghan, 1992). SPH has proved stable for
arbitrarily distributed nodes. This is due to the use of integral representation of field func-
tions, which pass the differentiation operations on the field function to the weight function.
Therefore, it reduces the requirement on the order of consistency on the approximated

1238_Frame_C04 Page 53 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

field function, and is actually very similar to that of weak form formulations. The difference
is that the weak form operation is implemented in the stage of shape function construction,
and not in the stage of creating the system equation. The major problem with the SPH
methods is the treatment of boundaries of the problem domain and the boundary conditions.

Formulation based on weak forms can produce a stable set of algebraic system equations

and gives discretized system equations that produce much more accurate results. There-
fore, weak form is often preferred by many to obtain an approximated solution. This book
uses weak form formulations to form discretized system equations for our MFree methods
for mechanics problems of solids, structures, as well as fluids.

The consistency requirement on the approximation functions for field variables is very

different for strong form and weak form formulation. For a 2

k

th-order differential gov-

erning system equation, the strong form formulation assumes the field variable possesses
a consistency of 2

k

th order. The weak form formulation, however, requires a consistency

of only

k

th order.

There are basically two major categories of often-used principles for constructing weak

forms. They are variational methods and weighted residual methods. There are different
forms of variational methods. Galerkin formulation may be the most widely used
approach for establishing system equations. It is, of course, applicable to deriving MFree
equations. Hamilton’s principle is often employed to produce approximated system equa-
tions for dynamic problems, and is also applicable to MFree methods. The minimum total
potential energy principle has been a convenient tool for deriving discrete system equa-
tions for FEM and also for many other types of approximation methods. The weighted
residual method is a more general and powerful mathematical tool that can be used for
creating discretized system equations for many types of engineering problems. It has been
and is still used for developing new MFree methods. All these approaches are adopted in
this book for creating discretized system equations for various types of MFree methods.

4.2

Hamilton’s Principle

Hamilton’s principle is one of the variational principles based on the energy principle; it
states, “Of all possible time histories of consistent displacement states which satisfy

∗

(a) the compatibility conditions

(4.1)

(b) the essential (displacement or kinematical) boundary conditions

(4.2)

(c) the conditions at initial (

t

1

) and final time (

t

2

)

(4.3)

the history corresponding to the actual solution makes the

Lagrangian functional

a mini-

mum.” These conditions are also called admissible conditions. Mathematically, Hamilton’s
principle states

(4.4)

where

L

is the Lagrangian functional. For a system of solids and structures it can be defined as

L

=

T

−

Π

s

+

W

f

(4.5)

∗

See Section 5.1 for the definitions of consistency and compatibility.

δ L dt

t

1

t

2

∫

0

=

1238_Frame_C04 Page 54 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

where

T

is the kinetic energy,

Π

is the strain energy, and the

W

f

is the work done by the

external forces.

The kinetic energy is defined by

(4.6)

where

Ω

stands for the whole volume of the solid. For solids and structures of elastic

materials, the strain energy of the system can be expressed as

(4.7)

The work done by the external forces can be obtained by

(4.8)

where

Γ

t

stands for the boundary of the solids on which traction forces are prescribed

(see

Figure 3.2).

The procedure for solving a problem using Hamilton’s principle is described as follows:

1. Construct a shape function to approximate the field functions using their values

at nodes in the problem domain. The approximated field function should be
consistent and admissible, i.e., it satisfies the conditions of Equations 4.1 to 4.3.

2. Calculate the kinetic energy

T

, the strain energy

Π

, and the work done by the

external forces

W

f

in terms of the approximated field functions.

3. Form the Lagrangian functional

L

using Equation 4.5, and then substitute it into

Equation 4.4. Use integration by parts and carry out the spatial and temporal
integrations that lead to a set of differential equations with respect only to time.

4. Solve the set of differential equations with respect to time, using standard pro-

cedures to obtain the dynamic field.

5. For static problems, drop the time-related terms, and step 3 leads to a set of

algebraic equations, which can be solved using a standard algebraic equation
solver to obtain the static field.

At step 1, the shape function is constructed in FEM using a mesh of elements. In MFree
methods, techniques without using elements are required.

Note that in order to use Hamilton’s principle, a functional has to be obtained, which

may not be possible for all problems. The advantage of using Hamilton’s principle is that
we do not have to know the strong form of the system equation.

4.3

Constrained Hamilton’s Principle

There are cases when the approximated field function does not satisfy the condition,
Equations 4.1 and 4.2, on parts of the problem domain including parts of the boundaries,
discrete curves, and points of locations. Hamilton’s principle has to be modified or con-
strained for such situations.

T

1
2

---

ρu˙

T

u˙

d

Ω

Ω

∫

=

Π

s

1
2

---

εεεε

T

σσσσ dΩ

Ω

∫

=

W

f

u

T

b

d

Ω

Ω

∫

u

T

t dΓ

Γ

t

∫

+

=

t

1238_Frame_C04 Page 55 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

Consider the following set of

k

conditions that the approximated field function cannot

satisfy:

(4.9)

where

C

is a given matrix of coefficients.

Our purpose now is to seek the stationary point of the Lagrangian functional subjected

to the constraint of Equation 4.9. There are basically two methods often used to modify
the functional that accommodates these constraints. They are the method of Lagrange
multipliers and the penalty method.

4.3.1

Method of Lagrange Multipliers

In the method of Lagrange multipliers, the modified Lagrangian is written as follows:

(4.10)

where

λλλλ

is a vector of the Lagrange multipliers given by

(4.11)

These Lagrange multipliers are unknown functions of independent coordinates in the
domain

Ω

. The modified Hamilton’s principle then seeks the stationary of the following

functional:

(4.12)

Note that since the Lagrange multipliers are unknown functions, the total number of
unknown field functions of the whole system is then increased. In the process of seeking
discretized system equations, these Lagrange multipliers must also be approximated in a
manner similar to the field functions. Therefore, the total number of

nodal unknowns

in

the discretized system equation will also be increased. The method of Lagrange multipliers
will, however, rigorously enforce the constraints. The penalty method, introduced in the
following subsection, does not increase the number of unknowns.

4.3.2

Penalty Method

In examining the constraint equations, Equation 4.9, we construct the following functional:

(4.13)

C u

( )

C

1

u

( )

C

2

u

( )
M

C

k

u

( )





























0

=

=

L˜

L

λλλλ

T

C u

( )dΩ

Ω

∫

+

=

λλλλ

T

λ

1

λ

2

…

λ

k

{

}

=

δ L˜ dt

t

1

t

2

∫

0

=

C

T

ααα

αC

α

1

C

1

2

α

2

C

2

2

… α

k

C

k

2

+

+

+

=

1238_Frame_C04 Page 56 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

where

ααα

α

is a diagonal matrix given by

(4.14)

where

α

1

,

α

2

,…,

α

k

are penalty factors. They can be given as functions of coordinates, but

usually they are assigned positive constant numbers. In any case,

C

T

ααα

α

C

will always be

non-negative, and hence zero is the minimum of the functional

C

T

ααα

α

C

. It would be zero

only if all the conditions in Equation 4.9 are fully satisfied. Therefore, the following
stationary condition of the functional

C

T

ααα

α

C

guarantees the best satisfaction of the con-

straint equations, Equation 4.9:

(4.15)

Performing the variation using the chain rule, we have

(4.16)

which leads to the following minimization condition:

(4.17)

If

α

i

=

0, the essential boundary condition is not enforced at all, because any C

i

will still

satisfy the ith equation in Equation 4.17. If

α

i

goes to infinity, the essential boundary

condition is fully enforced, because C

i

must be zero in order to satisfy Equation 4.17.

The above analysis shows that C

T

ααα

αC should be the right functional for our constraint.

The modified Lagrangian is then written as follows:

(4.18)

The serves only to counter the 2 that will be produced in the later variational operation. The
important difference between the penalty factor

α

i

and the Lagrange multiplier

λ is that

the penalty factor is a given constant (no variation is allowed), whereas the Lagrange
multiplier is a variable.

Because

α is a known constant, there will be no increase in the number of unknowns

in the system. The question is how to choose the penalty factor. To impose the constraint
fully, the penalty factor must be infinite, which is not possible in practical numerical
analysis. Therefore, in the penalty method these constraints cannot be satisfied exactly,
but only approximately. In general, the use of a larger penalty factor will enforce the
constraint. The problem is that if the penalty factor is too small, the constraints will not

ααα

α

α

1

0

0

0

0

α

2

0

0

0

0

O 0

0

0

0

α

k

=

δ C

T

ααα

αC

(

)

0

=

δ C

T

ααα

αC

(

)

2C

T

ααα

αδC

2

δC

T

ααα

αC

0

=

=

=

ααα

αC u

( )

α

1

C

1

u

( )

α

2

C

2

u

( )

M

α

k

C

k

u

( )





























0

=

=

L˜

L

1
2

---

C

T

u

( )ααα

αC u

( )dΩ

Ω

∫

+

=

1
2

--

1238_Frame_C04 Page 57 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

be properly enforced, but if it is too large, numerical problems will be encountered. A
compromise should be reached. Some kind of formula that is universally applicable would
be useful. To find such a formula, one may need to determine the factors that affect the
selection of the penalty factor in the actual event of solving the discretized system equations.

Note that use of the penalty method is a routine operation even in the FEM for imposing

essential boundary conditions including single-point and multi-point constraints.

4.3.3

Determination of Penalty Factor

Because discretization errors can be comparable in magnitude to the errors due to the
poor satisfaction of the constraint, Zienkiewicz (1989) has suggested using the following
formula for FEM analysis:

α = constant(1/h)

n

(4.19)

where h is the characteristic length, which can be the ratio of the element size to the
dimension of the problem domain, and n is the order of the elements.

In extending this formula to MFree methods, we suggest that h be the ratio of the nodal

spacing to the dimension of the problem domain, and n

= 1. The constant in Equation 4.19

should relate to the material property of the solid or structure. It can be 10

10

times Young’s

modulus.

This book prefers the following simple method for determining the penalty factor:

α = 1.0 × 10

4–13

× max(diagonal elements in the stiffness matrix)

(4.20)

In most of the examples reported in later chapters using penalty methods, the foregoing
equation is adopted.

It has also been suggested to use

α = 1.0 × 10

5

∼8

× Young’s modulus

(4.21)

for some examples, which works well.

Note that trials may be needed to choose a proper penalty factor. More discussion on

this is provided in Example 6.12 in Section 6.4.5, when we deal with nonlinear problems.

4.4

Galerkin Weak Form

The Galerkin weak form can be directly derived using Hamilton’s principle for problems
of solid mechanics. By using Equations 4.5 to 4.8, the Lagrangian L becomes

(4.22)

Substituting the above into Equation 4.4, we have

(4.23)

L

1
2

---

εεεε

T

σσσσ dΩ

u

T

b

d

Ω

u

T

t

d

Γ 1

2

---

ρu˙

T

u˙

d

Ω

Ω

∫

+

Γ

t

∫

+

Ω

∫

+

Ω

∫

–

=

δ

1
2

---

εεεε

T

σσσσ dΩ

u

T

b

d

Ω

u

T

t

d

Γ 1

2

---

ρu˙

T

u˙

d

Ω

Ω

∫

+

Γ

t

∫

+

Ω

∫

+

Ω

∫

–

dt

t

1

t

2

∫

0

=

1238_Frame_C04 Page 58 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

Moving the variation operation into the integral operations, we obtain

(4.24)

Changing the order of operation does not affect the results because all these operations
are linear, and therefore it does not matter which one operates first. The integrand in the
first integral term can be written as follows using the chain rule of variation:

(4.25)

Because the two terms in the foregoing equation are all scalars, their transposes are still
themselves. For the last term in Equation 4.25, we should have

(4.26)

Using the constitutive equation of Equation 3.8 and the symmetric property of the matrix
of material constant c, we should have

(4.27)

Therefore, Equation 4.25 becomes

(4.28)

Let’s now look at the last term in Equation 4.24. We first move the time integration into

the spatial integration; we should have

(4.29)

Again, this change of the order of integral operation does not affect the results because
they are linear operations. The time integration in Equation 4.29 can be changed as follows,
using the chain rule of variation and then the scalar property:

(4.30)

Note that the variation and time differentiation are also exchangeable, i.e.,

(4.31)

Integrating by parts with respect to the time, we have

(4.32)

1
2

---

δ εεεε

T

σσσσ

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ 1

2

---

δ ρu˙

T

u˙

(

)dΩ

Ω

∫

+

Γ

t

∫

+

Ω

∫

+

Ω

∫

–

dt

t

1

t

2

∫

0

=

δ εεεε

T

σσσσ

(

)

δ εεεε

T

σσσσ εεεε

T

δ σσσσ

+

=

εεεε

T

δ σσσσ

εεεε

T

δ σσσσ

(

)

T

δ σσσσ

T

εεεε

=

=

δ σσσσ

T

εεεε

δ cεεεε

( )

T

εεεε

δ εεεε

T

c

T

εεεε

δ εεεε

T

c

εεεε

δ εεεε

T

σσσσ

=

=

=

=

δ εεεε

T

σσσσ

(

)

2

δ εεεε

T

σσσσ

=

1
2

---

δ ρu˙

T

u˙

(

)dΩ

Ω

∫

dt

t

1

t

2

∫

1
2

---

δ ρu˙

T

u˙

(

)dt

t

1

t

2

∫

Ω

∫

d

Ω

=

δ ρu˙

T

u˙

(

)dt

t

1

t

2

∫

ρ

δ u˙

T

u˙

u˙

T

δ u˙

+

[

]dt

t

1

t

2

∫

2

ρ

δ u˙

T

u˙

[

]dt

t

1

t

2

∫

=

=

δ u˙

T

u˙

[

]dt

t

1

t

2

∫

d

δ u

T

dt

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

du

dt

-------

dt

t

1

t

2

∫

=

d

δ u

T

dt

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

du

dt

-------

dt

t

1

t

2

∫

δ u

T

d

2

u

dt

2

---------

–

dt

δ u

T

du

dt

-------

t

1

t

2

+

t

1

t

2

∫

=

1238_Frame_C04 Page 59 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

Invoking now the condition given in Equation 4.3, we know that the u has already satisfied
the conditions at the initial (t

1

) and the final time (t

2

). Therefore,

δu

T

has to be zero at t

1

and t

2

(no variation can exist for any given constant value). Therefore, the last term in

Equation 4.32 vanishes, which gives

(4.33)

Therefore, Equation 4.29 becomes

(4.34)

We now switch the order of integration to obtain

(4.35)

Equation 4.24 now becomes

(4.36)

To satisfy the above equation for all possible choices of u, the integrand of the time
integration has to vanish, which leads to

(4.37)

This is the well-known Galerkin weak form. For static problems, it simply reduces to

(4.38)

Equation 4.37 can also be viewed as the principle of virtual work, which states that if a

solid body is in equilibrium, the total virtual work performed by all the stresses in the
body and all the external forces applied on the body should vanish when the body is
subjected to a virtual displacement. The virtual work can be viewed as an alternative
statement of the equilibrium equation. In the situation given in Equation 4.37, we can view
that the solid is subjected to a virtual displacement of du. The first term in Equation 4.37
is the virtual work done by the internal stress in problem domain

Ω; the second term is

the virtual work done by the external body force; the third term is the virtual work done
by the external tractions on the natural boundaries; and the last term is the virtual work
done by the inertial forces. Therefore, using the principle of virtual work, we can actually
write out Equation 4.37 or 4.38 straight away.

d

δ u

T

dt

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

du

dt

-------

dt

t

1

t

2

∫

δ u

T

d

2

u

dt

2

---------

–

dt

t

1

t

2

∫

=

1
2

---

δ ρu˙

T

u˙

(

)dΩ

Ω

∫

dt

t

1

t

2

∫

ρ

δ u

T

d

2

u

dt

2

---------

–

dt

t

1

t

2

∫

Ω

∫

d

Ω

=

1
2

---

δ ρu˙

T

u˙

(

)dΩ

Ω

∫

dt

t

1

t

2

∫

ρ

δ u

T

u

˙˙

[

]dΩ

Ω

∫

dt

t

1

t

2

∫

–

=

δ εεεε

T

σσσσdΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

ρδ u

T

u

˙˙d

Ω

Ω

∫

–

Γ

t

∫

+

Ω

∫

+

Ω

∫

–

dt

t

1

t

2

∫

0

=

δ εεεε

T

σσσσ dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

ρδ u

T

u

˙˙d

Ω

Ω

∫

+

Γ

t

∫

–

Ω

∫

–

Ω

∫

0

=

δ εεεε

T

σσσσ dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

Γ

t

∫

–

Ω

∫

–

Ω

∫

0

=

1238_Frame_C04 Page 60 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

The Galerkin weak form can be formulated in many other ways. One alternative is to

formulate it via the method of weighted residuals, which is introduced in the following
section.

By using the stress–strain relation of Equation 3.8, and then the strain–displacement

relation of Equation 3.5, Equations 4.5 and 4.38 can be explicitly expressed as follows in
terms of displacement vector u:

(4.39)

This is the Galerkin weak form written in terms of displacement, and therefore it is
convenient to use, as the displacement is to be approximated in FEM or MFree methods.
For static problems, Equation 4.39 simply reduces to

(4.40)

The above two equations of Galerkin weak form are very handy in application to

problems of solid mechanics, because one does not need to perform integration by parts.
The discretized system equation can be derived very easily using approximated displace-
ments. The Galerkin procedure used in MFree methods is as follows:

1. Approximate the displacement at a point using MFree shape functions and nodal

displacements at the nodes surrounding the point. The approximation should
be consistent and has to satisfy Equations 4.1 and 4.2.

2. Substitute the approximated displacements into Equation 4.37, and factor out

the variation of the nodal displacements, leading to a set of differential equations
with respect only to time.

3. Solve the set of differential equations with respect to time, using standard pro-

cedures to obtain the dynamic field.

4. For static problems, use Equation 4.38 instead, which leads to a set of algebraic

equations that can be solved using standard algebraic equation solvers to obtain
the static field.

This Galerkin procedure is applied repeatedly in the following chapters for all kinds of
mechanics problems of solids and structures.

Note that in using the Galerkin procedure one need not know the strong form of the

system equation.

4.5

Constrained Galerkin Weak Form

For cases when the approximated field function does not satisfy the condition, Equation 4.1
or 4.2, on parts of the problem domain including parts of the boundaries and discrete points
of locations, we should use the constrained Hamilton’s principle to derive the constrained
Galerkin weak form. The procedure is the same as in Section 4.4, except the modified
Lagrange is used for the formulation. The following presents the final expressions.

δ Lu

(

)

T

c Lu

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

ρδ u

T

u

˙˙d

Ω

Ω

∫

+

Γ

t

∫

–

Ω

∫

–

Ω

∫

0

=

δ Lu

(

)

T

c Lu

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

Γ

t

∫

–

Ω

∫

–

Ω

∫

0

=

L˜

1238_Frame_C04 Page 61 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

4.5.1

Galerkin Weak Form with Lagrange Multipliers

For dynamic problems,

(4.41)

For static problems, it simply reduces to

(4.42)

4.5.2

Galerkin Weak Form with Penalty Factors

For dynamic problems,

(4.43)

For static problems, it simply reduces to

(4.44)

4.6

Minimum Total Potential Energy Principle

Another very often used energy principle in FEM is the minimum total potential energy
principle. This principle states that for a structure system that is at equilibrium, the total
potential energy in the system must be stationary for variations of admissible displace-
ments. This principle works in a very straightforward manner in the following three simple
steps:

1. Approximate the field function (displacement) in terms of the nodal variables,

say, d, which is a vector that collects all the nodal displacements in the problem
domain.

2. Express the total potential energy in terms of the nodal variables d. For solids

and structures of elastic materials, the total potential energy can be expressed as

Π = Π

s

− W

f

(4.45)

δ Lu

(

)

T

c Lu

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

Γ

t

∫

–

Ω

∫

–

Ω

∫

δ λλλλ

T

C u

( )dΩ

λλλλ

T

δ C u

( )dΩ

Ω

∫

–

Ω

∫

–

ρδ u

T

u

˙˙d

Ω

Ω

∫

+

0

=

δ Lu

(

)

T

c Lu

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

δ λλλλ

T

C u

( )dΩ

λλλλ

T

δC u

( )dΩ

Ω

∫

–

Ω

∫

–

0

=

Γ

t

∫

–

Ω

∫

–

Ω

∫

δ Lu

(

)

T

c Lu

(

)dΩ

δ u

T

b

d

Ω

δ u

T

t

d

Γ

δ C u

( )

T

ααα

αC u

( )dΩ

ρδ u

T

u

˙˙d

Ω

Ω

∫

+

Ω

∫

–

0

=

Γ

t

∫

–

Ω

∫

–

Ω

∫

δ Lu

(

)

T

c Lu

(

)dΩ

δu

T

b

d

Ω

δu

T

t

d

Γ

δ C u

( )

T

ααα

αC u

( )dΩ

Ω

∫

–

0

=

Γ

t

∫

–

Ω

∫

–

Ω

∫

1238_Frame_C04 Page 62 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

3. Use the following stationary conditions to create a set of discretized system

equations:

(4.46)

The number of the equation created is exactly the total number of the nodal
variables.

4.7

Weighted Residual Method

The weighted residual method is a very simple but powerful mathematical tool to obtain
approximated system equations of weak form. The concept of the weighted residual
method is straightforward and applicable, in principle, for most of the partial differential
equations that govern engineering problems, including mechanics of solids, structures,
and fluids.

Let us consider the general form of the partial differential equations for solid mechanics.

It can be rewritten in a general concise functional form:

f

(u(x, y, z))

= 0

(4.47)

For example, for solid mechanics problems formulated in Section 3.2.3, we have

(4.48)

In general, it is difficult to obtain the exact solution u(x, y, z) that satisfies Equation 4.47.
We therefore seek a u(x, y, z) that satisfies Equation 4.47 in a weighted integral sense over
a quadrature domain:

(4.49)

where

is a vector of weight functions corresponding to the equations in f. We hope

that u(x, y, z) is a good approximation of the exact solution. This is the idea of the weighted
residual approach. It is, indeed, very simple.

Although the method is simple and, in principle, it is applicable to most partial differ-

ential equations, different ways of implementation will lead to solutions of different
accuracies. Variational methods and finite difference methods can all be derived using the
weighted residual method by properly choosing the weight functions.

Note that function f can be either static or dynamic. The weighted residual procedure

will produce these two types of discrete systems accordingly.

∂Π

∂d

-------

∂Π

∂d

1

--------

∂Π

∂d

2

--------

M





























0

=

=

f u

x, y, z

(

)

(

)

L

T

cLu

b

ρu˙˙

–

+

=

W

T

f u

x, y, z

(

)

(

)dΩ

Ω

∫

0

=

)

W

)

1238_Frame_C04 Page 63 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

The procedure for solving a problem using the weighted residual method is as follows:

1. Construct a shape function to approximate the field function using the field

variables at the nodes in the domain with a certain order of consistency.

2. Construct weight functions. When the Dirac delta function is used, the weighted

residual method leads to the collocation method. If the shape function con-
structed in step 1 is also used as the weight function, the weighted residual
method leads to the Galerkin weak form, which was derived in the previous
section. In cases where the weight function constructed is different from the
shape function used, the method is often termed the Petrov–Galerkin method,
which is used in this book frequently.

3. Substitute the shape and weight functions into Equation 4.49 using integration

by parts, and carry out the integrations using boundary conditions, which leads
to a set of differential equations with respect only to time.

4. Solve the set of differential equations using standard procedures to obtain the

dynamic field.

5. For static problems, step 3 leads to a set of algebraic equations, which can be

solved using standard algebraic equation solvers to obtain the static field.

In the strong form used in Equation 4.49, there may be higher-order differentiations, say,

an order of 2k. However, the approximated field function needs a consistency of only kth
order, because the operation of integration by parts transfers another kth order of differen-
tiation to the weight function. This reduction of the requirement of order of consistency on
the approximated field function is very important to create equation systems that provide
accurate solutions. On the other hand, due to the translation of the differentiation to the
weight function, the weight function used has to be differentiable up to the kth order. Other
than this, there is no compulsory requirement on the weight function used. We are free to
use weight functions that give us convenience in the formulation. The choice of weight
function will also affect to a certain degree the accuracy of the results.

Note that to use the weighted residual method, the strong form of the system equations

needs to be known, but we do not have to know whether there is a functional available
for the problem.

It is very important to mention here that the quadrature domain

Ω used in Equation 4.49

does not have to be the entire domain of the problem. When a local quadrature domain
is used, the weighted residual method can be termed a local weighted residual method. This
local weighted residual method is used very often in later chapters, in the form of the
local Petrov–Galerkin method where the test and trial functions are chosen independently.

Note also that if the approximated field function contains the exact solution, the

weighted residual method will produce the exact solutions as long as there is no numerical
error in the computation. This important feature is useful in testing MFree methods
developed based on the weighted residual method.

4.8

Weighted Residual Method with Constraints

When we must satisfy constraint equations, such as Equation 4.9, in addition to the strong
form of governing system equations, such as Equation 4.47, the weighted residual method
can be written by simply adding the weighted residual of the penalized constraint equations:

(4.50)

W

T

f u

x, y, z

(

)

(

)dΩ

W

T

ααα

α C u

( )dΩ

Ω

∫

+

Ω

∫

0

=

)

)

1238_Frame_C04 Page 64 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

where

ααα

α is a diagonal matrix given by Equation 4.14. Note that

in the first integral is

a vector of weight functions corresponding to the equations in f, and

in the second

integral is a vector of weight functions corresponding to the equations in C. All these
weight functions can be, but do not have to be, the same.

4.9

Points to Note

NOTE 1

In both the method of Lagrange multipliers and the penalty method, the addi-

tional functionals need not be the domain type. If the unsatisfied conditions are on a
boundary or curves in the domain, the integrals should be changed to curve integrals; if
they are at points, the integrals should be changed to a summation over those points.

NOTE 2

In using the method of Lagrange multipliers, we form the constrained

Lagrangians by adding a functional for the constraint to the original Lagrangians, as
shown in Equation 4.10. Variation leads to a set of simultaneous system equations that is
enlarged with Lagrange multipliers as additional unknowns. This procedure works for
static problems. For dynamic problems, however, one may need special ways to deal with
the enlarged set of dynamic system equations, as both the stiffness and mass matrices are
nonpositive. A more practical way to obtain a set of well-behaved system equations is to
formulate separately the system equations with the original Lagrangian and the constraint
equations using their own weak form to obtain two separate discretized sets of equations.
The discretized constrained equations are decomposed to obtain a set of orthogonal vectors
that are then used to produce a condensed set of system equations using orthogonal
transformation techniques, which at the same time ensures the satisfaction of the con-
straints. We demonstrate this procedure in Section 11.1 in detail when we deal with free-
vibration problems using the EFG method for plates.

4.10 Remarks

We make the following two remarks:

• The weak forms to be used in MFree methods are the same as those used in

FEM. In FEM, one seldom uses the constrained principles and weak forms. In
MFree methods, there are many methods that use the constrained weak forms.

• The procedures used in applying these weak forms in MFree methods will be

slightly different from those in FEM, because of the difference in the forms of
the shape functions. The integration domain may no longer be the union of the
element, and it may overlap depending on the MFree method used.

Because of the uniqueness of the MFree method in creating shape functions, the devel-

opment of efficient methods to create MFree shape functions is the most important issue in
MFree methods. Once a robust and efficient method is developed to approximate admis-
sible field functions, a weak form presented in this chapter can then be used to derive the

W

)

W

)

1238_Frame_C04 Page 65 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC

discretized system equations. The next chapter is therefore devoted entirely to methods
for creating MFree shape functions, and the properties of the MFree shape functions created.

In following chapters, we use Hamilton’s principle, the Galerkin weak form, the mini-

mum total potential energy principle, and the local Petrov–Galerkin method to formulate
various MFree methods for mechanics problems of solids and structures. The local
Petrov–Galerkin method is used to formulate various MFree methods for mechanics
problems of fluids.

1238_Frame_C04 Page 66 Wednesday, June 12, 2002 4:48 PM

© 2003 by CRC Press LLC