Czichos H., Saito T., Smith L.E. (Eds.) Handbook of Metrology and Testing

Подождите немного. Документ загружается.

1028 Part E Modeling and Simulation Methods

gradient then takes the form

L =

∗

∗−1

∗

p−1

∗−1

. (18.61)

The slip system can be addressed using inﬁnitesimal de-

formation theory. Assuming that plastic deformation is

only caused by the sum of primary slips, the evolution

of the plastic part is composed of

p−1



(α)

, (18.62)

where γ

stands for the amount of shear strain of the

slip system (α). The plastic strain rate is then given by

the symmetric part of the above equation. Elastic strain

covers the rest of the strain rate, and the elastic rate

is related to the stress rate via elasticity moduli. The

resolved shear stress acting on the system (α) can be

expressed using the stress theorem as

(α)



(α)



. (18.63)

It is then necessary to connect the shear strain rate and

the resolved shear stress. A power law formula involv-

ing the resolved shear stress and strain rate was pro-

posed by Hutchinson [18.52], Pan and Rice [18.53]as

(α)



(α)



1/m−1



(α)



. (18.64)

Here the function g(γ

) is the resistance, a kind of drag

stress, under constant shear strain rate. The evolution of

g(γ

) is governed not only by the shear slips themselves

but also by the interactions of different slip systems.

This implies that the rate of g(γ

)isexpressedby

(α)



αβ



(β)



. (18.65)

–0.05 0.04 0.13

Strain ε (%)

0.23 0.32 0.41

Fig. 18.14 Strain concentration behavior with different grain boundary geometries

The coefﬁcient h

αβ

addresses the effect of slip (β)on

another slip (α). This coefﬁcient h

αβ

is suitably deﬁned

using a function involving total shear strain, critical

shear stress, interaction parameters and so on. From

these considerations, the total strain rate at the macro-

scopic level is related to the microscopic deformation as

follows: (1) The stress is decomposed into the resolved

shear stress of a particular slip system, (2) the local

stress then acts on the shear slip, and ﬁnally (3) each

slip is summed, giving the macroscopic strain. Solving

the above equation for the stress rate, we obtain

= E

ijkl

−



(α)

ijkl

(α)

(18.66)

using the inﬁnitesimal deformation theory. In the above

equation

(α)

stands for the symmetric part of

(18.62) and therefore the plastic strain rate under in-

ﬁnitesimal deformation. The effect of the crystal grain

size has also been considered, taking the strain gradient

into account [18.54].

Bottom

Top

t (t

)

δ (δ

,δ

)

Fig. 18.13

Grain bound-

ary sliding using

a contact/slide

element

Part E 18.4

Continuum Constitutive Modeling 18.4 Crystal Plasticity 1029

18.4.2 Grain Boundary Sliding

Most metallic materials are composed of large numbers

of crystal grains, each of which is located in a partic-

ular neighborhood deﬁned via grain boundaries. The

orientation of each crystal is believed to be random; the

macroscopic deformation behavior of polycrystalline

materials is reﬂected in this randomness. The material

isotropy of polycrystalline materials is caused by the

competing/misﬁtting effects of crystal slip on neigh-

borhoods. Since the grain boundary can be regarded as

a discontinuous plane under strain, it is particularly im-

portant to ensure that the mechanical analysis and the

boundary sliding model successfully describes not only

polycrystal deformation but also damage evolution at

elevated temperatures [18.55].

One simple way to idealize grain boundary slid-

ing is to use a special contact/slip model for the

grain boundary. A stiffness is introduced as shown

in Fig. 18.13, where the contact model is readily ap-

plied to the contact/slide element in the ﬁnite element

technique [18.56]. The local coordinate system can be

deﬁned at the grain boundary such that the connecting

force is decomposed into the normal force (t

) and shear

force (t

, t

). In the case of linear stiffness, we have

⎛

⎜

⎝

⎞

⎟

⎠

⎡

⎢

⎣

0 k

00k

⎤

⎥

⎦

⎛

⎜

⎝

