- Email: [email protected]
Elastic eigenstates in finite element modelling of large anisotropic elasticity
Computer methods in applied mechanics and engineering Comput.
EISEVIER
Elastic eigenstates
Methods
Appl.
Mech. Engrg.
160 (1998) 325-335
in finite element modelling anisotropic elasticity
Pawei Dhiewski Institute of Fundamental
of large
*, Piotr Rodzik
Technological Research, Swietokrzyska 21, 00-049 Warsaw, Poland
Received
11 April
1997; accepted
28 August
1997
Abstract FE algorithm for large anisotropic hyper-elastic deformations is presented. This algorithm is based on the concept of the local reference configuration. The spectra1 decomposition of the elastic modulus tensor into eigentensors of elastic deformations has @ 1998 Elsevier Science S.A. All rights been utilized in the constitutive equations. Some numerical results are also presented. reserved.
1. Introduction
The idea of description of elastic behaviour by means to the elastic eigenstates is due to (lord Kelvin) Thomson [l]. Recently, it has been popularized and promoted by a number of researchers (see e.g. Rychlewski [2-4], Pipkin [6], Sutcliffe [5], Mehrabadi and Cowin [7,8]). The spectral decomposition of modulus tensor has many advantages. From the mathematical point of view it assumes a fundamental role in determination of the nonlinear isotropic tensor functions. Furthermore, from the physical point of view, the eigenstates play an important role in the analysis of the material instability (acoustic tensor), wave propagation and in many other problems of mechanics and physics. The apparent advantages of the eigenstates formulation prompts the present authors to undertake a more detailed investigation of the description of anisotropic hyper-elastic materials by means of the elastic eigenstates. Most of the widely used models for large elastic deformations of solids falls into a category of the hypo-elastic or isotropic hyper-elastic models. It should be emphasized that hypo-elastic models disregard the most fundamental law of physics, i.e. the energy conservation. The hypo-elastic models are postulated in incremental form and do not guarantee the existence of potential, which could be identified with energy. As a consequence, the work done on closed deformation loop may not balance. In other words, it is one more attempt to implement the idea of a perpetum mobile! On the other hand, most of the hyper-elastic models, like neo-Hookian, Mooney-Rivlin’s or Ogden’s models, are not suitable for modelling the elastic anisotropy because, by definition, they are devoted to an isotropic behaviour and, usually, they do not conserve energy when tensors of elastic moduli become anisotropic. Therefore, in the following considerations we shall implement the conserving energy hyper-elastic approach to anisotropic material. In the present paper our interest is focused on the anisotropic hyper-elastic models. In Section 2 a thermodynamic foundation of the finite thermoelasticity is discussed. Section 3 is devoted to the
*Corresponding
author.
00457825/98/$19.00 @ 1998 Elsevier PII:SOO45-7825(97)00295-8
Science
S.A. All rights
reserved.
P Diuiewski,
326
P Rodzik/Cotnput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
fundamental constitutive equations for hyper-elastic materials. Subsection 3.1 presents the foundations of the theory of constitutive modelling based on the elastic eigenstates. Section 4 is devoted to the theoretical foundation of the numerical algorithm implemented here to the finite element code. The detailed equations are presented in relation to the plain problem. In Section 5 a comparison of the obtained results for the present and neo-Hookian models are discussed. 1.1. Reference 1.1.1.
configurations
Total Lagrangian formulation
In this most widely used approach the reference configuration is identified with the Lagrangian configuration. Such an approach is frequently used to model the elastic behaviour of materials. In the case of the elastic-plastic deformations, the concept of the Lagrangian configuration does not seem to be a good idea, since a relatively large plastic deformation can lead to a strong resealing of material constants, e.g. elastic moduli. There exists a lot of methods to overcome this problem. For example, instead of the Green strain measure one can assume the logarithmic measure. 1.1.2.
Updated Lagrangian formulation
To avoid the description favouring a certain configuration many authors prefer the usage of the socalled updated Lagrangian scheme. This scheme has also many shortcomings. For example, the results obtained depend on the length of the deformation step after which the updating is made. Moreover, strains calculated by means of the updated Lagrangian scheme depend not only on the initial and current configurations but also on the deformation path (cf. [9,10])! Besides, for the closed deformation loops the strain measure updated step by step along the closed path does not vanish. It means that such updated strain measure should be treated as a state parameter rather than as a geometric measure of the configuration change. So, it should be concluded that the updated Lagrangian scheme cannot be treated as a unique constitutive modelling tool for large deformations. It is rather a numerical technique being more or less sensible constitutive modelling itself. Due to the above shortcomings of the Lagrangian formulation many attempts have been made to adapt spatial formulation to the description of hyperelastic and/or hyper-elastic-plastic deformations. Let us mention a few of them. 1.1.3. Incremental formulations
Generally, in incremental formulations the constitutive relations are postulated between the objective time derivatives of stress and the elastic strain velocity tensors. In the literature, various time derivatives are employed: e.g. the material, corrotational and Lie derivatives. The currently used approaches fall in one of the following groups: (1) The models that reduce to the total Lagrangian scheme. In this case the elastic energy balance is preserved. For example: formulations based on the Lie derivative are often nothing more as only the total Lagrangian formulation rewritten in terms of the current configuration, cf. Marsden and Hughes [14]. So that the use of the Lie derivative for anisotropic materials retains all the shortcomings of the total Lagrangian formulation. (2) The formulations that do not reduce to the total Lagrangian. In this case the elastic energy is not preserved, e.g. for the Zaremba-Jauman derivative. Recently, special attention was given to the logarithmic strain measure and its derivatives (cf. [13,11,12]). Using this relation, the time derivative of strain is replaced by a finite difference. Such an approach has also a number of shortcomings. For instance, considering a loop in the deformation space, it is easy to show that for finite deformation steps the energy conservation is violated. Obviously, decreasing the length of updating step the energy balance can be reached in an asymptotic manner. In this place it is worth recalling that the fundamental advantage of the total Lagrangian formulation is that it preserves the elastic energy regardless of the length of subsequent deformation steps. 1.1.4. Spatial formulation based on the local reference configuration The formulation implemented here is a combination of the total Lagrangian formulation for elastic deformation and spatial formulation for plastic deformations. In the case of the elastic deformation this
P Dtuiewski, P Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (19%) 325-33
321
formulation retains the properties of the total Lagrangian, while, for plastic deformations, this formulation does not have the shortcomings of the Lagrangian formulation. The idea of the local reference configuration was introduced by Teodosiu [15], Mandel [16] and Rice [17] (cf. also [IS]).
2. Free energy density For elastic as well as for many elastic-plastic materials the local form of the balance laws for the mass, momentum, moment of momentum and energy can be assumed as follows: p+pdivu=O
(1)
diva+pb-pti=O
(2)
a-UT=0
(3)
-pti+u:Vv-divq+ph=O
(4)
where p, u, b, u, u, q, h denote the mass density, Cauchy stress tensor, body force density, velocity, internal energy density, heat flux, and heat source density, respectively. The balance equations can be completed by the entropy inequality pq+div
(
f
>
ph -r”O
(5)
where n is the entropy density. Using (4), the inequality (5) can be rewritten in the following form: -p$-pv?+u:Vv-;gradT>O
(6)
where $J is the free energy density defined as I,!J= u - VT. For elastic materials as well as for many elastic-plastic materials the constitutive equation for the free energy density can be stated as follows: ICI= 1Cl(&,T)
(7)
where Ze is the elastic strain tensor given as ze = ; (FZF,
- 1)
In the case of elastic materials the elastic deformation by the following equation: F = F,Fo
(8) tensor, Fe, is related to the deformation
gradient (9)
where F. is a tensor responsible for residual stresses in the reference configuration. This tensor describes the rotation and stretch of material with respect to the reference configuration. Usually this tensor is ignored; nevertheless, in many physical problems of elasticity the global reference configuration (without residual stress) does not exist, e.g. for dislocated crystals. Using the standard methods of thermodynamics the following identification can be introduced: 77=-z
a*
u = Fea-f$-
(10) e
FzdetF;’
(II)
The last equation can be also rewritten in the form A
A a*
u=pd^Ee
(12)
where G is a counterpart of the second Piola-Kirchhoff stress tensor determined in relation to the local relaxed isoclinic configuration (cf. [15]). The mentioned configuration is also called the intermediate
P DhQewski, l? Rodzik/Comput.
328
configuration formula:
Methods Appl. Mech. Engrg. 160 (1998) 325-335
(cf. [17]). The considered stress tensor is related to the Cauchy stress tensor by the following
a=F,&FTdetF;’
3. Constitutive
(13)
model
In this paper, by a linear elastic model we mean a linear constitutive
equation
G=&:e
(14)
where &. has been defined by (8), while fi is a fourth-order tensor of elastic constants. In general, this tensor has 36 independent components. Consider the constitutive model (14) in terms of a current configuration. Then, the above equation can be rewritten as -1
(15)
where Dkl”‘* = F,“,F,‘,FrMF&,EKLMN
detF;’
It is easy to note that the transformation rule (16) can be cumbersome in practical applications. Moreover, in such a case the strain tensor as well as the product Dkfmnemn are also to be calculated. The constitutive equation (14) can also be rewritten in the incremental form I
I
1 ;=D:(m+de)
1
(18)
where (w + de) = &F,’
(19)
-
$ = ir + (o
+ d,)a
+ a(cc, + d,)= + tr(cc, + d,)a
(20)
w and d, denote, respectively, the rate of rotation and elastic stretch of material. 3.1. Elastic eigenstates
Without loss of generality the elastic modulus tensor 5 can be decomposed
into six eigenstates
where the eigen-‘vectors’ +, are second-order tensors, while six scalars Db determine the elastic moduli for six sequential eigenstates. With respect to the required positive value of the elastic energy, the matrix of the components of fi is positively determined. It leads to the real (noncomplex) six eigenvalues. In the case of elastic isotropy the five eigenstates correspond to the arbitrary chosen five mutually independent pure shearings while the sixth corresponds to the volume compression. The spectral decomposition of elastic modulus tensor has many advantages. For example, the compliance tensor being the inverse of the elastic modulus tensor is obtained directly as
6 1 6-l =~7Tiib&ib b=l Db
(22)
P Dtuiewski, II Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325~.US
Substituting (21) into (14) and next (14) into (13) we obtain the following constitutive the Cauchy stress tensor
329
equation for
(23)
where ab = F, iibFz det F,’
(24)
D, = 6b 1 &b = 2 tr(ab - iib)
(25) (26)
The tensors ab denote eigenstates deformed, while the scalars &b represents the strain of eigenstates. Using (23) only the transformation of eigenstates has to be done from which the Cauchy stress tensor is built up in simple manner. On the other hand, if the transformation rules for the eigenvalues and eigenvectors are assumed in the form of (24), (25), then it is easy to prove that the transformed eigenvectors satisfy the following canonical form: 6 Dijwvl
= c
D,
a:” urn
(27)
h=I which means that the transformation
rules (24) and (25) hold the elastic eigenstates indifferently.
4. Finite element method It is worth emphasizing that the anisotropic hyper-elastic deformations are usually considered with the use of the initial (Lagrangian) configuration. Nevertheless, FEM as a method of integration of differential equations can be employed for solving the problem stated in the arbitrary chosen configuration. Therefore let us consider the balance of momentum formulated in the current configuration (cf. (2)). The balance law leads to the following form of the principle of virtual work:
.I ”
6upi,du+
Ji3VTu
:rrdu+ubdu=&ds
(28)
”
where Su denotes the virtual displacements. displacements
Assuming that the displacements
are determined
by nodal
6u = wsu the following equation for integration
(29) over finite elements is obtained
J
Wpridu+~Wbdu+~VTW~du=~W~ds
u
(30)
where W is the weighting function. The above equation is independent of the constitutive equations on the Cauchy stress, so it is valid also for anisotropic hyper-elastic materials. Generally, for hyperelastic anisotropic materials, we can assume here that the Cauchy stress tensor depends on the elastic deformation tensor. According to (9) we find F,’
= Fo(l - Vlviui)
(31)
where Ni denotes the shape function for ith node. In the index notation the above equation reads -’ N
F en - FtM (SC - VnNi$‘g,M)
(32)
P Dluiewski,
330
El Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
where urn is the vector of the nodal displacement determined in the direction of mth curvilinear coordinate, while g,M denotes the shifters between the vector components obtained in two mutually different vector basis in the Euclidean space (cf. [20]). For determination of the tangent stiffness matrix the important role takes the derivative aF,/aUi, where ui denotes the vector of displacements in the ith node. The derivative can be determined on the basis of the well-known formula
-
=
dUi
According to (32) we find
where NIUY
Un =
(35)
(cf. (31)). In our case the quasistatic problem is under consideration, and body forces are omitted, i.e. pii
where, by assumption, the inertial
=o
(36)
b=O
(37)
In such a case, equation set (30) can be considered as a nonlinear set P(u)
=.f
(38)
where P =
J’
VTWadv
(39)
f=jW& au
(40)
This set can be solved by means of the Newton-Raphson method. Then tangent stiffness matrix determined as K = dP/du holds the following general relation: K=
.I
VTW
a(~ detF,) det F;’ dv
(41)
dU
V
In the particular case, when the elastic behaviour is governed by (14), we find d
(a’j det Fe) auf
detF;’
=
FL@”
+ @“FL
+ F~g~,D’imn
>
[grV,Ni
+ U”VkVINi]
(42)
where the elastic modulus tensor Dilk’ is determined by (16). In many physical problems it is assumed that the stress tensor depends not only on the elastic deformation tensor but also on many other state variables, e.g. on the plastic deformation tensor, defect densities, plastic curvature, temperature, etc. The finite element formulation for elasticity is treated then as the point of departure for considerations on more complex systems. Therefore, in some cases, it is useful to treat the dependence (42) as a complex relation
a(a’jdetF+_) detF_l e
~F$N
6
=;N ek
g;aki + aikgi + g,,
c
D, a; LZ;“’
b&l
and (44)
P Diuiewski,
fl Rodzik/Comput.
Methods
Appl.
331
Mech. Engrg. 160 (1998) 32.5-335
The constitutive model presented above has been implemented to the FEM program FEAP originally written by R.J. Taylor. The procedure written by P. Wriggers for neo-Hookian material has been adopted by the present authors to the description of an anisotropic hyper-elastic material. In Appendix (5) the detailed equations resulting from (43), (44) are discussed for plain problems. EXAMPLE
1. In the first example the range of the solution convergency for the considered model was investigated. The plain finite element mesh has been assumed in the form of a square disk (see Fig. l(a)). Mesh contained 400 nine-node elements. The material was assumed to be isotropic (rubber) with the following elasticity constants:
mo =
I
r8.53
2.13
2.13
0.01
2.13 0.0
8.53 2.13 0.0
2.13 8.53 0.0
0.0 I ’ lo6 MPa 3.2
(45)
In the numerical calculations this matrix was decomposed a1 = [0.0,0.0,0.0, LOIT a2 = [-0.408, -0.408,0.816,0.0]T u3 = [0.707, -0.707,0.0,0.0]T u4 =
[-0.577, -0.577, -0.577,0.0]T
into four eigenvectors
and eigenvalues
D1 = 3.2. lo6 MPa Dz = 6.4. lo6 MPa D3 = 6.4. lo6 MPa D4 = 12.8. lo6 MPa
1OoooO
-
(46)
hyperelastic
anisotropic
model
- - - - neohookian model
P 80000600004OOOC
__--
2oooo 0 -2oOOo-40000. 0.5
/’ +’
I’
+ end of solution convergence 1
1.5
2
l/l, (4 Fig. 1. Example I: (a) boundary conditions and FEM compression (c); (d) elogation/loading diagram.
discretization;
(b,c) deformed
disk for maximal
tension
(b) and maximal
332
P Diuiewski, P Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
The disk has undergone the constant tension and compression, respectively (cf. Fig. l(a-c). Fig. l(b,c) shows the configurations obtained for the maximal single loading step after which the convergency still held. Fig. l(d) shows the elogation/loading diagrams obtained for the present model and for the neoHookian one. The points of the loss of convergency for the maximal loading step are marked by crosses for tension and for compression. It is seen that the disk stiffness for the present model changes larger than for the neo-Hookian one. This was the main reason why the convergence of solution was lost for the present model before the neo-Hookian one (cf. Fig. l(d)).
EXAMPLE 2.The second example shows the anisotropic response of the model. The mesh, boundary conditions and loading scheme are presented in Fig. 2(a). At the beginning, it was assumed that the material is isotropic like in Example 1. Fig. 2(b) shows predicted deformation. In the second and the third passes it is assumed that the material is reinforced in horizontal and in vertical direction, respectively,
(cf. (45)). Predicted deformations anisotropy is visible.
of disk are shown in Fig. 2(c,d). A strong influence of the material
(b)
(4 Fig. 2. Example II: (a) boundary conditions mesh for anisotropically reinforced material:
and FEM discretization; (b) deformed mesh for isotropic in horizontal direction (c) and in vertical direction (d).
material;
(c,d) deformed
I? Diuuiewski, rl RodziklComput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
333
5. Final remarks The spectral decomposition of the tensor of elastic moduli is the essential part of presented approach. Proposed model was formulated for finite element method with the use of the spatial formulation based on local reference configuration. In spite of the authors expectation, the use of the spectral decomposition does not give a significant reduction of the computing time nor the numerial code. Nevertheless, from the theoretical point of view the spectral decomposition has many advantages. Therefore, some of the present results can be useful in the modelling of other, more advanced problems where the eigenstates take the important role, e.g. the strongly nonlinear elastic behaviour, elastic-plastic instability, wave propagation, coupling problems, etc.
Appendix A For plain problems the Cauchy stress tensor u and the elastic strain tensor E, reduce into four independent components. For the orthonormal coordinate system their vectorial representation can be written as [(+I = [&I, (Tyy,a=, U”lT
[&] = [&XX, &YY, &ZZ,qT
(A.11
then the elastic modulus tensor D reduces into 16 nonzero components DXXXY D XXX,’ DXXYY DXXZZ [q
=
1
;;;;I
;g
DYYZZ DZZZZ
;g
DWX
DXYYY
DXY=
DXYXY
(A.21
I
Then the tensor D can be decomposed 4 D=xD,,ab
into four eigenstates
@ ab
(A.31
b==l
where the representation
of the eigenstates deformed have a vectorial form
(A.41 and deformed eigenstates are
(A-5)
For the plain problem the invariant of deformation 1
-lx
FeX= -detF,
FeyY
gradient can be found as
(‘4.6)
-1 F:=-&F:Y 1
_‘Y
Fex=-detF,_lY eY =
F
G4.7)
e TX
1
detF,
FX e/t’
WV (A.9)
P Dluiewski, P Rodzik/Comput.
334
Using (43) and (A.6)-(A.9)
Methods Appi. Mech. Engrg. 160 (1998) 325-335
the nonzero components
of da/S,
take the form
4
wi
eY
=eYy
det F,)
8(uxx
c =
4
det F,) aF:x
d(dy
4 = 2F;Xu”y
-
2F~xaXX + c
a(#’
det F,)
detF,)
c
4
Db a: (arFzy
det F,)
aF,; a(azz det F,)
aF:Y
-
= 4 c
D,
a: (a? F.& - arF&)
Db ar (aFFiy
= 2F,Y,aYy - 2F,XyaYY+ c
- ar Fty)
(A.17)
b&l
4 c Db a? (arFzy = b=,
(A.18)
- a: Fey)
4 c Db a: (a; F& - a~:Fz!) = b=l c
4
4 =
(A.19)
Db a: (a: F& - al;*F,&)
(A.20)
D, a? (a: Fzy
(A.21)
-$Ty)
&,
4 = F;yu*Y
- F&uYy +
c
a: (az;‘Fe’y - a: Fzy)
(A.22)
Db a: (a: F& - a? F&)
(A.23)
Db
b=l
det F,)
4 = F&dY
- F& uxx + c
WY a(u”y
(A.16)
4
TX d(+Y
(A.15)
b=l
c det F,)
Db a; (ayb F& - ayF&)
2F~xvxy + c b=1
d(uzz detF,)
a(+
(A.14)
4 =2F,;@’
d(azz det F,)
ex
(A.13)
- a; F&)
aFeY, d(gzZ
(A.12)
= b=,
aF:Y
d(@‘detF,)
(a: F& - arFz’)
4 - af F$) c Db ar (a? Fzy = b=l
aF:Y
a(#
Dbar
b=1
det F,) aF&
(A.11)
Db ar (a: Fzx - aFF&)
&I
@tFY
a(+
(A.lO)
b=l
det Fe)
d(rXX
D, a+; (a: Fzy - a: F’zy)
- 2F,X, (T’Y + c
detFe)
b=l = FX &y eX
%Y
a(UxydetFe) a%
_ FY #,’ eX
+ 4 c
Db
a: (a? F& - aTFix)
(A.24)
b=l 4 =FY
eY
uxx
-
F~yuxy + c b=I
Db a; (ayF&
- a? Fzy)
(A.25)
I? Dtuiewski,
I? Rodzik/Comput.
According to (44) the nonzero components
ac& _ - F$
Methods Appl. Mech. Engrg. 160 (1998) 32533.5
of dF,/aUi
335
are
(V,Ni Ftx, + VyNi F&)
(A.26)
-aF& _- F$ (VxNi F,& + V,Ni F&) au;
(A.27)
-
aq
= Fi (V-xNi Fzy + Vy Ni F,Yy)
(A.28)
aFJy = F; (V_xNi F,*y + V,Ni FeyV) ax;
(A.29)
dFX --LX = Fi (VxNi Fty + VyNi F,Yy)
(A.30)
dFX eY
(A.31)
au; ag
= F; (VXN; Fty + V,Ni Ffy)
aF? -..AC = Fs (VxNi F,X, + V,Ni FeyX)
auf
aFti., _ FG aup
(VxNi F& + V,Ni F&)
(A.32) (A.33)
References [l] W. Thomson (Lord Kelvin), On six principal strains of an elastic solid, Phil. Trans. R. Sot. 166 (1856) 495498. [2] J. Rychlewski, CEIIlNOSSST’TUV (mathematical structures of elastic materials) [in Russian], Inst. Mech. Probl. USSR Acad. Sci. Preprint No. 217 (1983) 113. [3] J. Rychlewski, On Hooke’s law, PMM 48 (1984) 420-435 (see also English translation, Prikl. Matem. Mekhan. 48, 303-314). [4] J. Rychlewski, Unconventional approach to linear elasticity, Arch. Mech. 47 (1995) 149-171. (51 S. Sutcliffe, Spectral decomposition of the elastic tensor, J. Appl. Mech. 59 (1992) 762-773. [6] A.C. Pipkin, Constraints in linearly elastic materials, J. Elasticity 6 (1976) 179-193. [7] M.M. Mehrabadi and S.C. Cowin, Eigentensors for linear anisotropic materials, Quart. J. Mech. Appl. Math. 43 (1990) 1541. [S] S.C. Cowin and M.M. Mehrabadi, On the structure of the linear anisotropic elastic symmetries. J. Mech. Phys. Solids 40 (1992) 1459-1472. [9] K. Heiduschke, Why for finite deformations, the updated Lagrangian formulation is obsolete, in: D.R.J. Owen and E. Ofiate. eds., Computational Plasticity: Foundations and Applications. Proc. 4th Int. Conf., Barcelona 36 April, 1995 (Pineridge Press, Swansea, UK, 1995) 2165-2176. lo] K. Heiduschke, Computational aspects of the logarithmic strain space description, Int. J. Solids Struct. 33 (1996) 747-760. 111 B. Raniecki and H.V Nguyen, Isotropic elastic-plastic solids are finite strain and arbitrary pressure, Arch. Mech. 36 (1984) 687-704. 121 R. Hill, Aspects of invariants in solid mechanics, Adv. Appl. Mech. 18 (1978) l-75. 131 A. Hoger, The stress conjugated to logarithmic strain, Int. J. Solids Struct. 23 (1987) 1645-1656. 141 J.E. Marsden and T.J.R. Hughes, Mathematical Foundations of Elasticity (Prentice-Hall, Englewood Cliffs. NJ, 1983). 151 C. Teodosiu, A dynamic theory of dislocations and its applications to the theory of the elastic-plastic continuum, in: A. Simmonds et al., eds., Fundamental Aspects of Dislocation Theory Nat. Bur. Stand. Spec. Publ. 317 II (1970) 837-876. 161 J. Mandel, Plasticite classique et viscoplasticite, Lecture Notes Int. Centre of Mech. Sci., Udine (Springer-Verlag. Berlin. 1972). 171 J.R. Rice, Inelastic constitutive relations for solids: an internal variable theory and its application to metal plasticity, J. Mech. Phys. Solids 19 (1971) 433455. 181 P Dtuiewski. On geometry and continuum thermodynamics of structural defect movement, Mech. Mater. 22 (1996) 2341. 191 P Dtuiewski and H. Anttinez, Finite element simulation of dislocation field movement, CAMES 2 (1995) 141-148. 20) J.L. Ericksen, Tensor fields, in: C. Truesdell and R.A Toupin, eds., The nonlinear Field Theories. Handbuch der Physik, Vol. III/l (ed. S. Fluge, Springer, Berlin, 1960). 211 O.C. Zienkiewicz and R.J. Taylor, The Finite Element Method, 4th edition, Vol. 2 (McGraw-Hill, London, 1991).
EISEVIER
Elastic eigenstates
Methods
Appl.
Mech. Engrg.
160 (1998) 325-335
in finite element modelling anisotropic elasticity
Pawei Dhiewski Institute of Fundamental
of large
*, Piotr Rodzik
Technological Research, Swietokrzyska 21, 00-049 Warsaw, Poland
Received
11 April
1997; accepted
28 August
1997
Abstract FE algorithm for large anisotropic hyper-elastic deformations is presented. This algorithm is based on the concept of the local reference configuration. The spectra1 decomposition of the elastic modulus tensor into eigentensors of elastic deformations has @ 1998 Elsevier Science S.A. All rights been utilized in the constitutive equations. Some numerical results are also presented. reserved.
1. Introduction
The idea of description of elastic behaviour by means to the elastic eigenstates is due to (lord Kelvin) Thomson [l]. Recently, it has been popularized and promoted by a number of researchers (see e.g. Rychlewski [2-4], Pipkin [6], Sutcliffe [5], Mehrabadi and Cowin [7,8]). The spectral decomposition of modulus tensor has many advantages. From the mathematical point of view it assumes a fundamental role in determination of the nonlinear isotropic tensor functions. Furthermore, from the physical point of view, the eigenstates play an important role in the analysis of the material instability (acoustic tensor), wave propagation and in many other problems of mechanics and physics. The apparent advantages of the eigenstates formulation prompts the present authors to undertake a more detailed investigation of the description of anisotropic hyper-elastic materials by means of the elastic eigenstates. Most of the widely used models for large elastic deformations of solids falls into a category of the hypo-elastic or isotropic hyper-elastic models. It should be emphasized that hypo-elastic models disregard the most fundamental law of physics, i.e. the energy conservation. The hypo-elastic models are postulated in incremental form and do not guarantee the existence of potential, which could be identified with energy. As a consequence, the work done on closed deformation loop may not balance. In other words, it is one more attempt to implement the idea of a perpetum mobile! On the other hand, most of the hyper-elastic models, like neo-Hookian, Mooney-Rivlin’s or Ogden’s models, are not suitable for modelling the elastic anisotropy because, by definition, they are devoted to an isotropic behaviour and, usually, they do not conserve energy when tensors of elastic moduli become anisotropic. Therefore, in the following considerations we shall implement the conserving energy hyper-elastic approach to anisotropic material. In the present paper our interest is focused on the anisotropic hyper-elastic models. In Section 2 a thermodynamic foundation of the finite thermoelasticity is discussed. Section 3 is devoted to the
*Corresponding
author.
00457825/98/$19.00 @ 1998 Elsevier PII:SOO45-7825(97)00295-8
Science
S.A. All rights
reserved.
P Diuiewski,
326
P Rodzik/Cotnput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
fundamental constitutive equations for hyper-elastic materials. Subsection 3.1 presents the foundations of the theory of constitutive modelling based on the elastic eigenstates. Section 4 is devoted to the theoretical foundation of the numerical algorithm implemented here to the finite element code. The detailed equations are presented in relation to the plain problem. In Section 5 a comparison of the obtained results for the present and neo-Hookian models are discussed. 1.1. Reference 1.1.1.
configurations
Total Lagrangian formulation
In this most widely used approach the reference configuration is identified with the Lagrangian configuration. Such an approach is frequently used to model the elastic behaviour of materials. In the case of the elastic-plastic deformations, the concept of the Lagrangian configuration does not seem to be a good idea, since a relatively large plastic deformation can lead to a strong resealing of material constants, e.g. elastic moduli. There exists a lot of methods to overcome this problem. For example, instead of the Green strain measure one can assume the logarithmic measure. 1.1.2.
Updated Lagrangian formulation
To avoid the description favouring a certain configuration many authors prefer the usage of the socalled updated Lagrangian scheme. This scheme has also many shortcomings. For example, the results obtained depend on the length of the deformation step after which the updating is made. Moreover, strains calculated by means of the updated Lagrangian scheme depend not only on the initial and current configurations but also on the deformation path (cf. [9,10])! Besides, for the closed deformation loops the strain measure updated step by step along the closed path does not vanish. It means that such updated strain measure should be treated as a state parameter rather than as a geometric measure of the configuration change. So, it should be concluded that the updated Lagrangian scheme cannot be treated as a unique constitutive modelling tool for large deformations. It is rather a numerical technique being more or less sensible constitutive modelling itself. Due to the above shortcomings of the Lagrangian formulation many attempts have been made to adapt spatial formulation to the description of hyperelastic and/or hyper-elastic-plastic deformations. Let us mention a few of them. 1.1.3. Incremental formulations
Generally, in incremental formulations the constitutive relations are postulated between the objective time derivatives of stress and the elastic strain velocity tensors. In the literature, various time derivatives are employed: e.g. the material, corrotational and Lie derivatives. The currently used approaches fall in one of the following groups: (1) The models that reduce to the total Lagrangian scheme. In this case the elastic energy balance is preserved. For example: formulations based on the Lie derivative are often nothing more as only the total Lagrangian formulation rewritten in terms of the current configuration, cf. Marsden and Hughes [14]. So that the use of the Lie derivative for anisotropic materials retains all the shortcomings of the total Lagrangian formulation. (2) The formulations that do not reduce to the total Lagrangian. In this case the elastic energy is not preserved, e.g. for the Zaremba-Jauman derivative. Recently, special attention was given to the logarithmic strain measure and its derivatives (cf. [13,11,12]). Using this relation, the time derivative of strain is replaced by a finite difference. Such an approach has also a number of shortcomings. For instance, considering a loop in the deformation space, it is easy to show that for finite deformation steps the energy conservation is violated. Obviously, decreasing the length of updating step the energy balance can be reached in an asymptotic manner. In this place it is worth recalling that the fundamental advantage of the total Lagrangian formulation is that it preserves the elastic energy regardless of the length of subsequent deformation steps. 1.1.4. Spatial formulation based on the local reference configuration The formulation implemented here is a combination of the total Lagrangian formulation for elastic deformation and spatial formulation for plastic deformations. In the case of the elastic deformation this
P Dtuiewski, P Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (19%) 325-33
321
formulation retains the properties of the total Lagrangian, while, for plastic deformations, this formulation does not have the shortcomings of the Lagrangian formulation. The idea of the local reference configuration was introduced by Teodosiu [15], Mandel [16] and Rice [17] (cf. also [IS]).
2. Free energy density For elastic as well as for many elastic-plastic materials the local form of the balance laws for the mass, momentum, moment of momentum and energy can be assumed as follows: p+pdivu=O
(1)
diva+pb-pti=O
(2)
a-UT=0
(3)
-pti+u:Vv-divq+ph=O
(4)
where p, u, b, u, u, q, h denote the mass density, Cauchy stress tensor, body force density, velocity, internal energy density, heat flux, and heat source density, respectively. The balance equations can be completed by the entropy inequality pq+div
(
f
>
ph -r”O
(5)
where n is the entropy density. Using (4), the inequality (5) can be rewritten in the following form: -p$-pv?+u:Vv-;gradT>O
(6)
where $J is the free energy density defined as I,!J= u - VT. For elastic materials as well as for many elastic-plastic materials the constitutive equation for the free energy density can be stated as follows: ICI= 1Cl(&,T)
(7)
where Ze is the elastic strain tensor given as ze = ; (FZF,
- 1)
In the case of elastic materials the elastic deformation by the following equation: F = F,Fo
(8) tensor, Fe, is related to the deformation
gradient (9)
where F. is a tensor responsible for residual stresses in the reference configuration. This tensor describes the rotation and stretch of material with respect to the reference configuration. Usually this tensor is ignored; nevertheless, in many physical problems of elasticity the global reference configuration (without residual stress) does not exist, e.g. for dislocated crystals. Using the standard methods of thermodynamics the following identification can be introduced: 77=-z
a*
u = Fea-f$-
(10) e
FzdetF;’
(II)
The last equation can be also rewritten in the form A
A a*
u=pd^Ee
(12)
where G is a counterpart of the second Piola-Kirchhoff stress tensor determined in relation to the local relaxed isoclinic configuration (cf. [15]). The mentioned configuration is also called the intermediate
P DhQewski, l? Rodzik/Comput.
328
configuration formula:
Methods Appl. Mech. Engrg. 160 (1998) 325-335
(cf. [17]). The considered stress tensor is related to the Cauchy stress tensor by the following
a=F,&FTdetF;’
3. Constitutive
(13)
model
In this paper, by a linear elastic model we mean a linear constitutive
equation
G=&:e
(14)
where &. has been defined by (8), while fi is a fourth-order tensor of elastic constants. In general, this tensor has 36 independent components. Consider the constitutive model (14) in terms of a current configuration. Then, the above equation can be rewritten as -1
(15)
where Dkl”‘* = F,“,F,‘,FrMF&,EKLMN
detF;’
It is easy to note that the transformation rule (16) can be cumbersome in practical applications. Moreover, in such a case the strain tensor as well as the product Dkfmnemn are also to be calculated. The constitutive equation (14) can also be rewritten in the incremental form I
I
1 ;=D:(m+de)
1
(18)
where (w + de) = &F,’
(19)
-
$ = ir + (o
+ d,)a
+ a(cc, + d,)= + tr(cc, + d,)a
(20)
w and d, denote, respectively, the rate of rotation and elastic stretch of material. 3.1. Elastic eigenstates
Without loss of generality the elastic modulus tensor 5 can be decomposed
into six eigenstates
where the eigen-‘vectors’ +, are second-order tensors, while six scalars Db determine the elastic moduli for six sequential eigenstates. With respect to the required positive value of the elastic energy, the matrix of the components of fi is positively determined. It leads to the real (noncomplex) six eigenvalues. In the case of elastic isotropy the five eigenstates correspond to the arbitrary chosen five mutually independent pure shearings while the sixth corresponds to the volume compression. The spectral decomposition of elastic modulus tensor has many advantages. For example, the compliance tensor being the inverse of the elastic modulus tensor is obtained directly as
6 1 6-l =~7Tiib&ib b=l Db
(22)
P Dtuiewski, II Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325~.US
Substituting (21) into (14) and next (14) into (13) we obtain the following constitutive the Cauchy stress tensor
329
equation for
(23)
where ab = F, iibFz det F,’
(24)
D, = 6b 1 &b = 2 tr(ab - iib)
(25) (26)
The tensors ab denote eigenstates deformed, while the scalars &b represents the strain of eigenstates. Using (23) only the transformation of eigenstates has to be done from which the Cauchy stress tensor is built up in simple manner. On the other hand, if the transformation rules for the eigenvalues and eigenvectors are assumed in the form of (24), (25), then it is easy to prove that the transformed eigenvectors satisfy the following canonical form: 6 Dijwvl
= c
D,
a:” urn
(27)
h=I which means that the transformation
rules (24) and (25) hold the elastic eigenstates indifferently.
4. Finite element method It is worth emphasizing that the anisotropic hyper-elastic deformations are usually considered with the use of the initial (Lagrangian) configuration. Nevertheless, FEM as a method of integration of differential equations can be employed for solving the problem stated in the arbitrary chosen configuration. Therefore let us consider the balance of momentum formulated in the current configuration (cf. (2)). The balance law leads to the following form of the principle of virtual work:
.I ”
6upi,du+
Ji3VTu
:rrdu+ubdu=&ds
(28)
”
where Su denotes the virtual displacements. displacements
Assuming that the displacements
are determined
by nodal
6u = wsu the following equation for integration
(29) over finite elements is obtained
J
Wpridu+~Wbdu+~VTW~du=~W~ds
u
(30)
where W is the weighting function. The above equation is independent of the constitutive equations on the Cauchy stress, so it is valid also for anisotropic hyper-elastic materials. Generally, for hyperelastic anisotropic materials, we can assume here that the Cauchy stress tensor depends on the elastic deformation tensor. According to (9) we find F,’
= Fo(l - Vlviui)
(31)
where Ni denotes the shape function for ith node. In the index notation the above equation reads -’ N
F en - FtM (SC - VnNi$‘g,M)
(32)
P Dluiewski,
330
El Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
where urn is the vector of the nodal displacement determined in the direction of mth curvilinear coordinate, while g,M denotes the shifters between the vector components obtained in two mutually different vector basis in the Euclidean space (cf. [20]). For determination of the tangent stiffness matrix the important role takes the derivative aF,/aUi, where ui denotes the vector of displacements in the ith node. The derivative can be determined on the basis of the well-known formula
-
=
dUi
According to (32) we find
where NIUY
Un =
(35)
(cf. (31)). In our case the quasistatic problem is under consideration, and body forces are omitted, i.e. pii
where, by assumption, the inertial
=o
(36)
b=O
(37)
In such a case, equation set (30) can be considered as a nonlinear set P(u)
=.f
(38)
where P =
J’
VTWadv
(39)
f=jW& au
(40)
This set can be solved by means of the Newton-Raphson method. Then tangent stiffness matrix determined as K = dP/du holds the following general relation: K=
.I
VTW
a(~ detF,) det F;’ dv
(41)
dU
V
In the particular case, when the elastic behaviour is governed by (14), we find d
(a’j det Fe) auf
detF;’
=
FL@”
+ @“FL
+ F~g~,D’imn
>
[grV,Ni
+ U”VkVINi]
(42)
where the elastic modulus tensor Dilk’ is determined by (16). In many physical problems it is assumed that the stress tensor depends not only on the elastic deformation tensor but also on many other state variables, e.g. on the plastic deformation tensor, defect densities, plastic curvature, temperature, etc. The finite element formulation for elasticity is treated then as the point of departure for considerations on more complex systems. Therefore, in some cases, it is useful to treat the dependence (42) as a complex relation
a(a’jdetF+_) detF_l e
~F$N
6
=;N ek
g;aki + aikgi + g,,
c
D, a; LZ;“’
b&l
and (44)
P Diuiewski,
fl Rodzik/Comput.
Methods
Appl.
331
Mech. Engrg. 160 (1998) 32.5-335
The constitutive model presented above has been implemented to the FEM program FEAP originally written by R.J. Taylor. The procedure written by P. Wriggers for neo-Hookian material has been adopted by the present authors to the description of an anisotropic hyper-elastic material. In Appendix (5) the detailed equations resulting from (43), (44) are discussed for plain problems. EXAMPLE
1. In the first example the range of the solution convergency for the considered model was investigated. The plain finite element mesh has been assumed in the form of a square disk (see Fig. l(a)). Mesh contained 400 nine-node elements. The material was assumed to be isotropic (rubber) with the following elasticity constants:
mo =
I
r8.53
2.13
2.13
0.01
2.13 0.0
8.53 2.13 0.0
2.13 8.53 0.0
0.0 I ’ lo6 MPa 3.2
(45)
In the numerical calculations this matrix was decomposed a1 = [0.0,0.0,0.0, LOIT a2 = [-0.408, -0.408,0.816,0.0]T u3 = [0.707, -0.707,0.0,0.0]T u4 =
[-0.577, -0.577, -0.577,0.0]T
into four eigenvectors
and eigenvalues
D1 = 3.2. lo6 MPa Dz = 6.4. lo6 MPa D3 = 6.4. lo6 MPa D4 = 12.8. lo6 MPa
1OoooO
-
(46)
hyperelastic
anisotropic
model
- - - - neohookian model
P 80000600004OOOC
__--
2oooo 0 -2oOOo-40000. 0.5
/’ +’
I’
+ end of solution convergence 1
1.5
2
l/l, (4 Fig. 1. Example I: (a) boundary conditions and FEM compression (c); (d) elogation/loading diagram.
discretization;
(b,c) deformed
disk for maximal
tension
(b) and maximal
332
P Diuiewski, P Rodzik/Comput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
The disk has undergone the constant tension and compression, respectively (cf. Fig. l(a-c). Fig. l(b,c) shows the configurations obtained for the maximal single loading step after which the convergency still held. Fig. l(d) shows the elogation/loading diagrams obtained for the present model and for the neoHookian one. The points of the loss of convergency for the maximal loading step are marked by crosses for tension and for compression. It is seen that the disk stiffness for the present model changes larger than for the neo-Hookian one. This was the main reason why the convergence of solution was lost for the present model before the neo-Hookian one (cf. Fig. l(d)).
EXAMPLE 2.The second example shows the anisotropic response of the model. The mesh, boundary conditions and loading scheme are presented in Fig. 2(a). At the beginning, it was assumed that the material is isotropic like in Example 1. Fig. 2(b) shows predicted deformation. In the second and the third passes it is assumed that the material is reinforced in horizontal and in vertical direction, respectively,
(cf. (45)). Predicted deformations anisotropy is visible.
of disk are shown in Fig. 2(c,d). A strong influence of the material
(b)
(4 Fig. 2. Example II: (a) boundary conditions mesh for anisotropically reinforced material:
and FEM discretization; (b) deformed mesh for isotropic in horizontal direction (c) and in vertical direction (d).
material;
(c,d) deformed
I? Diuuiewski, rl RodziklComput.
Methods Appl. Mech. Engrg. 160 (1998) 325-335
333
5. Final remarks The spectral decomposition of the tensor of elastic moduli is the essential part of presented approach. Proposed model was formulated for finite element method with the use of the spatial formulation based on local reference configuration. In spite of the authors expectation, the use of the spectral decomposition does not give a significant reduction of the computing time nor the numerial code. Nevertheless, from the theoretical point of view the spectral decomposition has many advantages. Therefore, some of the present results can be useful in the modelling of other, more advanced problems where the eigenstates take the important role, e.g. the strongly nonlinear elastic behaviour, elastic-plastic instability, wave propagation, coupling problems, etc.
Appendix A For plain problems the Cauchy stress tensor u and the elastic strain tensor E, reduce into four independent components. For the orthonormal coordinate system their vectorial representation can be written as [(+I = [&I, (Tyy,a=, U”lT
[&] = [&XX, &YY, &ZZ,qT
(A.11
then the elastic modulus tensor D reduces into 16 nonzero components DXXXY D XXX,’ DXXYY DXXZZ [q
=
1
;;;;I
;g
DYYZZ DZZZZ
;g
DWX
DXYYY
DXY=
DXYXY
(A.21
I
Then the tensor D can be decomposed 4 D=xD,,ab
into four eigenstates
@ ab
(A.31
b==l
where the representation
of the eigenstates deformed have a vectorial form
(A.41 and deformed eigenstates are
(A-5)
For the plain problem the invariant of deformation 1
-lx
FeX= -detF,
FeyY
gradient can be found as
(‘4.6)
-1 F:=-&F:Y 1
_‘Y
Fex=-detF,_lY eY =
F
G4.7)
e TX
1
detF,
FX e/t’
WV (A.9)
P Dluiewski, P Rodzik/Comput.
334
Using (43) and (A.6)-(A.9)
Methods Appi. Mech. Engrg. 160 (1998) 325-335
the nonzero components
of da/S,
take the form
4
wi
eY
=eYy
det F,)
8(uxx
c =
4
det F,) aF:x
d(dy
4 = 2F;Xu”y
-
2F~xaXX + c
a(#’
det F,)
detF,)
c
4
Db a: (arFzy
det F,)
aF,; a(azz det F,)
aF:Y
-
= 4 c
D,
a: (a? F.& - arF&)
Db ar (aFFiy
= 2F,Y,aYy - 2F,XyaYY+ c
- ar Fty)
(A.17)
b&l
4 c Db a? (arFzy = b=,
(A.18)
- a: Fey)
4 c Db a: (a; F& - a~:Fz!) = b=l c
4
4 =
(A.19)
Db a: (a: F& - al;*F,&)
(A.20)
D, a? (a: Fzy
(A.21)
-$Ty)
&,
4 = F;yu*Y
- F&uYy +
c
a: (az;‘Fe’y - a: Fzy)
(A.22)
Db a: (a: F& - a? F&)
(A.23)
Db
b=l
det F,)
4 = F&dY
- F& uxx + c
WY a(u”y
(A.16)
4
TX d(+Y
(A.15)
b=l
c det F,)
Db a; (ayb F& - ayF&)
2F~xvxy + c b=1
d(uzz detF,)
a(+
(A.14)
4 =2F,;@’
d(azz det F,)
ex
(A.13)
- a; F&)
aFeY, d(gzZ
(A.12)
= b=,
aF:Y
d(@‘detF,)
(a: F& - arFz’)
4 - af F$) c Db ar (a? Fzy = b=l
aF:Y
a(#
Dbar
b=1
det F,) aF&
(A.11)
Db ar (a: Fzx - aFF&)
&I
@tFY
a(+
(A.lO)
b=l
det Fe)
d(rXX
D, a+; (a: Fzy - a: F’zy)
- 2F,X, (T’Y + c
detFe)
b=l = FX &y eX
%Y
a(UxydetFe) a%
_ FY #,’ eX
+ 4 c
Db
a: (a? F& - aTFix)
(A.24)
b=l 4 =FY
eY
uxx
-
F~yuxy + c b=I
Db a; (ayF&
- a? Fzy)
(A.25)
I? Dtuiewski,
I? Rodzik/Comput.
According to (44) the nonzero components
ac& _ - F$
Methods Appl. Mech. Engrg. 160 (1998) 32533.5
of dF,/aUi
335
are
(V,Ni Ftx, + VyNi F&)
(A.26)
-aF& _- F$ (VxNi F,& + V,Ni F&) au;
(A.27)
-
aq
= Fi (V-xNi Fzy + Vy Ni F,Yy)
(A.28)
aFJy = F; (V_xNi F,*y + V,Ni FeyV) ax;
(A.29)
dFX --LX = Fi (VxNi Fty + VyNi F,Yy)
(A.30)
dFX eY
(A.31)
au; ag
= F; (VXN; Fty + V,Ni Ffy)
aF? -..AC = Fs (VxNi F,X, + V,Ni FeyX)
auf
aFti., _ FG aup
(VxNi F& + V,Ni F&)
(A.32) (A.33)
References [l] W. Thomson (Lord Kelvin), On six principal strains of an elastic solid, Phil. Trans. R. Sot. 166 (1856) 495498. [2] J. Rychlewski, CEIIlNOSSST’TUV (mathematical structures of elastic materials) [in Russian], Inst. Mech. Probl. USSR Acad. Sci. Preprint No. 217 (1983) 113. [3] J. Rychlewski, On Hooke’s law, PMM 48 (1984) 420-435 (see also English translation, Prikl. Matem. Mekhan. 48, 303-314). [4] J. Rychlewski, Unconventional approach to linear elasticity, Arch. Mech. 47 (1995) 149-171. (51 S. Sutcliffe, Spectral decomposition of the elastic tensor, J. Appl. Mech. 59 (1992) 762-773. [6] A.C. Pipkin, Constraints in linearly elastic materials, J. Elasticity 6 (1976) 179-193. [7] M.M. Mehrabadi and S.C. Cowin, Eigentensors for linear anisotropic materials, Quart. J. Mech. Appl. Math. 43 (1990) 1541. [S] S.C. Cowin and M.M. Mehrabadi, On the structure of the linear anisotropic elastic symmetries. J. Mech. Phys. Solids 40 (1992) 1459-1472. [9] K. Heiduschke, Why for finite deformations, the updated Lagrangian formulation is obsolete, in: D.R.J. Owen and E. Ofiate. eds., Computational Plasticity: Foundations and Applications. Proc. 4th Int. Conf., Barcelona 36 April, 1995 (Pineridge Press, Swansea, UK, 1995) 2165-2176. lo] K. Heiduschke, Computational aspects of the logarithmic strain space description, Int. J. Solids Struct. 33 (1996) 747-760. 111 B. Raniecki and H.V Nguyen, Isotropic elastic-plastic solids are finite strain and arbitrary pressure, Arch. Mech. 36 (1984) 687-704. 121 R. Hill, Aspects of invariants in solid mechanics, Adv. Appl. Mech. 18 (1978) l-75. 131 A. Hoger, The stress conjugated to logarithmic strain, Int. J. Solids Struct. 23 (1987) 1645-1656. 141 J.E. Marsden and T.J.R. Hughes, Mathematical Foundations of Elasticity (Prentice-Hall, Englewood Cliffs. NJ, 1983). 151 C. Teodosiu, A dynamic theory of dislocations and its applications to the theory of the elastic-plastic continuum, in: A. Simmonds et al., eds., Fundamental Aspects of Dislocation Theory Nat. Bur. Stand. Spec. Publ. 317 II (1970) 837-876. 161 J. Mandel, Plasticite classique et viscoplasticite, Lecture Notes Int. Centre of Mech. Sci., Udine (Springer-Verlag. Berlin. 1972). 171 J.R. Rice, Inelastic constitutive relations for solids: an internal variable theory and its application to metal plasticity, J. Mech. Phys. Solids 19 (1971) 433455. 181 P Dtuiewski. On geometry and continuum thermodynamics of structural defect movement, Mech. Mater. 22 (1996) 2341. 191 P Dtuiewski and H. Anttinez, Finite element simulation of dislocation field movement, CAMES 2 (1995) 141-148. 20) J.L. Ericksen, Tensor fields, in: C. Truesdell and R.A Toupin, eds., The nonlinear Field Theories. Handbuch der Physik, Vol. III/l (ed. S. Fluge, Springer, Berlin, 1960). 211 O.C. Zienkiewicz and R.J. Taylor, The Finite Element Method, 4th edition, Vol. 2 (McGraw-Hill, London, 1991).