finite element models and
constitutive models for large deformation
elasto-plasticity
ERIK VAN DER GIESSEN
MODELS IN NONLINEAR THERMOMECHANICS
finite element models and
constitutive models for large deformation
elasto-plasticity
finite element models and
constitutive models for large deformation
elasto-plasticity
PROEFSCHRIFT
'■ x - : ..S
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, prof.dr. J.M. Dirken, in het openbaar te verdedigen ten overstaan van een commissie door het College van Dekanen daartoe aangewezen, op donderdag
24 september 1987 te 14.00 uur
door
ERIK VAN DER GIESSEN geboren te Rotterdam werktuigkundig ingenieur
9 Prome+h.e.MPoIeinl w |
y
MODELS IN NONLINEAR THERMOMECHANICS
van
Erik van der Giessen
1. Nadat experimenteel onderzoek in de technische mechanica door de opkomst van numerieke methoden in de laatste decennia sterk op de achtergrond is geraakt, is thans de verdere ontwikkeling van numerieke methoden afhankelijk van nauwgezet experimenteel onderzoek.
2. Hoewel dualiteit universeel herkend wordt als een sleutelbe grip in de mechanica, wordt het, dit fysisch begrip geheel dekkende, wiskundige begrip dualiteit in de wiskundige be schrijving ervan veelal onbenut gelaten.
3. De spanningstensor van Cauchy is principieel een gemengde tensor (par. 4.1 van dit proefschrift).
4. De plastische spin in een elastisch-plastisch materiaal is altijd nul, tenzij de thermodynamische toestand mede bepaald wordt door een essentieel niet-symmetrische tensoriële toe standsgrootheid .
5. Certificatie van gebruikers van commerciële eindige- element-enprogrammatuur i.v.m. modelvorming, gebruik en interpretatie is symptoombestrijding. De oorzaak zou worden opgeheven door de benodigde fysische kennis d.m.v. bijvoorbeeld kennis systemen in de programmatuur op te nemen.
6. De ontwikkeling van eindige-elementenprogramma1s in de pro grammeertaal C, zoals tegenwoordig regelmatig wordt voor gesteld, zou een staps terugwaarts betekenen in de trend naar hiërarchisch gestructureerde, onderhoudsvriendelijke en flexibele programmatuur.
7. Bij doctoraal examens zou aan de gangbare beoordeling van verslaggeving en theoretische en experimentele vaardigheden, een beoordeling m.b.t. informatica-aspecten moeten worden toegevoegd.
8. De verlenging van ronde proefstaven bij grote elastisch-plastische torsievervormingen wordt in hoofdzaak veroorzaakt door anisotropie. In verband hiermee, vormt de experimentele bepaling van deze grootheid een belangrijk middel bij de beoordeling van plasticiteitsmodellen of bij de bepaling van de modelparameters.
9. De meest efficiënte manier om het getal % in decimale vorm te
onthouden is d.m.v. de breuk 355/113.
This thesis considers two subjects concerning modelling in thermomechanics, namely finite element models of continua and constitutive models for large deformation elasto-plasticity. The common mathematical language used is that of modern linear algebra.
In the first part, finite element models are- discussed for the discrete description of the thermomëchanical behaviour of continua. The approach constitutes an alternative to commonly adopted Galerkin methods in conjunction with, numerical integra tion. In addition to the nodal point quantities governing the interaction between elements, internal state quantities are defined in 'a finite number óf so-called material sampling points within the element. The fully algebraic general thermomëchanical equations of the discrete model are formulated directly using the principles of virtual power and of virtual heat flow.
The application to problems of large deformation elasto-plasticity both according to a material and a spatial approach is discussed. Two-dimensional problems of incompressible viscous fluid flow and of heat flow illustrate the theory. The so-called redundant pressure method is presented as a solution method for handling the incompressibility constraints.
In the second part of this thesis,.constitutive models for large deformation elastó-plasticity are presented. Starting, point is the general constitutive theory for materials possessing. a natural reference state. The key proposition involved is that.the plastic deformation process is governed by 'a specific stress tensor which characterises the current stress state in the geometry in the natural reference state. Using . a strictly deterministic approach, implying a generalized normality condi tion for the plastic rate of deformation, constitutive equations are formulated first for nonharderiing anisotropic elastic-plastic materials, taking due account of possibly large elastic strains.
Next, based upon the concept of a plastically induced orientational structure, a relatively simple and reasonably specific constitutive theory is presented which accounts ex plicitly for the progressive development of anisotropy during plastic deformation processes. With one scalar and one essential ly nonsymmetric so-called micro-stress tensor, which is dual to the tensor governing the evolution of this orientational structure, as additional internal local state variables charac terizing the thermodynamic reference state, a wide range of induced anisotropic behaviour is covered. Particular.examples of the characteristic material functions involved7are presented. .
A thermodynamically consistent, large deformation, kinematic. hardening model emerges from the theory as a special case. Various aspects of the theory are demonstrated by consideration of large simple shear deformation processes. The implications of the theory on the notion of a plastic spin are amply discussed.
page
1 INTRODUCTION 1
2 BASIC NOT(AT)IONS 4 1 Vectors and linear maps 4
2 Duality 5 3 Tensors as linear maps 6
4 Dual and transpose tensors 7 5 Inverse tensors, invariants and deviators 9
6 Physical space 9
3 DISCRETE THERMOMECHANICS OF CONTINUA 12
1 Continuum models of matter 12 2 A discrete formulation 13
1 Geometry 13 2 Spatial interactions: Flow of mass and entropy 15
3 The local state: The concept of material
sampling points 18 1 Specification of material sampling points 20
2 Material sampling points in complex elements 22 4 Determinism in finite element thermomechanics 24 5 The discrete equations of thermomechanics 25
1 Material derivatives in spatial descriptions 25
2 Deformations of finite elements 26 3 The balances of mass and of entropy 28 4 The principles of virtual power and heat flow 29
5 The energy balance 36 6 Constitutive equations 37
3 A few applications 37 1 Large deformation elasto-plasticity 37
2 Incompressible Newtonian flow 40
1 Governing equations 40 2 Trying to solve the problem ... 42
3 The Redundant Pressure Method 44
4 Applications 48 3 Heat transfer 52
1 Governing equations 52
2 An example 56 4 Concluding remarks 59
ELASTO-PLASTICITY 61 1 Some preliminary remarks 62
2 The theory of the natural reference state 64 1 Continuous media and the natural reference state 64
1 The concept of a local natural reference
state 64 2 Kinematics 67
2 Thermodynamics 71 1 The local state 71 2 Entropy flow processes 74
3 Irreversible deformation processes 79 3 Summary of the constitutive theory 84 3 The natural reference state and elasto-plasticity 86
1 A medium with a natural reference state as
a model of an elastic-plastic material 86 2 Constitutive equations for time-independent
plasticity 86 3 Constitutive equations for isotropic materials 91
4 The natural reference vs. intermediate state 93
5 A simple illustration 94 4 A constitutive theory for induced anisotropic
hardening 96 1 The thermo-elastic reference state 97
2 The plastically induced orientational structure 99
3 Thermodynamics 101 4 Anisotropic hardening described by present model 106
1 A dichotomy 106 2 Explicit anisotropic hardening 107
3 Implicit anisotropic hardening 108 5 On constitutive relations for the internal
state variables 112 1 Constitutive relations for the plastic
state variables 112 2 Constitutive relations for the induction
tensor 117 6 Summary of the constitutive theory 119
7 Illustration of the theory 121 8 A special case: Kinematic hardening 127
1 Constitutive equations 127 2 Large simple shear: A case study 131
5 Comparison with other approaches 136
APPENDIX: Element matrices for quadratic
triangular elements 140
REFERENCES 143
1. INTRODUCTION
Industrial forming processes used for manufacturing usually complex shaped products rest on the possibility of the materials used — usually metals — to sustain very large deformations. Also for instance/ in ductile materials, considerable elastic-plastic deformations may occur in preceeding failure. In order to improve product quality and control of forming processes, and knowledge on the. failure phenomena etc., physical mathematical models are needed for the simulation of the thermomechanical processes involved.
Today, appropriate models for that purpose are conveniently established in two steps. Firstly, the actual corpuscular material is replaced with a continuous medium and secondly, this continuum is replaced with an assembly of finite elements of the continuum. The mathematical description of the macroscopic behaviour of the actual discrete, corpuscular material is then obtained in terms of the discrete finite element model. The macroscopic material behaviour is taken into account by means of constitutive models considered appropriate in terms of the continuum model. The finite element model of the continuous medium must ensure that the general basic physical laws of thermomechanics are satisfied.
In this work we shall deal with both finite element models for continuum thermomechanics and constitutive models for large deformation elasto-plasticity. Both subjects are definitely not new branches within contemporary applied mechanics; nevertheless they are subjected to rapid developments. Constitutive modelling of large deformation plasticity in particular seems to be still in its infancy.
Restricting the discussion to continuum thermomechanical problems, two approaches to finite element methods may be distinguished which differ in basic conception. In the first -and almost common approach — the structure or body is regarded as an essentially continuous system for which the class of admissible continuous fields for the physical quantities involved is replaced with a restricted, finite set of continuous functions defined in subdomains of the system. Subsequently, the basic laws of thermomechanics are enforced in an approximate sense by means of a mathematically motivated principle of Galerkin. In this view, the finite element method is a mere mathematical tool for obtaining approximate solutions of the governing differential field equations.
In a second, alternative approach the structure or body is discretized from the onset. Finite element models within such an approach is the subject of Chapter 3 of this work. Each finite element is considered a distinct thermomechanical system, the interactions of which with adjacent elements are taken into account by means of pertinent nodal point quantities, while the
(internal) state quantities are introduced at material sampling points within the element. It will be shown that by means of the physically motivated principle of virtual work and virtual heat flow, the basic laws of thermomechanics can directly be formulated for this discrete model in a consistent way. The final governing system of general finite element equations possess a
fine, computational-friendly hierarchic structure, which no longer demands tedious numerical integration.
The application of the finite element theory presented to problems of large deformation elasto-plasticity will be dis cussed. Further, the method will be demonstrated by means of ■ applications to two-dimensional problems of heat flow and of slow stationary flow of incompressible viscous fluids. For the latter problems a new solution method will be presented which deals with the troublesome incompressibility constraints.
The first attempts at modelling elastic-plastic material behaviour at large deformations (the word deformation implying both strain and rotation) within a continuum framework date back to the 1960's. Nevertheless, even today, this topic is characterized by considerable disagreement over proposed ap proaches and consensus seems to have been achieved only in so far as only the most elementary aspects of the behaviour are concerned.
As compared with classical elasto-plasticity at very small strains, the mathematical description of large deformation elasto-plasticity is conceptually much more sophisticated due to the essential geometric nonlinearitiès involved and the increased complexity of the processes taking place at large deformations. This first feature demands the constitutive theory to be formulated in a general nonlinear continuum thermomechahical framework. In order of increasing complexity, one may distinguish between the following models of the actual behaviour of materials
(polycrystalline metals in particular) at large deformations:
- plastic flow without hardening, with initial anisotropy (perfect plasticity),
- hardening with invariable directions of anisotropy, - development of anisotropy (texture).
(This list is far from being complete; in this work, for instance time-independent phenomena and the effect of microstructural defects are excluded).
Starting out from a constitutive theory based upon standard thermodynamic concepts and on the concept of a natural reference state, we shall first present in Chapter 4 a constitutive description of non-hardening, large deformation èlasto-plasticity with initial anisotropy. Subsequently we will discuss a theory which accounts for the development of anisotropy in the course of the deformation process (induced anisotropic hardening). This theory is based upon the concept of a plastically induced orientational structure which attempts to incorporate the development of preferred directions into the continuum model. A kinematic hardening model emerges from the theory as a special case.
Various aspects of the theories proposed will be illustrated by considering the behaviour in large simple shear — a deformation process which involves both large strains and large rotations. The description of plastically induced anisotropy will be discussed in quite some detail. Considerable attention will also be paid to the controversial matter concerning the notion of plastic spin.
notions and notations of modern linear algebra used throughout this work. The mathematical tools involved appear to be most convenient in expressing the nature of discrete, finite element models. Furthermore, as tools in constitutive modelling, they provide a mathematically more consistent yet more explanatory alternative to the currently popular polyadic notation of vectors and tensors.
2. BASIC NOT fAT)IONS
Many of the arguments are carried in the language of linear algebra (Bowen and Wang, 1976; Dodson and Poston, 1977). The intrinsic (coordinate-free) notation of vectors and tensors combines generality with clearness and concision. In this chapter we recapitulate some of the notions involved, mainly for reasons of notation.
2.1 Vectors and linear maps
Vector spaces will be denoted by italics and vectors (as
elements of vector spaces) by bold face letters, e.g. UÊ£7. Only finite dimensional vector spaces will be considered. Rn is the
standard real n-space. The coordinate vector of a.vector usU in a
chosen basis {fi-^U), i.e. the (column) matrix [u1] of components
of u=u1/ii is denoted by u. The summation convention applies for
repeated diagonal indices.
A (linear) map A from a vector space U into a vector space
v>
A : Ü —> V : u i—> v:=Au ;
is again a vector. The vector space of all linear maps U—»V is
denoted by Lin(U;V) and dim Lin (U,V) =dim U dim V, The image (or
range) and kernel (or null space) of A are denoted by Im A and Ker A respectively. In general,
dim U > dim Im A < dim V,
dim U = dim Im A + dim Ker A.
The rank of the matrix A corresponding to A is given by Rank A := dim Im A.
An inner product space U is a vector space with an inner
product map U*U—>R, which we will denote by a dot, e.g. u*v
(u,ve/7) .
Subspaces Ulr..,Un of some vector space may be combined to
form other spaces through the operation of summing, e.g. n
U := ^ ui' (2.2)
i = l
If the intersection of any two of the subspaces is the null space, such a sum is a direct sum
n
U := © C/j. (2.3)
i = l
Inversely, a particular subspace of the vector space U defined by
(2.2) or (2.3) can be selected by one of the (canonical) (2.1)
projections
n
Pk : U - > £7k : u= ^ U j ^ Pk( u ) = uk. i = 1
In the remainder of this chapter, the vector space X is an
inner product space of dimension dim X=n with {gi> (i=i,..,n) a given basis in X.
2.2 Duality
The dual space U* of a vector space U is the vector space Lin(U;R) of all linear functional^
u* : U —» R : u i—> u*u. (2.4)
The suggestive notation <u*,u):=u*u, termed the scalar product,
is also used. The dual of U* is U**=U and dim C7*=dim U. If VcU,
the subspace VxcU* given by
V1 := {u*| <u*,v>=0, VveV}
is termed the orthogonal complement of 7.
Let C7j_, .., C7n be subspaces of U such that U is their direct sum, Eq. (2.3), and let U1*, . . ,Un* be their dual spaces. The
scalar product on U then is the sum of the scalar products on Uit
<u,u*> = < I u j , ^U j* ) = J (Uj,^*). (2.5) i = 1 j = 1 i = 1
The basis {g4} in X induces a dual basis (g1) in X* through the scalar product
(g'^gj) = {g^g1) = V j (2.6)
where S1* is the Kronecker delta symbol. Furthermore, the metric on X defined by the inner product, induces two isomorphisms
G+ : X —» X*., Gf : X* —> X (2.7)
such that
<G*x ,y ) = <x ,G^y > = x -y , Vx ,y eX ;
<Gfx*,y*> = <x*;Gfy*> = x*«y* , Vx*,y*eX* (2.8) This particular notation is used because G^, and G* may be viewed
upon as operators which respectively "lower and raise the indices".
The partial derivative of a (nonlinear) function f : C7x • •
xV-->i? with respect to e.g. vectors ue£7 at some (u0, . . ,v0) eUx • • xV is a linear map
3f
— (u0,..,v0) : ü —> R (2.9)
9u
i.e. a vector in Lin(U;R)=U* (the point of application (u0,..,v0) will be omitted if there is no possibility of confusion) . With {/ij } and {/i1} denoting dual bases in U and U*, we have in terms of components
. 3f 8f =
u = u'/ij , — = — fil (2.10)
3u 3U1
2.3 Tensors as linear maps
The tensor product, denoted by ®, of vector spaces Ux,..,Un
is again a vector space, namely
U^®*'®Un := Lin{Ux* , . . , C7n* ;JR) .
The dimension of this tensor space is dim(C710«'®C7n) = dim ü1«,dim
Un. Tensor spaces of the type
[7®« «®C7 0 C7*®« «®[7* lc times h times
are said to be contravariant of order k and covariant of order h,
or of type (k h), and are denoted by khE7 (with proper extension to other orderings of the tensor product). With {gj} and (g1} being dual bases in X and X*, the set of tensor products
{g. 0. .gj ®g ^..fcg h; il,. . ik, il t- . , i h = l,..,n}
forms a basis for the tensor space k^X. The isomorphisms defined
in (2.7) can be expressed as
G ^ g ^ g ^ g ' , gij=gi«gj; G ^ g ^ g ^ g j , gij=gi-gj.
By virtue of the natural isomorphism
Lin (U;V) ~ V®U* (2.11)
for any two vector spaces U and V, the linear map A : U—»V can be
identified with the tensor AeV®U*. The particular order of the tensor product in (2.11) is adopted for the matter of notation: For a map A : X—»X we have in terms of the components of A in the induced tensor basis in * XX , A = A ^ g ^ g ^ ,
x = x'gj i—» y := Ax = A1jxkgi®gigk = A^xJgj
which shows that due to the convention (2.11) gi®g3gk can be interpreted as g4 (g] ,gk) , i.e. the two adjoining base vectors
operate on one another (this is readily extended to the case
where X itself is a tensor space). The identification (2.11)
simplifies the notation and the readability of tensorial expressions.
Furthermore we shall use the natural isomorphism
(£7®V)* ~ U*®V* (2.12)
to identify the dual tensor space of C70V with L7*®V*.
2.4 Dual and transpose tensors
The dual map or tensor corresponding to the map A : U—*V is
the map
A* : V* -» u* : v* i—> A*v* := v*A (2.13)
i m p l y i n g t h a t
<A*v*,u> = <v*,Au> , V ue£7, v*GV*. ( 2 . 1 4 )
Note that, while AeV®U*, A* is an element of the tensor space [7*®V and. not of the dual space (V®U*) *=V*®U. With (MiK (M1) and
{V[)i { P1} dual bases in U and V, we have in terms of components •
A=AijPi®/zj, A*=A*1J|ii®»j : A*jj = A3i
If further B : V—»JV, the dual of the composite map C = BA : U—»(V
—»)F/ is C* = A*B* : W*—> (V*—») U*. If B:U*—»V* then the tensor
BeV*®U=(V®£7*) * is dual to the tensor AeV®U* and the following property holds <A,B) = <A* ,B*) . In the special case that the map A can be represented by a simple tensor, i.e. by a tensor product
A=p®/i* (veU, fi*eU*) , its dual is the tensor A*=/H*®P.
The transpose of a map A : X-*Y between two inner product
spaces is the map AT : Y—>X such that (cf. Eq. (2.14)) ATy x = y A x , V xeX, yzY
(the dot in the left-hand and right-hand side indicate the inner products . on X and Y respectively). Hence, (Y®X*)T = X®5T* . With
{£i>/ (£*} and {i)i}, ( V > dual bases in X and Y respectively, we
have in terms of components :
A=Aijl|i®£J, A^A^jti®!1 : AT»j = i ^ A1^ "1
respectively. If further B : Y—*Z, the transpose of the composite
map C = BA : X—>(Y—>)Z is CT = ATBT : Z-^ (Y—>)X.
The dual and transpose of the map A : X—»F are related through isomorphisms between the dual spaces of the form (2.7). This is shown in the following diagram:
X
«-G,Y
X'
,*T
From this
, * diagram * we deduce e.g. A
=G\A*GJ,
just the that theA*=A*; this Notice
matrix A* of A* in any bases {^eX} and {ii^Y) is
transpose of the matrix A of A in these bases, i.€ does not hold in general for the matrix A~
A map A : X—>X is symmetric if AT
AT = -A (notice that both A and AT are*elements in X®X*). If A is nonsymmetric it can be decomposed uniquely into a symmetric part A+ and an antisymmetric part A" according to
T of A'.
A and antisymmetric if
Af = |(A+A') = |(A-AT)
such that A=A++A". The spaces (X®X*)+ and (X®X*)~ of all symmetric and antisymmetric maps X—»X respectively are subspaces of X®X* and have the following dimensions
dim(X®X*)+ = |n(n+i), dim(X®X*)~ = |n(n-i).
if dimX=n (dimX®X*=n2). In particular, we have the decomposition
x®x* = (x®x*)
+e(x®x*)".
A map A : U—»C7* is called self-dual here if A* = A (notice
that for these kinds of maps A,A*eC7*®L7*) . The isomorphisms defined in (2.7-8) are examples of self-dual maps. Notice that in contrast with symmetric maps as discussed above, self-dual maps
have symmetric corresponding matrices, i.e. A = A*, in any dual bases in U. A self-dual map A is positive (semi-) definite if
2.5 Inverse tensors. invariants and deviators
The map A : X—»X is invertible if the inverse map A- 1 : X—»X
exists such that A" 1Af=AA" ^ 1 , with I the identity map or unit
tensor on X, !=&* ^q®qi .
The inverse exists if the determinant of A, defined through
det A := det[A'j]
as the determinant of the matrix [A*j] of its components A * j , is non-zero. The trace of A is defined by tr A:=Ak k. With dim X=n,
tr A and det A may be used as the first and the last of the n
invariants of A (scalar valued characterizations of A independent
of the basis in X) . The second invariant of A is jtr A2 (A2=AA) ;
it may be shown that tr A2=<A,A*).
The deviatoric part, denoted by a superscript d, of a map A
: X—>X is the map Ad= A - ( V 3 t r A) I. The vector space of all
deviator maps X—-»X is denoted by (X®X*)d and possesses the
properties
r*X d
(X®X*)a c X®X' dim (X®X*)d= n-i
2.6 Physical space
The physical space is a p-dimensional (i<p<3) Euclidean
point space (manifold) ^p with (inner product) translation space
Xp (Bowen and Wang, 197 6 ) . (We shall suppress the superscript p
if the dimension is of no importance.) To any point xe? we may attribute a position vector X ( X ) E X with respect to a fixed point Oe^ in a unique way.
An atlas in .yp with coordinates
basis {gj} at x according to {x1} (i p) , induces a 3x 9 i: = 3x
r
t hso that gj is tangent to the ix n coordinate curve at x. The dual
basis is induced through the scalar product (2.6). In the particular case of a Cartesian coordinate system the dual bases will be denoted by {ej} and {e1}.
A most convenient geometric interpretation of the dual basis is to regard it as a set of hyperplanes in 9 (see Fig. 1 ) . A dual
vector in X* can then be viewed upon as a pattern of parallel hyperplanes and the scalar product used to construct X can be thought of as "the number of surfaces of the dual vector that is pierced by the vector when they operate on one another". Furthermore, a tensor in X®X* can be thought of as a "structure consisting of a plane and an arrow". (Of course a plane can also be represented by its normal depicted again as an arrow. In many texts on continuum mechanics, this has led people to tacitly identifying X with X* (the coordinate-free notation of vectors and tensors thus obtained is known as the polyadic notation). However "duality" is considered here too much an essential notion in physics to be obscured in the mathematical description).
x = const
x — const
Fig. 1. Basis vectors and dual basis vectors for a
coordinate system {x
1}.
Various physical quantities observed when studying phenomena in <f can be represented by vectorial quantities. Mathematically,
such quantities can be given a vector structure in the algebraic sense. It is almost common in treatises on mechanics to regards these vectors as elements of the vector space X of y. It has
hardly ever been realized that from a formal linear algebra point of view physically different vectors should reside in different, properly defined vector spaces (this applies in particular to physically dual vectors). Fortunately however, isomorphisms may be defined between each of these vector spaces and the vector space X (the physical vector spaces may be thought of as being "modelled on X") . In order to avoid an extremely involved and obscure notation, we will do setup properly the necessary vector spaces (e.g. the velocity space V of velocity vectors at some
point xe^) but we shall continually consider them isomorphic to X without any further designation (e.g., V~X and for the space F
dual to V, F~X*) . (The visualization of a force vector (in the
algebraic sense) i&F as a set of parallel hyperplanes, although
quite opposite to the traditional point of view advocated in all elementary and most advanced textbooks on mechanics, is consistent with physics. This will be clear immediately if one thinks of forces derived from a potential; also, in quantum-mechanics, the momentum of a particle is characterized by a pattern of de Broglie wave's surfaces of equal phase (Misner et
al., 1973).) Hence, any operation on a (physical) vector space
can be reduced to an operation on X or X* (or tensor spaces of X)
and, vice versa, any operation defined on X (or more general k hX) can be applied to any vector space modelled on X (khX) using the same symbolism.
In studying continua in y we will have to deal with fields
of vectors and tensors (e.g. the velocity field v(x) or v(xl); although formally v(x) represents a vector only for some xe^, we shall simply write v(x)eV~X). The most important operator on fields we will need here is the gradient operator. For (kh)~ tensor fields it is defined here to be a map
Vx : V - > V ® * * =kh+i * : u — * Vxu
with Vxu the value of the usual gradient of u(x)EkhX at some point xe^,
Vxu = u1 1 ' ' lkj l , , j h.n gi x® « > ® gi l c® g: lQ> «®gJ h® gn; ' l » • • ( ' k / J 1 ' • • ' ^ h ' n ~ l / • • / P
; denoting covariant differentiation. (If there is no possibility of confusion regarding the point of application, we will omit the subscript x). Further, we define here the dual gradient operator
according to (2.13) by
Vx* : X*®khX=1khX —» (khX)* : wi—> Vx*w.
For vector (i.e. (10)-tensor) fields, Vxel 2X and Vx*e1 1 1X can be expressed as
3 d
V„ = dK — r e i ® ej® ek, V „ * = 8 .j — r ei® ei® ek.
x k 8 x3 i x ■ 3 xk J
The divergence of a vector field v(x)eX at x is defined by
3. DISCRETE THERMOMECHANICS OF CONTINUA
The almost common view today on finite element discre tizations for continuum flow and deformation problems regards the finite element method as a mere mathematical tool for obtaining approximate solutions of the governing differential equations
(e.g. Bathe, 1982; Carey and Oden, 1983; Zienkiewicz, 1977). Within each element an approximation of the actual solution is sought within the class of interpolation polynomials of typically low degree on the basis of the nodal values. Commonly, an associated functional, obtained usually by Galerkin's method, then provides the system of generally nonlinear algebraic matrix equations for the unknown nodal quantities.
Except for the simplest elements, the matrices obtained involve rather complicated integrals which are commonly evaluated by numerical integration. Although it seems in the spirit of this view on finite element techniques to aim at minimization of the errors introduced in the numerical integration process, it has been frequently noted that a reduced order of integration may actually improve the accuracy and convergence of the approximate solution.
The author shares the opinion expressed for instance by Besseling (1983a,b) that it might be doubted whether the mathematical formulation of the physical problem in terms of functional analysis without taking due note of the underlying physical concepts is indeed a fruitful starting point. As distinct from the above viewpoint, the finite element approach presented here is perhaps best characterized as a direct discrete description of continuum thermomechanics. Rather than starting out from the standard continuum description in terms of differential field equations, we return in this approach to the most primitive of continuum thermomechanical concepts and we derive, by consideration of the thermomechanical behaviour of finite elements of the continuous medium, the discrete equations of thermomechanics.
The essential novel feature of the present approach is the concept of so-called material sampling points. This concept is believed to complete the analogy between continuum field theory and the discrete, finite element theory (Van der Giessen, 1986e).
The main line of thinking will be exemplified by considering spatial quadratic triangular elements that may be used e.g. in a spatial description of two-dimensional thermomechanical pro cesses.
3.1 Continuum models of matter i
There is a profound doubt whether a macroscopic description of actually corpuscular matter on the basis of the physical behaviour of the constituent particles, if ever successful, will bear any practical significance. Therefore, the actual material
is commonly modelled by a continuous medium, in terms of which a mathematical description of the macroscopic behaviour may then be obtained.
Let the body under consideration be replaced with a body $ of a continuous medium, being a connected compact set (or manifold) of elements X, the material points of the continuous medium. The configuration of the body at some instant t is defined by the bijection fè—»^ : x=x(X,t) or X=X(x,t), xe^, and
can be represented by a domain Bc^. The local configuration is defined as the infinitesimal neighbourhood Nx of the point xeS corresponding to the infinitesimal volume element Jf^cÈ around the
material point X.
The motion of any material point X with position x(x(X,t))eX is described by the velocity vector v=3x/8t (denote the vector space of all velocities at that point by V~X) . In a material description of the behaviour of the continuum, all physical
quantities are referred to the material points X, whereas in a
spatial description they are referred to the spatial points
x=x(X,t). The material rate of change q of any physical quantity q=q(x,t) at XE? in a spatial description is given by
3q
q = —+(Vq)v. (3.1) 3t
3.2 A d i s c r e t e formulation
3 . 2 . 1 Geometry
The second step in constructing the finite element model of the actual body is to replace the continuum body Si with a
connected set of finite elements of the continuum. One can distinguish between material and spatial finite element descrip
tions. In a material description, each finite element is a distinct material body $ec$ that is being transported through space (the current configuration of $e is a (closed) region Bec^). In a spatial description, each element is a fixed (closed) region BgC^ and the continuous medium is being transported through the elements.
The characteristic feature of finite elements is that they have a relatively simple geometry which may be characterized by a small number of parameters. An important class of element geometries may be defined by a regular transformation
x = Tc(^)xc , c = i,. .,ce (3.2)
between the position vector x of any point xeBe and its coordinates f1(i=i,..,p) in a suitably chosen local coordinate
system in Be based upon polynomials Tc of typically low degree. (The vector space of all polynomials on p-dimensional space of degree up to and including d is denoted by Pp d. Its dimension is given by ( p + d ) l dim Pp d = p ! d ! p + d^ I d •) (3.3)
The xc in (3.2) are the position vectors of the ce geometric
nodal points xceBe which serve as basis points for the polynomial
transformation. If it is required that the transformation from if1
to x is governed by a complete polynomial of some order d (called
the geometric degree of the element), then the number of
geometric nodes must be chosen according to
ce = dim p Pd. (3.4)
Both the polynomials Tc and the location of the geometric nodes
oh the element must satisfy certain conditions in order to ensure continuity of the transformations (3.2) at interelement bound aries. Usually these compatibility conditions imply that the geometric nodes must be specified as boundary points common with the adjacent elements. The class of element geometries in y2
generated by (3.2) comprises e.g. of the popular triangular and quadrilateral elements both with straight or (quadratically) curved edges.
Thus, the geometry of any element for some chosen set {rc ; c = i,..,ce} is actually determined by the set {xc} of
geometric nodes; we will write Be A= { xc} . Notice that in a
material description xc=x(Xc,t) with Xc the corresponding
material geometric node in Be. The union
E
BA := U Be A
e = l
over all E elements is the discrete configuration of the body.
The subset of the nodes defining the discrete boundary of the
body is denoted by 9BA. Its complement BA/3BA contains all
internal geometric nodes.
The so-called simplex element constitutes the discrete
analogue of the infinitesimal volume element (cf. with Cauchy's famous tetrahedron in classical treatises on stresses). A simplex in yp is defined as the convex hull of p+i points xv (the
vertices of the simplex) in ^p. For simplexes in y? it is
convenient to use as local coordinates, the P + I natural coordinates £* defined so as to have the properties
P+ 1 .
f1( xw) = 81 v / J r1 = i ; v = i , . . ,P.
i = 1
The geometry of simplex elements is determined by a (complete) linear transformation (3.2) based upon these natural coordinates
(notice that indeed dim Pp 1= p + i ) . On the basis of simplex
elements, so-called complex elements of geometric degree d may be
designed by using complete polynomial transformations of order d in terms of £* . The number of geometric nodes for complex elements are hence given by (3.4). Figure 2 shows a typical simplex element and complex element for d=2 in two-dimensional space.
Complex elements (or simplex elements in the special case d=i) possess a number of beautiful geometric properties which are a result of the intrinsic symmetry properties of the underlying
(a)
"
2(b)
0 = geometrie node
Fig. 2. Two-dimensional simplex element (a) and complex element of degree two (b).
simplex regions. These symmetry properties can be formulated as follows: The geometry of a simplex in 5^p is left unchanged by any of the (p+i)I permutations of the natural coordinates (£')
(i = l/..,P+l) . For example, a simplex in 92 has 3! symmetries, the
3 even (cyclic) permutations corresponding to "rotations" in the plane and the 3 odd permutations corresponding to "reflections" with respect to the lines !2=(3, £3={1 and £l=S2 (this is the
reason for having located the geometric nodes x4,x , x6 for the complex element in Fig. 2.b as shown).
In addition to their geometric elegance, complex elements have the practical advantage that any region can be covered to any desired degree of accuracy. Furthermore they turn out to be most convenient in automatic mesh generation procedures.
3.2.2 Spatial interactions: Flow of mass and entropy
The two primitive concepts in the macroscopic description of thermomechanical processes are the flow of mass and the flow of entropy. The flow of the continuous medium through physical space is governed by the velocity field v(x,t)eV. In a similar way the flow of entropy through the< moving medium is governed by the entropy flux field h(x,t)etf (H denotes the vector space of
entropy velocities at some point x) . The interactions of any subbody with its surrounding are represented by the fields of velocity and entropy flux along with their energetically dual
fields.
As the flow of matter and the flow of entropy are fully similar processes, it is natural to set up a similar description of both processes simultaneously.
The essence of the finite element theory is now contained in the observation that knowledge of the flow of mass and of entropy is adequate in a phenomenological sense if it is available in a finite number of points of the continuum. An approximation of the
velocity and entropy flux fields may be obtained by interpolating locally, i.e. within each element, between the values in these points.
Within each element we select say ne of such nodal points xn, in each one of which we define a velocity vector vn and an entropy flux vector hn. Approximations of the fields within the element are most conveniently obtained by utilizing polynomials (of typically low degree) in terms of the local coordinates if1 as interpolation functions,
v(x):= (Dn(^)vn , h(x):= <Dn(^)hn. (3.5)
In order to meet the fundamental continuum requirements of continuity of flow of matter and of entropy, continuity of v(x) and h(x) is to be ensured at interelement boundaries. Like in the description of the element geometry (Section 3.2.1), this condition usually implies that the nodal points must be specified as boundary points common with the adjacent elements. It is natural to let c e of the nodal points coincide with the geometric nodes. A so-called isoparametric formulation is obtained by taking ne: = ce and <J)n (£ *) :=rn (£ ') (n=i,..,ne). In the case of a subparametric formulation, ne>ce and. for ^"(f1) higher order polynomials are adopted than for rc(£i). (It is well-known (e.g. Bathe, 1982) that in general a superparametric formulation, in which case ne<ce, does not yield convergent solutions upon successive mesh refinements.)
For complex elements in ^p, complete polynomial approxima tions of degree d for v(x) and h(x) are obtained by taking
ne = dim Pp d
(cf. Eq.(3.4)). In case of a subparametric formulation, the precise location of the ne-ce nodes on the element is easily established by invoking its geometric symmetry properties. Figure 3 shows a spatial quadratic simplex element for the spatial description of two-dimensional thermomechanical processes. The element geometry is defined by its ce=p+i=3 vertices (see Fig. 2.a). For a quadratic interpolation for velocity and entropy flux,
[4>n]* = [f1(2r1-l) f2(2f2-l) f3(2f3-l) 4f1f2 4f2^3 4r3f1], the (subparametric) element must contain ne=dim P22=6 nodal points. For continuity and symmetry, these nodes are located at the three vertices and at the midpoints of the three edges.
Denoting the vector spaces of velocities and entropy fluxes at nodal point n by Vn and Hn respectively, we can define the element velocity space Ve and the element entropy flux space He
g e o m e t r i e nodal point nodal point
Fig. 3. Two-dimensional spatial quadratic simplex element.
Vc
n= 1
He : - ffi Hn.
n= l
( 3 . 6 )
The velocity space VA and the entropy flux space HA for the
entire finite element model BA can be constructed either through a sum over all E elements,
:= J V, • H* ■= I H e ' (3.7)
e = l e = l
or immediately through a direct sum over all M nodal points
Vf := 0 Vr n = l N
e
n = l HA :- © Hn.The rate of exchange of energy of an element with its surroundings (either any adjacent elements or the surroundings of the body) accompanying the flow of mass and of entropy, is composed of the mechanical power WE e of the forces<exerted on the element boundary and the rate of heat supply Qe through the boundary. For each element, this rate of exchange of energy now defines the element vector spaces of nodal forces Fe and of nodal
contact temperatures Te as the dual spaces of Ve and He
respectively, ^ e = ^ e Te - He / 1=» t h r o u g h t h e s c a l a r p r o d u c t s ( i n W=Js"1=Nms~1) WEe = < te, ve> Qe = < Te, he> , VeE Ue, , hee He, tee Fe TeE Te. ( 3 . 8 )
Dual to the construction (3.6) of the element spaces from the nodal spaces, the vector spaces Fe and Te can be decomposed into the corresponding nodal vector spaces Fn and Tn respectively according to
Fe = e' Fn , Te = e" Tn. (3.9)
n = l n = 1
We note that F = Vn* and Tn=Hn*, and that both vector spaces are isomorphic to X .
The corresponding vector spaces for the entire model are defined similar to Eqs. (3.7):
E E FA := 2. Fe i TA : = 2 Te'
e = 1 e = 1
For the discrete body BA as a whole then, the externally supplied mechanical work WE and rate of heat Q are given by
E WE , .= ^ WE e = < t , V > e = l E Q := I Qe = <T'h> e = 1 / / ve\^A , hsHA , t e FA; T Ê TA.
At the boundary 3BA we have to specify either the nodal velocities {nodal entropy fluxes} or the nodal forces {contact temperatures}. Therefore we split VA {HA} according to
VA = VA° © VAC , HA = HA° ® HA'
into a subspace of prescribed ('°') velocities {entropy fluxes} and its direct complement of calculable ('c') velocities {entropy fluxes}. The subspace FA°cFA {TA°c!TA} dual to VA° {HA0} is the vector space of nodal reaction forces {reaction temperatures}. The direct complement of FA° {and TA 0} ,
FA = FA° 9 FAC , TA = TA° 9 TA C C i m C
is the subspace FA C { TA c} o f prescribed nodal forces {contact temperatures} dual to VA {HA C}.
3.2.3 The local state; the concept of material sampling points Under the influence of the interactions with the surrounding elements, the thermodynamic state of the continuous medium contained within an element will be subject to continuous change during the process. As usual we assume the existence of a local thermodynamic state: Any infinitesimal volume element of the
continuum is a distinct thermodynamic system, the state of which is determined by local state variables which refer to this
element only.
Similar to the assertion that knowledge of the flow of mass and of entropy is sufficient in a phenomenological sense if it is available in finite number of points, we now assert that it is also sufficient to have the information on the thermodynamic state available in a finite number of points. These points may be regarded as points in which the material and its state is locally sampled, for which reason we will refer to them as material sampling points. In a sense, the material sampling points in the
discrete description based upon this finite resolution approach
replace the familiar material points in continuum theory.
In addition to its nodal points, we now attribute to each finite element a number, say se, of such material sampling points xs, in each one of which we specify the pertinent local variables. (Typically, in a spatial description, spatially fixed points will be specified as sampling points, while in a material description, material points Xs will be used such that xs=x(Xs,t)). Denoting with qs a typical local state variable as an element of some vector space Qs at sampling point x8, the thermodynamic state of the element is represented by the set of vectors
se
qe e 0e := e Qs. (3.10)
s= l
Similarly, the thermodynamic state of the discrete body BA as a whole, is represented by the set of vectors
E
qA e QA ' = © 0e- (3.11)
e = 1
If desired an approximation of the field q=q(x,t) inside the element may in turn be obtained again by interpolation between the sampling point values qs (s=i,..,se).
To each material sampling point we furthermore attribute a
volume fraction vs such that the volume Ve of the element is at any instant given by
*e
Ve
=
I
vs.s= 1
The value of some quantity Qe for any element on the basis of the corresponding specific (i.e. per unit volume) sampling point quantities can now be evaluated simply as
se
Qe = I *ss"sv i
s= 1
analogous to the expression ,fqdV in continuum field theory. For instance, the mass of an element is given by
me := ^ Ps»s' (3.12) s = l
with p8 the mass density at sampling point xs.
It is essential to realize that the concept of material sampling points is fundamentally different from that of nodal points. In the nodal points, external (typically vectorial)
quantities are defined that represent the interactions between adjacent elements. In the material sampling points on the other hand, internal quantities are defined, such as mass density and
specific entropy, which characterize the local state of the material within the element. Being physically entirely different entities, nodal points and material sampling points will play an entirely different role in this finite element approach.
3.2.3.1 Specification of material sampling points
Given an element with its nodal points, the question arises of how many sampling points are to be applied per element and where they are to be located. As far as the number of sampling points is concerned, it should be realized that there is no point in adopting a description of the state of higher accuracy than that induced through the discrete description of the flow of mass and of entropy. On the other hand however, all possible local changes in mass and entropy induced by the nodal velocities and nodal entropy fluxes must bé represented. Now, the local changes of mass density and entropy are primarily determined by the local divergence of the velocity and the entropy flux. The divergence of velocity, Asv, and entropy flux, Ash, at any sampling point xs can be obtained simply from the field approximations (3.5):
Asv := tr(V„ v(x)) , Ash := tr(Vx h(x)). (3.13)
s ■- s
Denoting the vector spaces of velocity divergence and entropy divergence at sampling point x^ by DVS and DHS respectively
(notice that DVszDHszR), the expressions (3.13) define maps
Ve -» DVS : ve .—> Asv:=Aseve
Ase : . (3.14)
He -> DHS : he ^-* Ash:=Asehe.
(Notice that we use the same symbol As e for both maps, since Ve and HB are isomorphic.) Combining the divergence spaces for all
sampling points in the element through direct sums according to (3.10),
S S
DVe := ©e DVS , DHe : = ®* DHs, (3.15)
s = l I s = l
Ve —> DVe : ve i—» Aev:=Aeve
Ae : . (3.16)
H9 -» Me : he ►—* Aeh:=Aehe
from the element velocity {entropy flux} space onto the element velocity {entropy} divergence space. With dim We= d i m DHe=se it
follows from the general properties (2.1.a) that the number of sampling points should be chosen such that
se > dim In Ae (3.17)
Now, dim Im Ae is the number of independent divergence parameters determined by the pne nodal velocities {entropy fluxes}. The motion of an element in yp characterized by the nodal velocities
include p rigid body translations and Ip(p-i) rigid body rotations which do not contribute to the velocity divergence
(i.e. local volume changes); similar observations apply to the entropy fluxes (though the interpretation may not be so obvious). Further we observe that in yp, one divergence parameter is defined by lp(p+i) quantities (notice that this is equal to the dimension of a symmetric second-order tensor space in yp) . (This statement may be demonstrated by considering volume changes: An elementary volume element (a simplex!) in ^p is defined by p+i points and its volume changes are defined by the changes in
length of its (p+l)!/2!(p_i)!=fp(p+I) edges.) Hence we conclude that the number of sampling points must satisfy
p np- p - | p ( p - l > 2n
se > — S - _ = _ i _ 1 . ( 3.1 8 ) ip(p+i) p+i
Again we emphasize that there is no point in adopting a more detailed description of the state of the material within the element than is compatible with the description of the spatial interactions. That is, one should aim at a minimal number of sampling points satisfying the lower bound in (3.18). For example, the minimum number of sampling points that satisfies condition (3.18) for a quadratic, isoparametric quadrilateral element (ne=ce=8) is se=5.
In general, it is not possible to specify precisely the position of the material sampling points within an element together with the associated volume fractions. The only guideline available is the aim to optimally sample the material contained within the element. With this in view, the conditions one can put
forward in general regarding the position of the sampling points and the size of the volume fractions are:
a. Each sampling point is located in the interior of the
element.
b. The position of the sampling points respects the
c . For each e l e m e n t , t h e i d e n t i t y
Ps"s = S Pd V' * e
s = l Be
where p=p(x) is the approximation of the mass density field given by p(x)=ps(£*)ps with 9s(fc1) denoting the
interpolation functions for the point quantities, must hold for all ps (s=i,..,se).
The motivation for the first requirement lies in the very concept of material sampling points: In each sampling point, the material contained within the element is locally sampled; obviously, the collection of local states at the sampling points cannot properly represent the state of the element as a whole if the sampling points are located on the element boundary. (It is noted that this implies that field approximations on the basis of sampling point quantities will in general not be continuous at inter-element boundaries.) The second condition expresses the fact that in sampling the material, we do not wish to bias one sampling point over the other. The third condition is connected with the two main steps in constructing a finite element model of a body: The first step is to replace the actual corpuscular material with a continuous medium and in such a way that each infinitesimal volume element of the continuum is attributed a mass density equal to the macroscopic mass per unit volume of the actual material. This ensures an identical macroscopic motion of mass of the actual material and the continuum model. The third requirement now is the analogue of this condition in the second step of replacing the continuous body by a finite element assembly and ensures identical element mass according to the continuum and the finite element model.
3.2.3.2 Material sampling points in complex elements
Owing to the intrinsic properties of the underlying simplex elements, the conditions for the number of sampling points and their location in complex elements in general can be rendered more specific.
First of all, it is observed that in a dth-order complex
element, i.e. velocity and entropy flux are interpolated using complete polynomials of degree d, the number of independent divergence parameter contained in the interpolation functions (3.5) is exactly equal to the dimension of the polynomial space of degree d-i. Hence, according to (3.3), the minimum number of material sampling points in a complex element of order d satisfying (3.17) is equal to
se =
'p + d- r
direct substitution for p=i,2,3 successively.)
Secondly, in the case of complex elements, the third requirement stated in the previous section for the location of the material sampling points and their volume fractions can be replaced with a stronger condition. Defining first the element
inertia tensors of order k for the continuum model (superscript <p) and for the finite element model (superscript A) by
Ik*> := ƒ px ®- '®x dV B„ k t i m e s T A -"-k S= 1 Pss-"-s x »«.®x8 Vs k t i m e s
(zeroth-order inertia is the mass of the element) , the third requirement reads
cc. For any complex element of order d, the identity V p8
IkA = IkP
must hold for all ke[o,d-i], where p=p(x) is the approximation of the mass density field using complete polynomial interpolation functions of degree d-i for the sampling point quantities.
This condition, as a generalization of the original general one, accounts for the finite dimensions of a finite element (as opposed to the infinitesimal volume elements in continuum theory, in which case inertia of order larger than zero is not defined) . It should be realized that the present condition may be formulated only for complex elements: The reason for being able to exactly specify the upper bound d-i for the order of the inertia tensors is a consequence of the properties of the
under-O = geometric nodal point • =nodal point
X i m m a t e r i a l s a m p l i n g point
Fig. 4. The material sampling points for a spatial
quadratic simplex element.
lying simplex elements and of complete polynomials (in fact, both properties rely upon Pascal's triangle).
Now, the assertion made here is that the three requirements a, b and cc are sufficient to establish for any complex element the position of the se sampling points (pse sampling point coordinates) and the associated volume fractions. As an example,
Van der Giessen (1986e) has discussed the procedure for two-dimensional spatial quadratic simplex elements (see Fig. 3 ) ; the result is shown in Fig. 4. (An explicit expression for the corresponding map Ae, Eq. (3.16), is given in the Appendix.) All attempts to prove this assertion in general however have failed sofar.
3.2.4 Determinism in finite element thermomechanics
Throughout this work, we shall adopt a strict deterministic
approach to the thermomechanics of continua. This approach is in the spirit of the deterministic approach in classical dynamics, where, for any given initial position and any given initial velocity of a mass point, the subsequent motion is determined
from physical laws. For the thermomechanical behaviour of continua we adopt here a deterministic standpoint which goes beyond the deterministic laws of pure continuum dynamics as it concerns also continuum thermodynamics.
Besseling (1968) and Van der Giessen (1986a,b) have formulated a principle of determinism within the framework of a
field description of continuum thermomechanics. Here we postulate this principle in a form adapted to the present finite element formulation:
An initial velocity distribution veeV? and an initial entropy flux distribution heeHe for a finite element in some given configuration Be A, together with an initial thermo-dynamic state of the medium characterized by {qeeQe}/ prescribed independently of each other at some instant, shall suffice to determine the subsequent motion of the element and the accompanying local thermodynamic processes at its material sampling points.
It should be noted that in addition to the assumption that the local state at each sampling point is determined by the instan taneous values of the local state variables (cf. Sec. 3.2.3), this principle of determinism now asserts that the instantaneous rates of change of all local state variables at each sampling point are completely determined by these local state variables and the instantaneous element motion and entropy flux.
The validity of such a strict deterministic approach may be the subject of controversy, in particular in so far as irreversible thermodynamic processes are concerned. However, it should be observed that in fact determinism is at the basis of most theories of continuum thermomechanics (including the work of two great contributors in this field, viz. Ziegler (e.g. 1977) and Biot (e.g. 1977)) and that until this moment no evidence has been found for its falsity.
3.2.5 The discrete equations of thermomechanics
3.2.5.1 Material derivatives in spatial descriptions
According to (3.1), the material rate of change of any sampling point quantity qs in a spatial description is given by the sum of the local rate of change 3qs/3t at x8 and a convective contribution. This convective contribution is determined by the velocity vs at the sampling point xs and the local gradient of qs
(denoted symbolically by Vsq ) . Hence, the discrete analogue of (3.1) reads
qs = T~ + (VsS> vs. (3.20)
3t
The sampling point velocity v? can be evaluated by interpolation between the nodal velocities, Eq. (3.5.a). Denoting the sampling point velocity space by SVS (5VS~X), the inter polation defines a map
Bs e : Ve -» SVé : ve i—» vs:=Bseve. (3.21) Combining the velocity spaces for all sampling points in the
element by a direct sum of the form (3.10), we obtain for each element a map
Be : Ve _» SVe : ve —-* Ve:=Beve (3.22) from the nodal velocity space into the sampling point velocity
space SVe. This map for two-dimensional quadratic simplex elements is given in the Appendix.
For the evaluation of the local gradient Vsq at each sampling point, various procedures may be used. The first (and most obvious) procedure would be to calculate the gradient of the
field approximation within the element obtained from the sampling point values. However, this approach is not applicable when linear simplex elements are used since, according to (3.19), these elements contain one sampling point only. Furthermore, poor approximations of the gradient may result: For instance, when using quadratic simplex elements, the same value for the gradient will be obtained for all three sampling points in the element, since these points define exactly a linear distribution.
In the case of a linear simplex element, the gradient at its single (see Eq. (3.19)) sampling point is evaluated in a straightforward manner from the sampling point values in the p + i adjacent elements (a linear distribution in ^p is defined by the values in p+i points). This procedure, which seems more in the spirit of the present finite resolution approach, may be generalized in order to apply also to higher order complex elements. The key step is to select p+i points in the neighbourhood of the sampling point under consideration which span a simplex region containing this sampling point; the local gradient is then evaluated in the same way as in case of linear simplex elements. Figure 5 illustrates the procedure for
quad-X r material sampling point
□ =radditional basis point
Fig. 5. Evaluation of gradient at material sampling point o
from values at i, ii and iii.
ratic simplex elements. Let o be the sampling point at which to evaluate the gradient. In addition to the neighbouring sampling points i and 11 in the two adjacent elements, we specify a third basis point at the centroid of the element. The value at this point H I is obtained from a simple interpolation between the sampling point values in the element and the gradient at o is evaluated from the values at i, I.I and 111.
It should be realized that this procedure implicates that the gradient Vsq cannot be computed at the element level. Instead it is obtained in general through a linear map VS A from the system's vector space QL, Eq, (3.11), for the entire model:
VsS := VSA<IA , qA£ÖA- (3.23)
Obviously, several other procedures for the evaluation of the sampling point gradient may be designed. A similar issue arises in Galerkin-type methods using numerical integration and an extensive literature on this subject is already available (a recent contribution is due to Huetink (1986)).
3.2.5.2 Deformations of finite elements
At each point x of the continuum the rate of change dx'eV of a line element dxeX in a neighbourhood Nx is governed by a local linear map
L : X —> V : dx ■—* dx' :=L dx
characterized by the second order rate of deformation tensor L at
L := Vxv (3.24)
as the gradient of the velocity field. The vector space L : = Lin(X',V) of local deformation-rates is isomorphic to the (11)
tensor space X®X*.
In the finite element model, now consider at each material sampling point xs an infinitesimal neighbourhood Ns and an
arbitrary line element dxs in Ns. The rate of change of dxs is
then governed by a local rate of deformation tensor LseLs
according to
dxs':=Ls dxs.
The strain-rate ds and the spin os at the sampling point are
defined as usual by
ds: = |(Ls+Ls T) , us: = H Ls- Ls T) (3.25)
and are elements of the subspaces D$ and IVs of Ls satisfying
LS=DS®WS. By local application of (3.24) at each sampling point
on the basis of the velocity field approximation within the element according to (3.5.a), the sampling point rates of deformations are given by maps
Gs e : Ve -> Ls : ve ►—* Ls= Gs eve. (3.26)
These maps then define the maps
Ds e : Ve -» Ds : ve I—* ds= Ds eve
onto the sampling point strain-rate space.
Combining the vector spaces Ls and Ds for all sampling
points in an element through direct sums,
se se
Le := 9 Ls , De := 8 Ds, (3.27)
s = l s = l
we obtain for each element the following maps
(3.28) Ge : ^e -> Le : ve ►—» Le=Geve
De : Ve -» De • ve •—» de= Deve >
onto the element rate of deformation space Le and the element
strain-rate space Pe respectively. The vectors (in an abstract
algebraic sense) Le e ie an<* dee öe represent the rates of
deformation and strain respectively of each finite element. From the general properties (2.1) of linear maps it follows that for the maps in (3.28)
dim Le > dim Im Ge = dim Ve - dim Ker Ge, )
(3.29) Since the rigid body translations of the element are defined by
Geve=0 or veeXer Ge,
while all rigid body motions (translations and rotations) are given by
Deve=0 or veeJTer De, (3.30)
the inequalities (3.29) ensure that all possible rates of deformation and of strain of the element with velocity space Ve can be characterized by vectors in Le and De respectively. By observing that for finite elements in ^p we have
dim Le = p2se , dim Ker Ge = p,
dim De = |ptp+i)se , dim Ker De = èp(p+i>,
and dim Ve-pne, it readily follows that if se is taken according to (3.18), both inequalities (3.29) are indeed satisfied (the proof of (3.29.a) relies on the trivial fact that ne>p>l). (In fact, the inequalities (3.29.b) and (3.18) are identical, by virtue of the fact that Asv=tr Ls=tr ds.) Explicit expressions
for the maps Ge and De for two-dimensional quadratic simplex elements as shown in Figs. 3 and 4 are given in the Appendix.
3.2.5.3 The balances of mass and of entropy
In addition to the local mass density ps introduced already in Sec. 3.2.3, we define at each material sampling point the local entropy ss per unit mass. Their material rates of change due to the flow of mass and the flow of entropy are given by the well-known expressions
p8 := -p8(Asv) , psss := -(Ash) + ps;>8 (3.31) where Ó S is the local rate of entropy production per unit mass
due to the accompanying irreversible processes and which must satisfy at any instant the Second Law of thermodynamics
38 > 0 , Vs. (3.32)
Using then Eqs. (3.14), we find for each element
Ps : = - P sAs eve > PSSS := - As ehe+ psJs, (3.33) which would be referred to as the discrete balances of mass and
of entropy (although in the present formulation, they appear as evolution equations for mass density and entropy). The material rates p8 and s8 in a spatial description are given by expressions of the form (3.20-22). Conservation of mass for each element, i.e. me=0 with me given by (3.12), irrespective of the current mass distribution demands that the rate of change of the sampling point volume fractions satisfies
^s = »sAseve < Vs. (3.34)
3.2.5.4 The principles of virtual power and heat flow
The physical laws governing the motion of an element and (its thermodynamic counterpart) the heat flow through the element find expression in the principles of virtual power and of virtual heat flow (Besseling, 1979, 1983b).
In harmony with the principle of determinism put forward in Sec. 3.2.4, let us consider some element Be A of the body in some configuration BA. At that particular instant let us apply instantaneously an arbitrary, virtual velocity distribution
8veeVe and, independently of this, an arbitrary virtual entropy flux distribution 8hee#e (here and in the sequel, the symbol 8 merely indicates that we are dealing with virtual and, by definition, time independent quantities).
The principle of virtual work now states as a necessary
condition for any motion of the element that the virtual power 8We of the externally applied forces on the element be equal to the virtual rate of change of the kinetic energy, 8Te for all virtual rigid body motions. The virtual power comprises of i) the virtual power 8W'e of the forces at the element boundary due to the interactions with the adjacent elements (see Eq. (3.8.a)) and
ii) the virtual power 8W!e of the internal or body forces (e.g. due to gravity). Hence, using (3.30),
8WEe + S W ^ = 5Te , V 8veeüTer De. (3.35)
It is in the spirit of the present finite resolution approach now to sample both the internal virtual power 5W1 e and the virtual kinetic energy-rate 8Te at the material sampling points in the element. With vs denoting as before the actual velocity at some sampling point, the kinetic energy at that point is given by
Ts := l<msvs,vs> , ms: = psPs.
Since by virtue of (3.33-34) the "sampling point mass" ms is constant, the actual rate of kinetic energy is found to be given by the scalar product of the inertia force m8vg and the
instantaneous velocity vs : Ts=<msvs,vs>. With f8 denoting the body force per unit mass at the sampling point considered, the sampling point internal virtual power 8W1 s and virtual kinetic energy rate 8TS at a virtual velocity 5vs are now given by the expressions
«W1, := <Fs,8vs> , 8TS := <m8vs,8vs>
with Fs:=msfs. The material acceleration and virtual velocity at the sampling point, are obtained by interpolation between the nodal accelerations and virtual velocities (cf. Eq. (3.21)):