⎞

⎟

⎠

. (18.67)

The normal component is devoted to the opening and

closing (δ

) of the grain boundary, while the shear

component is related to relative sliding (δ

, δ

). The

normal stiffness k

is set to a large enough value that

it causes neither detaching nor overlapping with the

neighboring crystal. Permanent shear sliding is also al-

lowed and modeling is also performed on the basis of

the Mohr–Coulomb yield condition. The constitutive

model for the contact/sliding is then similar to the model

for the bulk crystal. It should be pointed out, however,

that stiffness parameters are rarely speciﬁed in conven-

tional test data, and it is even believed that large vari-

ations in stiffness occur which affect the real evolution

of damage. Tvergaard and coworkers [18.55]employed

viscous sliding at the grain boundary which they applied

to investigate tertiary creep and damage evolution.

18.4.3 Inhomogeneous Deformation

In this section two examples of simulation are de-

scribed: one considers the role of the geometry of the

grain boundary while the other investigates the devel-

0.47 0.48 0.49

Strain (%)

0.51 0.52

ya) 1×1×1

b) 2×2×2

c) 4×4×4

0.53

Fig. 18.15a–c Inhomogeneous deformation and surface roughness.

The formation of a wavy pattern is observed in a large-scale simu-

lation

opment of surface roughness. In order to evaluate the

effects of the grain boundary, we look at (a) combi-

nations of slip systems, (b) the contact geometry, and

stiffness (c) and the interactions of slip systems (a)

have relatively small effects on the variation in strain

concentration. When the stiffness is large, nodal dis-

crepancy is restrained and so the strain concentration is

enhanced to some extent, but this is not a signiﬁcant fac-

tor compared to the effect of geometry (b). Figure 18.14

Part E 18.4

1030 Part E Modeling and Simulation Methods

shows how the strain distribution relates to the inclined

geometry of the grain boundary [18.57]. The strain dis-

tribution varies greatly in models of highly inclined

crystals. Since the effect of a slip system is examined

in the central ﬁgure, where the grain boundary forms

a horizontally ﬂat surface, the differences observed in

the other ﬁgures reﬂect the effect of grain boundary

geometry. It follows from these ﬁgures that the strain

distribution around the grain boundary is strongly af-

fected by its geometry. Although the inﬂuence of the

boundary geometry is usually far less than the factors

related to mechanical behavior, the randomness of the

geometry dominates the variation in the strain, and the

variation in the strain may be related to damage evolu-

tion in polycrystalline materials.

The size of the crystal grain affects statistical prop-

erties related to mechanical behavior. Here we examine

the effect of size on the surface roughness. It is worth

noting that the real size is not taken into account be-

cause the crystal model does not reﬂect the real size.

Introduction of higher order variables such as strain

gradient or nonlocal variables covering a certain range

of physical distance is a promissing way to capture

the real size effect in crystal grains. Figure 18.15 de-

picts the deformed shape of a ﬁnite element mesh

under monotonic tension along the z-direction [18.57].

Grain boundary sliding is not introduced in the model.

8×8×24= 1536 cubic elements are located in the ob-

jective domain, and each element has random crystal

slip system in Fig. 18.15a.2×2×2(=8) elements com-

pose a block corresponding to a crystal grain, associated

with the same slip system used in Fig. 18.15b, and

4×4×4 (=64) elements form a block with the same

slip system used in Fig. 18.15c. The displacement on the

surface is magniﬁed to emphasize the surface state. The

variation in strain occurs due to competing/misﬁtting ef-

fects with a neighboring grain (block); in other words

the strain concentration increases when the slip sys-

tem has a large misﬁt with the neighborhood. It should

be pointed out that a wavy pattern is predicted by ev-

ery simulation model, and that the wavelength seems

to be proportional to the block (grain) size. In fact,

experimental results suggest that surface roughness is

enhanced during plastic deformation, and that the wave-

length of this is equal to the several times the average

crystal grain size. The simulated results qualitatively

agree with this tendency, although other kinds of un-

certainties, in for example mechanical parameters, grain

boundary stiffness and variation in grain sizes, apply to

real materials.

References

18.1 L.E. Malvern: Introduction to the Mechanics of

a Continuous Medium (Prentice-Hall, Englewood

1969)

18.2 P. Perzyna: Thermodynamic theory of viscoplastic-

ity, Adv. Appl. Mech. 11, 313–354 (1971)

18.3 A.K. Miller: Uniﬁed Constitutive Equations for Creep

and Plasticity (Elsevier, London 1987)

18.4 A.S. Krausz, K. Krausz: Uniﬁed Constitutive Laws

of Plastic Deformation (Academic, San Diego

1996)

18.5 J. Lemaitre, J.L. Chaboche: Mechanics of Solid Ma-

terials (Cambridge Univ. Press, Cambridge 1994)

18.6 G.A. Maugin: Thermomechanics of Plasticity and

Fracture (Cambridge Univ. Press, Cambridge 1995)

18.7 A.K. Miller: An inelastic constitutive model

for monotonic, cyclic, and creep deformation,

Trans. ASME J. Eng. Mater. Technol. 98,97–105

(1976)

18.8 M.C.M. Liu, E. Krempl: A uniaxial viscoplastic model

based on total strain and overstress, J. Mech. Phys.

Solids 27, 377–391 (1979)

18.9 S.R. Bodner, A. Merzer: Viscoplastic constitutive

equations for copper with strain rate history and

temperature effects, Trans. ASME J. Appl. Mech.

100, 388–394 (1978)

18.10 Y. Estrin, H. Mecking: An extension of the

Bodner-Partom model of plastic deformation, Int.

J. Plasticity 1,73–85(1985)

18.11 K.C. Valanis: A theory of viscoplasticity without

a yield surface, Arch. Mech. 23, 517–551 (1971)

18.12 D. Kujawski, Z. Mroz: A viscoplastic material model

and its application to cyclic loading, Acta Mech. 36,

213–230 (1980)

18.13 Y.F. Dafalias, E.P. Popov: A model for nonlin-

ear hardening materials for complex loading, Acta

Mech. 21, 173–192 (1975)

18.14 N. Ohno: A constitutive model for cyclic plasticity

with a nonlinear strain region, Trans. ASME J. Appl.

Mech. 49, 721–727 (1982)

18.15 O. Watanabe, S.N. Atluri: Constitutive modeling

of cyclic plasticity and creep using internal time

concept, Int. J. Plasticity 2,107–134(1986)

18.16 T. Inoue, N. Ohno, A. Suzuki, T. Igari: Evaluation

of inelastic constitutive models under plasticity-

creep interaction for 2.1/4Cr-1Mo steel, Nucl. Eng.

Des. 114, 295–309 (1989)

Part E 18

Continuum Constitutive Modeling References 1031

18.17 T. Inoue, F. Yoshida, N. Ohno, M. Kawai, Y. Ni-

itsu, S. Imatani: Evaluation of inelastic constitutive

models under plasticity-creep interaction in mul-

tiaxial stress state, Nucl. Eng. Des. 126, 1–11 (1991)

18.18 T. Inoue, S. Imatani, Y. Fukuda, K. Fujiyama,

K. Aoto, K. Tamura: Inelastic stress-strain response

for notched specimen of 2.1/4Cr-1Mo steel at 600

◦

Nucl. Eng. Des. 150, 129–139 (1994)

18.19 A.J.M. Spencer: Theory of invariants. In: Continuum

Physics, Vol. 1, ed. by C.A. Eringen (Academic, New

York 1971) pp. 239–353

18.20 A.J.M. Spencer: Isotropic invariants of tensor func-

tions. In: Application of Tensor Functions in Solid

Mechanics, CISM Courses and Lectures, Vol. 292, ed.

by J.P. Boehler (Springer, Berlin, Heidelberg 1987)

pp. 141–169

18.21 R. Hill: The Mathematical Theory of Plasticity (Ox-

ford Univ. Press, Oxford 1950)

18.22 Y. Tomita, A. Shindo: Onset and growth of wrin-

kles in thin square plates subjected to diagonal

tension, Int. J. Mech. Sci. 30, 921–931 (1988)

18.23 S. Imatani, T. Saitoh, K. Yamaguchi: Finite element

analysis of out-of-plane deformation in laminated

sheet metals based on an anisotropic plasticity

model, Mater. Sci. Res. Int. 1, 89–94 (1995)

18.24 A. Phillips, R. Kasper: On the foundation of ther-

moplasticity: An experimental investigation, Trans.

ASME J. Appl. Mech. 40, 891–896 (1973)

18.25 E. Shiratori, K. Ikegami: Studies of the anisotropic

yield condition, J. Mech. Phys. Solids 17,473–491

(1969)

18.26 A. Baltov, A. Sawczuk: A rule of anisotropic hard-

ening, Acta Mech. 1, 81–92 (1965)

18.27 J.F. Williams, N.L. Svensson: A rationally based

yield criterion for workhardening materials, Mec-

canica 6, 104–114 (1971)

18.28 D.W.A. Rees: The theory of scalar plastic deforma-

tion function, Z. Angew. Math. Mech. 63, 217–228

(1983)

18.29 D.C. Drucker: Relation of experiments to mathe-

matical theories of plasticity, Trans. ASME J. Appl.

Mech. 16, 349–357 (1949)

18.30 P. Mazilu, A. Meyers: Yield surface description of

isotropic materials after cold prestrain, Ing. Archiv.

55, 213–220 (1985)

18.31 S. Imatani, M. Teraura, T. Inoue: An inelastic

constitutive model accounting for deformation-

induced anisotropy, Trans. JSME A 55, 2042–2048

(1989)

18.32 N. Ohno, J.D. Wang: Kinematic hardening rules

with critical state of dynamic recovery, Int. J. Plas-

ticity 9, 375–390 (1993)

18.33 F.D. Fischer, M. Berveiller, K. Tanaka, E.R. Ober-

ainger: Continuum mechanical aspects of phase

transformations in solids, Arch. Appl. Mech. 64,

54–85 (1994)

18.34 F.D. Fischer, G. Reisner, E. Werner, K. Tanaka,

G. Cailletaud, T. Antretter: A new view on trans-

formation induced plasticity (TRIP), Int. J. Plasticity

16, 723–748 (2000)

18.35 Q.P. Sun, K.C. Hwang: Micromechanics modelling

for the constitutive behavior of polycrystalline

shape memory alloys, J. Mech. Phys. Solids 41,1–17

(1993)

18.36 T. Inoue, Z.G. Wang: Coupling between stresses,

temperature and metallic structures during pro-

cesses involving phase transformation, Mater. Sci.

Tech. 1, 845–850 (1985)

18.37 R.M. Bowen: Theory of mixture. In: Continuum

Physics, Vol. 3, ed. by C.A. Eringen (Academic, New

York 1976) pp. 1–127

18.38 T. Inoue, T. Yamaguchi, Z.G. Wang: Stresses and

phase transformations occurring in quenching of

carburized steel gear wheel, Mater. Sci. Tech. 1,

872–876 (1985)

18.39 P. Ding, T. Inoue, S. Imatani, D.Y. Ju, E. de Vries:

Simulation of the forging process incorporating

strain-induced phase transformation using the

ﬁnite volume method (Part I: Basic theory and nu-

merical methodology), Mater. Sci. Res. Int. 7,19–26

(2001)

18.40 S. Bhattacharyya, G. Kehl: Isothermal transforma-

tion of austenite under externally applied tensile

stress, Trans. ASM 47, 351–379 (1955)

18.41 M. Fujita, M. Suzuki: The effect of high pressure on

the isothermal transformation in high purity Fe-C

alloys and commercial steels, Trans. ISIJ 14, 44–53

(1974)

18.42 S.V. Radcliffe, M. Schatz: The effect of high pressure

on the martensitic reaction in iron-carbon alloys,

Acta Metall. Mater. 10, 201–207 (1962)

18.43 W.A. Johnson, F.R. Mehl: Reaction kinetics in pro-

cesses of nucleation and growth, Trans. AIME 135,

416–458 (1939)

18.44 C.L. Magee: The nucleation of martensite. In: Phase

Transformation, ed. by H.I. Aaronson (ASM Int.,

Metals Park 1968)

18.45 K. Shinagawa, H. Nishikawa, T. Ishikawa, Y. Hosoi:

Deformation-induced martensitic transformation

in type 304 stainless steel during cold upsetting,

Iron Steel 3, 156–162 (1990)

18.46 P. Ding, D.Y. Ju, T. Inoue, S. Imatani, E. de Vries:

Simulation of the forging process incorporating

strain-induced phase transformation using the ﬁ-

nite volume method (Part II: Effects of strain rate

on structural change and mechanical behavior),

Mater. Sci. Res. Int. 7, 27–33 (2001)

18.47 G.I. Taylor: Plastic strain in metals, J. Inst. Metals

62,307–324(1938)

18.48 C. Teodosiu: Dislocation modelling of crystalline

plasticity. In: Large Plastic Deformation of Crys-

talline Aggregates, CISM Courses and Lectures,

Part E 18

1032 Part E Modeling and Simulation Methods

Vol. 376, ed. by C. Teodosiu (Springer, Wien, New

York 1997) pp. 21–80

18.49 E.C. Aifantis: The physics of plastic deformation,

Int. J. Plasticity 3, 211–247 (1987)

18.50 V. Tvergaard: Inﬂuence of grain boundary sliding

on material failure in the tertiary creep range, Int.

J. Solids Structures 21, 279–293 (1985)

18.51 R.J. Asaro: Micromechanics of crystals and poly-

crystals,Adv.Appl.Mech.23, 1–115 (1983)

18.52 J.W. Hutchinson: Bounds and self-consistent esti-

mates for creep of polycrystalline materials, Proc.

R. Soc. London A 348,101–127(1976)

18.53 J.R. Rice: Inelastic constitutive relations for solids,

J. Mech. Phys. Solids 19, 433–455 (1971)

18.54 T.Inoue,S.Torizuka,K.Nagai,K.Tsuzaki,T.Ohashi:

Effect of plastic strain on grain size of ferrite trans-

formed from deformed austenite in Si-Mn steel,

Mater. Sci. Tech. 17, 1580–1588 (2001)

18.55 E. van der Giessen, V. Tvergaard: A creep rupture

model accounting for cavitation at sliding grain

boundaries, Int. J. Fracture 48, 153–178 (1991)

18.56 G. Beer: An isoparametric joint/interface element

for ﬁnite element analysis, Int. J. Num. Meth. Eng.

21,585–600(1985)

18.57 R. Kawakami, S. Imatani, R. Maeda: Effects of

crystal grain and grain boundary sliding on the

deformation of polycrystal, J. Soc. Mater. Sci. Jpn.

52, 112–118 (2003)

Part E 18

1033

Finite Elemen

19. Finite Element and Finite Difference Methods

Finite element methods (FEM) and ﬁnite differ-

ence methods (FDM) are numerical procedures for

obtaining approximated solutions to boundary-

value or initial-value problems. They can be

applied to various areas of materials measurement

and testing, especially for the characterization

of mechanically or thermally loaded speci-

mens or components. (Experimental methods

for these ﬁelds have been treated in Chaps. 7

and 8.)

The principle is to replace an entire contin-

uous domain of a body of interest by a number

of subdomains in which the unknown function

is represented by simple interpolation functions

with unknown coefﬁcients. Thus, the original

boundary-value problem with an inﬁnite number

of degrees of freedom is converted into a prob-

lem with a ﬁnite number of degrees of freedom

approximately.

19.1 Discretized Numerical Schemes

for FEM and FDM...................................1035

19.2 Basic Derivations in FEM and FDM ..........1037

19.2.1 Finite Difference Method (FDM) ......1037

19.2.2 Finite Element Method (FEM) .........1038

19.3 The Equivalence

of FEM and FDM Methods ......................1041

19.4 From Mechanics to Mathematics:

Equilibrium Equations

and Partial Differential Equations..........1042

19.4.1 Heat Conduction Problem

in the Two-Dimensional Case ........1043

19.4.2 Elastic Solid Problem

in the Three-Dimensional Case ......1043

19.5 From Mathematics to Mechanics:

Characteristic of Partial Differential

Equations ............................................1047

19.5.1 Elliptic Type .................................1048

19.5.2 Parabolic Type .............................1048

19.5.3 Hyperbolic Type ...........................1048

19.6 Time Integration for Unsteady Problems. 1049

19.6.1 FDM............................................1049

19.6.2 FEM ............................................1050

19.7 Multidimensional Case ..........................1051

19.7.1 Finite Difference Method...............1051

19.7.2 Finite Element Method .................1052

19.8 Treatment of the Nonlinear Case............1055

19.9 Advanced Topics in FEM and FDM ...........1055

19.9.1 Preprocessing ..............................1055

19.9.2 Postprocessing.............................1056

19.9.3 Numerical Error............................1056

19.9.4 Relatives of FEM and FDM ..............1057

19.9.5 Matrix Calculation

and Parallel Computations ............1057

19.9.6Multiscale Method........................1059

19.10 Free Codes ...........................................1059

References ..................................................1059

A ﬁnite element analysis includes the following basic

steps.

•

Discretization or subdivision of the domain of the

body of interest.

•

Selection of the interpolation functions to provide

an approximation of the unknown solution within

an element.

•

Formulation of the system of equations (using, for

example, the typical Ritz variational and Galerkin

methods).

•

Solution of the system of equations. If the sys-

tem of equations is solved, the desired physical

variables at the nodes can then be computed and

the results displayed in form of curves, plots, or

color pictures, which are more meaningful and

interpretable.

Finite element methods (FEM) use a complex system of

points called nodes, which make up a grid called a mesh.

This mesh is programmed to contain the material and

structural properties that deﬁne how the structure will

Part E 19

1034 Part E Modeling and Simulation Methods

react to certain loading conditions. Nodes are assigned

at a certain density throughout the material depending

on the anticipated stress levels of a particular area. Re-

gions that receive large amounts of stress usually have

a higher node density than those that experience little or

no stress. Points of interest may concern, for example,

fracture points of previously tested materials (Sect. 7.4)

or high-stress areas. The mesh acts like a spider web in

that from each node, there extends a mesh element for

each of the adjacent nodes. This web of vectors is what

carries the material properties to the object, creating

many elements.

As a typical example, consider a material or a body

in which the distribution of an unknown variable – such

as temperature or displacement – is required. The ﬁrst

step of any ﬁnite element analysis is to divide the actual

geometry of the structure using a collection of discrete

portions called ﬁnite elements. The elements are joined

together by shared nodes. The collection of nodes and

ﬁnite elements is known as the mesh (Fig. 19.1).

The variable to be determined in the analysis is as-

sumed to act over each element in a predeﬁned manner.

The number and type of elements is chosen to ensure

that the variable distribution over the whole body is

adequately approximated by the combined elemental

representations. After the problem has been divided into

discrete units, the governing equations for each element

are formulated and then assembled to give a system

equations that describe the behavior of the body as

a whole.

Finite difference methods (FDM) are also based on

the similar idea. In this method, the elements and mesh

are called grids and grid, respectively (Fig. 19.2). The

major differences of FDM from FEM are (1) Governing

partial differential equations are approximated directly

by ﬁnite difference approximation, not by interpolation

functions nor via the Galerkin method, (2) The dis-

cretized whole domain is not covered by a ﬁnite number

of interpolation functions, but is represented by a ﬁnite

number of sampling nodes. A ﬁnite difference analysis

can be summarized as follows.

•

Discretization or subdivision of the domain of the

body of interest

•

Selection of the ﬁnite difference functions to

provide an approximation of the derivatives of un-

known solution at a node

Body of

interest

Nodes

Elements

Mesh

Fig. 19.1 Collection of nodes and ﬁnite elements (mesh) in

FEM

Body of

interest

Nodes

Grids

Grid

Fig. 19.2 Collection of nodes and grids in FDM

•

Formulation of the system of equations at all the

nodes

•

Solution of the system of equations. If the system of

equations is solved, the desired physical variables

at the nodes can then be computed and the results

displayed in form of curves, plots, or color pictures,

which are more meaningful and interpretable, just

as FEM does.

In this chapter, both FEM and FDM are explained.

The contents include the basics necessary to under-

stand the formulation, more advanced topics, practical

information, and information on the use of commercial

software. This chapter is for nonspecialists in macro-

scopic mechanics and computational mechanics, so

only the key concepts are focused on in the discus-

sion. Since the reader might not be accustomed to tensor

representation, all derivations are shown without tensor

notations. The explanation will focus on a steady or un-

steady heat conduction problem in the one-dimensional

case, however, multidimensional and nonlinear cases

are also handled within the limited pages.

Part E 19

Finite Element and Finite Difference Methods 19.1 Discretized Numerical Schemes for FEM and FDM 1035

19.1 Discretized Numerical Schemes for FEM and FDM

FEM and FDM are classiﬁed into so-called discretized

numerical schemes. Discretized numerical schemes are

approximated ways to solve boundary problems, de-

scribed by partial differential equations, with complex

shaped domain with complicated boundary conditions,

which can usually not be solved analytically by hand

calculation. It should be noted that FEM and FDM are

not applicable to phenomena that cannot be represented

by partial differential equations.

In this section, we will focus on the process of dis-

cretized numerical schemes by using a simple example

of a two-coupled spring model.

The relationship between a displacement u and

a force F can be connected with a spring constant k as

in the following equation, based on the knowledge of

high school’s physics

ku = F . (19.1)

Fig. 19.3 Simple spring model

For the preparation to consider a more complicated

case, we’d like to extend (19.1) for a more general case

by introducing a spring model element, as shown in

Fig. 19.4.

Fig. 19.4 Spring model element

By introducing u

, F

, u

,andF

at nodes 1 and 2,

the spring model element can be represented with ma-

trix notation



k −k

−kk









. (19.2)

By setting u

= 0 as a boundary condition, this matrix

system is reduced to (19.1) in the following fashion



k −k

−kk









, (19.3)



−k











. (19.4)

Since the ﬁrst row cannot be used for calculation due to

the unknown value of F

, we can get the same equation

as (19.1) by omitting the ﬁrst row in (19.4)

= F

. (19.5)

This is just a warm-up explanation of the two-coupled

spring model below.

Now let’s think about the following two-spring

coupling model.

Fig. 19.5 Two-coupled spring model

The two-coupled spring model can be divided into two

spring-model elements as shown in Fig. 19.6.

By combining two matrices in Fig. 19.6,wegetthe

following system of equations

⎛

⎜

⎝

−k

000

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

00 0

0 k

−k

0 −k

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

→

⎛

⎜

⎝

−k

0 −k

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

(19.6)

Part E 19.1

1036 Part E Modeling and Simulation Methods

–k

Fig. 19.6 Two elements of the two-spring model

Since the inner force is always balanced, F

=0.

Then,

⎛

⎜

⎝

−k

0 −k

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

. (19.7)

However, the combined matrix above cannot be in-

verted due to the singularity of the matrix, which is

physically interpreted such that the position of the bal-

ancing springs in air is not determined explicitly. We

need a ﬁxed point to solve this matrix, which is called a

(ﬁrst) boundary condition. Suppose that node 1 is ﬁxed

as u

=0. In this case, the matrix becomes solvable in

the following way.

1. Substitute u

=0 into the matrix system

⎛

⎜

⎝

−k

0 −k

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

⎛

⎜

⎝

⎞

⎟

⎠

; (19.8)

2. Omit the ﬁrst column in the matrix

⎛

⎜

⎝

−k

⎞

⎟

⎠





⎛

⎜

⎝

⎞

⎟

⎠

; (19.9)

3. Since F

is an unknown value, the ﬁrst row is omit-

ted. Then we get a 2 × 2 matrix system that is now

solvable



−k









. (19.10)

By solving the matrix system of (19.10), we get the

following solutions, which is equivalent to the formula

of the two-coupled spring model 1/k =1/k

+1/k

. (19.11)

The procedure above can be summarized as follows.

1. Discretize the targeted domain into basic element

models.

2. Construct a balancing (equilibrium) equation (local

matrix) at each basic element model.

3. Combine all local matrices, add boundary condi-

tions, and then solve global the matrix.

The process explained above is common among

so-called discretized numerical schemes such as FEM

and FDM. Usually the spring constant above becomes

a spring coefﬁcient matrix in these methods. The only

difference between FEM and FDM can be summarized

in the way to derive the spring matrix part. We’ll treat

basic derivations in FEM and FDM in the following

section.

Part E 19.1

Finite Element and Finite Difference Methods 19.2 Basic Derivations in FEM and FDM 1037

19.2 Basic Derivations in FEM and FDM

The governing equation of steady heat conduction prob-

lems in the one-dimensional case is given by the

following equation, where u is a temperature, k is a heat

conduction coefﬁcient, and f is a heat source

∂

∂x

=−f . (19.12)

(0)= b

f = const.

k = const.

u(1) = a

01x

Fig. 19.7 One-dimensional steady heat conduction prob-

lem

Suppose that the boundary conditions (please recall

Sect. 19.1) are provided as follows

u(1) =a (u = a at x =1) ,

(0) =b



∂u

∂x

=b at x =0



. (19.13)

Here, u(x

) =a is called either the ﬁrst boundary con-

dition or the Dirichlet boundary condition, while the

boundary condition of u

) =b is named the second

boundary condition or Neumann boundary condition.

The problem represented by (19.12)and(19.13)is

deﬁned as a boundary value problem. (For unsteady

problems, please refer to Sect. 19.6.)

Of course, the analytical solution of this problem

can be given in the following form in this case, but

we would like to consider the approximated numerical

simulation by FDM and FEM in the following parts

u(x) =a +(x −1)b+



⎡

⎣



f (z)dz

⎤

⎦

dy . (19.14)

The exact solution function with respect to (wrt) the

x-coordinate, under the parameters of a =1, b = 0.2,

f =1, h = 0.5, and k = 1, is

u(x) =−0.5x

+0.2x +1.3 ,



u(0) =1.300, u(0.5) =1.275, u(1) =1.000



(19.15)

19.2.1 Finite Difference Method (FDM)

In FDM, the derivatives of the governing equation are

replaced by so-called ﬁnite difference approximations.

u(x–h)

n –1

x –h

u(x)

u(x+h)

n +1

x +h

Fig. 19.8 Finite difference approximation of the ﬁrst

derivative of u wrt x

To approximate ﬁrst derivative ∂u/∂x at node n in

Fig. 19.8, there exist three schemes in FDM. The central

difference scheme in (19.16) guarantees second-order

accuracy with respect to x, and ﬁrst-order accuracy in

the other two one-sided methods in (19.17)and(19.18)

[19.1].

Central difference scheme

∂u

∂x

= lim

h→0

u(x +h) −u(x −h)

≈

u(x +h)−u(x −h)

; (19.16)

Backward difference scheme

∂u

∂x

= lim

h→0

u(x) −u(x −h)

≈

u(x) −u(x −h)

;

(19.17)

Forward difference scheme

∂u

∂x

= lim

h→0

u(x +h)−u(x)

≈

u(x +h) −u(x)

(19.18)

u(x+h)–u(x)

n –1

x –h

n +1

x + h x

u(x)–u(x–h)

Fig. 19.9 Finite difference approximation of the second

derivative of u wrt x

Part E 19.2