- Email: [email protected]
On the force–displacement characteristics of finite elements for elasticity and related problems
Finite Elements in Analysis and Design 104 (2015) 35–40
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
On the force–displacement characteristics of finite elements for elasticity and related problems J.N. Reddy n, Arun Srinivasa Department of Mechanical Engineering, Texas A & M University, College Station, TX 77843-3123, USA
art ic l e i nf o
a b s t r a c t
Article history: Received 23 October 2014 Received in revised form 25 February 2015 Accepted 28 April 2015
The displacement based finite element models for hyperelastic materials are considered, and it is shown that —based on the invariance of the discretized strain energy to translations and rigid body rotations of the elements—the magnitude of the nodal forces can be written purely in terms of the strains along the edges composing the elements, while their directions are along the unit vectors corresponding to the edges. In other words, the nodal forces are “truss like” with one major difference: the magnitude of the forces depends upon the extensional strains of neighboring edges. In order to motivate the ideas, we first explicitly demonstrate this result for triangular finite element formulation of plane elasticity problems and then show that the result is valid for any finite element formulation of hyperelasticity. & 2015 Elsevier B.V. All rights reserved.
Keywords: Plane elasticity Triangular finite elements Extensional strains along edges Mathematical proof
1. Introduction The aim of this paper is to show a relatively simple result that has hitherto not been documented in the FEM literature, despite a long history of finite element analysis of structural problems (which were originally motivated by matrix structural analysis of trusses). The idea is stated simply thus: for hyperelasticity problems, the magnitude of the nodal forces can be written purely in terms of the strains along the edges composing the elements, while their directions are along the unit vectors corresponding to the edges, i.e., no matter what the element type, we can show that the problem can always be reduced to a truss-like problem (with the “links” of the truss being represented by straight lines connecting the nodes of each element), but with one difference—namely that the magnitude of each of the forces at a node (which are along each of the links connecting the node to other nodes, like in a truss) depends not only upon the strain in the link but also on the strains along the links connecting the other nodes. This provides a surprisingly simple picture of the finite element method but provides additional insight and new possibilities for generalization and solution of problems. In order to prove this result, we proceed in two steps: in the first step, we start with the traditional FEM formulation of the isotropic 2-D elasticity problem using the linear triangular element and show explicitly that the magnitude of the forces at the nodes of a constant strain triangular (CST) element depends only upon the extensional strains of the edges of the elements that share the node, and that the directions of the nodal forces are along the edges that have the node in common. In other words, the problem is like that of a generalization of n
Corresponding author. Tel.: þ 979 862 2417. E-mail address: [email protected] (J.N. Reddy).
http://dx.doi.org/10.1016/j.finel.2015.04.011 0168-874X/& 2015 Elsevier B.V. All rights reserved.
a “truss”, i.e. the forces are along the truss elements but their magnitudes dependent upon the strains in neighboring truss elements. We then show that further that the result is not specific to a constant strain triangle, but is true of all finite element formulations of elasticity problems (both linear and non-linear). The result follows from the fact that the discretized strain energy of the finite element representation has to (1) depend only upon nodal positions (once the shape functions are chosen) and (2) must be invariant to translations and rigid body rotations of the elements, together with well-established representation theorems using the integrity bases for isotropic functions. 2. Governing equations The governing equations for plane elasticity (plane stress or plane strain) problems are summarized here (see Reddy [1,2]). The presentation here is limited to linearized elasticity, where no distinction is made between the Cauchy stress tensor and the second Piola–Kirchhoff stress tensor and between the current coordinates x and the reference coordinates X. The domain of the elastic solid is denoted by Ω, which is an open bounded subset of R2 . The boundary of Ω is denoted by Γ ¼ ∂Ω ¼ Ω Ω, where Ω represents the closure of Ω. A typical point belonging to Ω is denoted as x. We employ the customary designations for the Sobolev spaces H s ðΩÞ and H s ðΓ Þ, where s Z0 [4,5]: Z Z H s ðΩÞ ¼ u : j Ds uj 2 dΩ o1 ; H s ðΓ Þ ¼ u : j Ds uj 2 dΓ o 1 Ω
Γ
ð1Þ
36
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
The corresponding norms are denoted as J J Ω;s and J J Γ ;s . Likewise the inner products associated with these spaces are denoted as ð; ÞΩ;s and ð; ÞΓ ;s respectively. The product spaces Hs ðΩÞ ¼ ½H s ðΩÞ2 are constructed in the usual way. The governing equations of a plane elastic continuum consist of equations of equilibrium, strain–displacement relations, and stress– strain relations. These equations are expressed in matrix form. Equations of equilibrium: DT σ þ f ¼ 0
ð2Þ
where (σ xx ; σ yy ; σ xy ) are the rectangular Cartesian components of the stress tensor σ , fx and fy are the components of the body force vector f (measured per unit volume) along the x- and y-directions, respectively, and 8 9 " # ( ) σ = > < xx > fx ∂=∂x 0 ∂=∂y T σ yy ; f¼ ; σ¼ D ¼ ð3Þ f 0 ∂=∂y ∂=∂x > > y :σ ; xy
Strain–displacement relations:
εxx ¼
∂ux ; ∂x
or
ε ¼ Du;
εyy ¼
∂uy ; ∂y
2εxy ¼
8 9 ε = > < xx > ε ¼ εyy ; > : 2ε > ; xy
∂ux ∂uy þ ∂y ∂x
( u¼
ux
ð4Þ
)
uy
ð5Þ
where (εxx ; εyy ; εxy ) are the rectangular Cartesian components of the strain tensor ε and ux and uy are the rectangular Cartesian components of the displacement vector u. Stress–strain (or constitutive) relations: 8 9 2 9 38 σ = ε = c11 c12 0 > > < xx > < xx > 7 σ yy ¼ 6 ð6Þ 4 c12 c22 0 5 εyy > > > :σ > ; 0 0 c33 : 2εxy ; xy or
σ ¼ Cε;
2
c11 6 C ¼ 4 c12 0
c12
0
3
0 7 5 c33
c22 0
ð7Þ
where cij (cji ¼ cij Þ are the elasticity (material) constants for an orthotropic medium with the material principal directions (x1 ; x2 ; x3 ) coinciding with the coordinate axes ðx; y; zÞ used to describe the problem. The cij can be expressed in terms of the engineering constants ðE1 ; E2 ; ν12 ; G12 Þ for an orthotropic material, say for plane stress problems, by c11 ¼
E1
;
c22 ¼
ð1 ν12 ν21 Þ c12 ¼ ν12 c22 ¼ ν21 c11 ;
E2
ð1 ν12 ν21 Þ c33 ¼ G12
ð8Þ
Eqs. (2)–(8) can be combined into a single equation in terms of the displacement vector u as DT CDu þ f ¼ 0
ð9Þ
Eq. (9) contains, for the two-dimensional case, two equations of equilibrium in terms of the two displacement components ux and uy.
where Bðw; uÞ is a bilinear form, F ðwÞ is a linear form, Z Bðw; uÞ ¼ h ðDwÞT CðDuÞ dx e I ZΩ wT f dx þ h wT t ds F ðwÞ ¼ h Ωe
Γe
Here h denotes the thickness of the elastic solid, ðÞT denotes the transpose of the enclosed quantity, V and W are vector spaces (e.g., appropriate subsets of the product (Sobolev) space H1 ðΩÞ ¼ H 1 ðΩÞ H 1 ðΩÞ). The quantity w δu represents the weight or test function. Unlike classical solutions that are defined unambiguously point-wise, weak solutions exist with respect to test functions and are therefore understood in the context of distributions. As a result, the weak solution (and its derivatives) is typically defined unambiguously in Ω up to a set of measure zero. 4. Finite element model The finite element model associated with Eq. (11) is obtained by restricting the solution space to a finite-dimensional subspace V h of the infinite-dimensional function space V, and the weight function to a finite-dimensional sub-space W h W. As a result, in the discrete problem we seek uh A V h such that Bðwh ; uh Þ ¼ F ðwh Þ
8 wh A W h
The strains can be computed as
ε ¼ Du ¼ DΨΔ BΔ;
σ ¼ CBΔ
To obtain the finite element model, we substitute Eq. (13) into the variational statement in Eq. (10) Z Z I e e e T e T 0 ¼ he ðδΔ ÞT BT CBΔ dx he ðδΔ ÞT Ψ f dx he ðδΔ ÞT Ψ t ds Ωe
or
The principle of virtual displacements [3] is equivalent to the following variational problem: find u A V such that
where
ð10Þ
ð14Þ
where D is defined earlier in Eq. (3) and B is a 3 2n matrix 2 ∂ψ 3 ψ2 ψn 1 0 ∂∂x 0 ⋯ ∂∂x 0 6 ∂x ∂ψ ψ2 ψn 7 1 6 7 0 ∂∂y ⋯ 0 ∂∂y B ¼ DΨ ¼ 6 0 ð15Þ 7 ∂y 4 ∂ψ ∂ψ ∂ψ ∂ψ 5 ∂ψ n ∂ψ n 1 1 2 2 ⋯ ∂y ∂x ∂y ∂x ∂y ∂x
e T
¼ ðδΔ Þ K Δ F Q
8wAW
ð12Þ
The quantity h in the definition of the subspaces V h and W h represents the average size of all the elements in a given finite element discretization. We assume that the domain Ω R2 is e discretized into a set of NE non-overlapping sub-domains Ω , h e NE called finite elements, such that Ω Ω ¼ ⋃e ¼ 1 Ω . The geometry of each element is characterized using the standard isoparametric ^ e to the physical bijective mapping from the master element Ω e element Ω . e Let u be approximated over Ω by the finite element interpolations (the element label e on the variables is omitted in the interest of brevity) ( ) ( ) ux w 1 ¼ δ ux u¼ ¼ ΨΔ; w ¼ δu ¼ ¼ ΨδΔ uy w 2 ¼ δ uy " # ψ1 0 ψ2 0 … ψn 0 Ψ¼ 0 ψ 0 ψ2 … 0 ψn 1 n oT Δ ¼ u1x u1y u2x u2y … unx uny ð13Þ
3. The principle of virtual displacements and the variational problem
Bðw; uÞ ¼ F ðwÞ
ð11Þ
e
e
e
Ωe
e
Γe
ð16Þ
Ke Δ ¼ Fe þ Q e e
Ke ¼ he
ð17Þ
Z
Z Ωe
BT CB dx;
Fe ¼ he
Ω
Ψ T f dx; e
Q e ¼ he
I Γe
ΨT t ds
ð18Þ
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
The element stiffness matrix Ke is of the order 2n 2n and the element load vector Fe and the vector of internal forces Q e is of order 2n 1, where n is the number of nodes in a Lagrange finite element. For the linear right-angled triangular (master) element (n ¼3) with horizontal side of a and vertical side of b, matrix B is given by 2 1 3 1 a 0 0 0 0 a 6 7 0 1b 0 0 0 1b 7 B ¼ DΨ ¼ 6 ð19Þ 4 5 1 1 1 1 b a 0 a b 0 Fig. 1 shows the nodal displacements and forces of the linear master triangular element. The stiffness matrix takes the form (area of the triangle, A ¼ ab=2 and h is the thickness of the body) 2c 3 þ c33 Þ 11 33 12 þ c332 ðc12 ab ca112 cab c332 cab a2 b b 6 7 6 ðc12 þ c33 Þ c33 þ c22 c12 c33 c33 c22 7 2 2 6 ab ab ab a2 a2 b b 7 6 7 c11 c11 c12 7 12 cab 0 0 abh6 6 a2 e ab 7 a2 K ¼ ð20Þ 6 7 c33 c33 33 2 6 cab ca332 0 0 7 ab a2 6 7 c33 c33 c33 6 c33 7 ab 0 0 7 2 2 6 ab b b 4 5 c12 c22 c12 c22 ab 2 0 0 2 ab b
The nodal forces are given by 2c þ c33 Þ 11 þ c332 ðc12 ab a2 b 6 6 ðc12 þ c33 Þ c33 þ c22 2 6 ab a2 b 6 c11 12 cab abh6 6 a2 ¼ 6 33 > > F 2y > 2 6 cab ca332 > > > 6 > > > > c 6 33 33 F > > 3x > cab > 6 b2 > > > 4 ; : F 3y > c12 c22 ab 2 b
ca112
33 cab
12 cab
ca332
c11 a2
c332 b 33 cab
0
0
0
c33 a2 c33 ab
c33 ab c33 2 b
c12 ab
0
0
0
38 9 12 cab > u1x > > 7> > > > 1> > > > c222 7 > > uy > > b 7 > 7> 2 > c12 7> = x 7 ab 7 2 > uy > 0 7 > > > 7> > > > u3 > > 0 7 > > 7> x > > 5> > > 3> ; : c22 u y 2 b
ð21Þ where Fix and Fiy are the horizontal (x-direction) and vertical (ydirection) forces, respectively, at node i, i ¼ 1; 2; 3 (see Fig. 2). For example, at node 1 we have abh c11 c33 1 ðc12 þ c33 Þ 1 F 1x ¼ uy þ 2 ux þ 2 2 ab a b c11 c33 c33 c12 2 u2x u2y 2 u3x u3y ab ab a b abh c11 1 c33 ¼ u u2x þ 2 u1x u3x 2 a2 x b c
i c12 1 33 uy u3y þ u1y u2y þ ab ab abh c11 c c ¼ εxx þ 12 εyy þ 33 γ xy ð22Þ 2 a a b abh c12 þ c33 1 c33 c22 1 þ 2 uy ux þ F 1y ¼ 2 2 ab a b c12 2 c33 2 c33 3 c22 3 ux 2 uy ux 2 uy ab ab a b c
c22
abh c12 1 c33 1 33 2 ux ux þ ux u3x þ 2 u1y u3y þ 2 u1y u2y ¼ 2 ab ab a b abh c12 c c εxx þ 22 εyy þ 33 γ xy 2 b b a
ð23Þ
where
εyy ¼ εð2Þ ¼
u3y u1y
;
εxx ¼ εð3Þ ¼
b u3x u1x u2y u1y þ γ xy ¼ γ 23 ¼ b a
u x3 u 3y
Side 1 Undeformed element Side 2
3
Deformed element
3
θ
b
2
ux2
1 1
u1x
u1y a
2
u 2y
x
Side 3
Fig. 1. The linear triangular element with the 6 displacement degrees of freedom.
b
9 8 F 1x > > > > > > >F > > > > 1y > > > > > > > < F 2x =
¼
y
37
u2x u1x ; a
From the strain transformation equations, we have
εn ¼ εxx sin 2 θ γ xy cos θ sin θ þ εyy cos 2 θ
ð24Þ
Fig. 2. The linear triangular element with the 6 displacement degrees of freedom.
with a sin θ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 2 a2 þ b
b cos θ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 þ b
Thus, defining the axial strain along the edges 2–3 as εð1Þ and using the above strain transformation rules, we can obtain the shear strains in terms of the three axial strains εðiÞ (i¼1, 2, 3) as a b
b a
γ xy ¼ εð3Þ þ εð2Þ
2
a2 þ b ð1Þ ε ab
ð25Þ
Therefore, we can express F1x and F1y solely in terms of the extensional (normal) strains in sides 1, 2, and 3: " !# 2 abh c11 ð3Þ c12 ð2Þ c33 a ð3Þ b ð2Þ a2 þ b ð1Þ ε þ ε þ ε þ ε ε F 1x ¼ 2 a a a b b ab " # 2 abh a2 þ b 1 c11 c33 ð3Þ ð1Þ ð2Þ c33 ε ðc12 þ c33 Þε a 2 þ 2 ε ¼ ð26Þ 2 2 a a ab b " !# 2 abh c12 ð3Þ c22 ð2Þ c33 a ð3Þ b ð2Þ a2 þb ð1Þ ε þ ε þ ε þ ε ε 2 a b b a b ab " # 2 2 abh a þ b c22 c33 ð2Þ 1 ð1Þ ð3Þ c c ε b þ ε ð þc Þ ε ¼ 33 12 33 2 2 b a2 b a2 b
F 1y ¼
Similarly, abhhc12 ð2Þ c11 ð3Þ i ε þ ε ; F 2x ¼ 2 a a
F 3x ¼
F 2y ¼
ð27Þ
i c33 hh 2 ð3Þ 2 2 a ε þ b εð2Þ ða2 þ b Þεð1Þ 2a ð28Þ
i c33 hh 2 ð3Þ 2 2 a ε þ b εð2Þ ða2 þb Þεð1Þ ; 2b
F 3y ¼
abhhc22 ð2Þ c12 ð3Þ i ε þ ε 2 b b ð29Þ
One can verify that the nodal forces are in equilibrium (i.e., F 1x þ F 2x þ F 3x ¼ 0 and F 1y þ F 2y þF 3y ¼ 0). If we define the unit vectorspalong ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the three edges as 2 e1 2 ¼ i; e1 3 ¼ j, e2 3 ¼ ð ai þ bjÞ= a2 þ b Þ (where i and j are the unit vectors along the coordinate axes, it is easy to see that the nodal forces at node 1 are trivially along the edges 1–2 and 1–3. Next, if we consider the nodal forces at node 2, we see that (using
38
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
the geometry of the triangular element) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2
a a2 þ b F 2y e2 3 F2 ¼ F 2x iþ F 2y j ¼ F 2x þ F 2y e1 2 þ b b
ð30Þ
and
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 b a2 þ b F3 ¼ F 3x iþ F 3y j ¼ F 3y þ F 3x e1 3 F 3x e2 3 a a
ð31Þ
In other words, the nodal forces at each node can be written in terms of the unit vectors along the edges. For the case of b ¼a, we have the result ah 2c33 εð1Þ ðc12 þc33 Þεð2Þ ðc11 þ c33 Þεð3Þ 2 ah 2c33 εð1Þ ðc22 þ c33 Þεð2Þ ðc12 þ c33 Þεð3Þ F 1y ¼ 2 ah c33 ah ð1Þ c12 εð2Þ þc11 εð3Þ ; F 2y ¼ 2ε εð2Þ εð3Þ F 2x ¼ 2 2 ah ah c22 εð2Þ þ c12 εð3Þ F 3x ¼ c33 2εð1Þ εð2Þ εð3Þ ; F 3y ¼ 2 2 F 1x ¼
Ψ e ðrei Þ ¼ Ψ ððrei Þn Þ;
ð32Þ
5. Strain energy function and nodal forces associated with a finite element discretization of hyperelasticity 5.1. Form of the strain energy function Having explicitly demonstrated that the nodal forces for right angled triangular elements for 2-D isotropic plane elasticity can be recast in the form of a “truss-like” problem, we now turn to the general problem of FEM formulations of hyperelasticity and prove that the result is not a special quirk of triangular elements but it follows from fundamental invariance properties of the strain energy function. Here, we need to generalize the notion of an “edge” as being “the straight line joining any two nodes of a finite element”. Thus, a four-node quadrilateral has six edges, four along the boundary and two along the diagonals. Similarly an eight-node cube will have 28 edges. Thus, we suppose that the elastic medium Ω has been h discretized into a finite element mesh Ω consisting of NN nodes and NE elements. It is evident that, for a Green elastic material, given the interpolation functions and the location of the nodes, the strain energy Ψ of the discretized domain can be determined entirely by the position vectors ri of the NN nodes, since these are the generalized coordinates, that (along with the fixed interpolation functions) uniquely determine the configuration of the deformed body. Thus, for a finite element model of a Green elastic body, we must have NE X
Ψ e ðrei Þ
ðrei Þn ¼ Q e rei þ ce
ð34Þ
ðrei Þn e
A routine calculation further reveals that the direction of the nodal force vectors is along the edges of the triangle. In other words, the FEM formulation of the plane problem in elasticity using right angled triangular elements resembles a truss formed by the nodes and edges of the elements. The nodal forces are along the truss elements (edges of the FEM formulation) while their magnitudes depend upon the extensional strains of all the neighboring truss elements. This is not a surprising result for this specific case, especially considering the special geometry chosen. However, we now show that for any finite element discretization of hyperelasticity, the above result holds.
Ψ¼
elements, i runs from 1 to 3. If we choose some other higher order element that i runs from 1 to the number of nodes per element. Thus the notation Ψe stands for the strain energy of the eth element, while the notation rei stands for the current position of the ith node of the eth element (see Fig. 3). Irrespective of how the discretized strain energy is computed, it has to meet a fundamental criterion: translations and rigid body rotations (in the case of infinitesimal deformations, it is the infinitesimal rotation) superimposed on a given configuration of any element (i.e. on any given state of strain) should not change the strain energy function . Thus the elemental strain energy Ψe must satisfy (for each element)
ð33Þ
is the position vector with superposed rigid-body where motion, Q is a proper orthogonal tensor, and ce is a constant vector. The above requirement is the counterpart of the notion of invariance under superposed rigid body motions [2] for all classical continua. We now state and prove one of the important results of this study. Theorem 1. A necessary and sufficient condition for the strain energy to meet the condition in Eq. (34) is that
Ψ e ¼ Ψ e ððdeij Þ2 Þ e dij
ð35Þ
J rei rej J
¼ is the distance between the ith and jth nodes of where the eth element (see Fig. 3). Proof. The fact that form in Eq. (35) is sufficient is a consequence of the fact that the distance between any two nodes of an element is unaffected by translations and rigid body rotations of the element. In other words, it is obvious from Eq. (34), the properties of orthogonal tensors (i.e., Q T Q ¼ I, the identity tensor and Q ¼ 7 1), and the fact that ce is a constant that 2 n n n n n n n e ¼ J rei rej J 2 ¼ rei rej rei rej dij
¼ Q rei rej Q rei rej
2
e ¼ rei rej Q T Q rei rej ¼ dij ð36Þ
y NN
26
•3
27 Ωe
2
•
1
•15
24 13 12
3•
2
e
•
•
15
1
1
2
r15 = r1e
3
4
r 12e
7
10
23 22 11
Global node numbers
Element node numbers
e¼1
where superscript “e” represents the element number, a lower case subscript denotes the node number within an element, that is, ðÞei stands for the enclosed quantity at the ith node of the eth e element Ω . Thus if we choose to use mesh of triangular
x Fig. 3. A finite element mesh of linear triangular elements showing the position of the 15th global node and its node number relative to the eth element Ωe. Thus the node, which has a global number of 15, is the 1st node of the eth element.
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
To prove necessity, we note that since Eq. (34) is valid for all c we first choose c ¼ 0 and obtain
assumed that
Ψ e ððrei Þn Þ ¼ Ψ e ðQ rei Þ
Ψ e ððdeij Þ2 Þ ¼
ð37Þ
In other words, Ψ should be an isotropic scalar function of the vectors rei . Now by virtue of well known representation theorems for isotropic functions of vectors (see [6]), the requirement in Eq. (37) implies that Ψe is a function of the integrity basis for its vector arguments. In the current case, this implies that the strain energy can depend only upon all possible inner products between the vectors, that is ðΨ Þn Ψ ððrei Þn Þ ¼ Ψ ðQ rei Þ ) Ψ ¼ Ψ ðrei rej Þ ði; j ¼ 1; …; nÞ e
e
e
e
NNE X NNE X i
39
Ψ eij ððdeij Þ2 Þ
ð44Þ
j4i
That is, for the peridynamics approaches, the strain energy is written as just the sum of the strain energies along each edge with no cross dependence. On the other hand for FEM formulations, the strain energy depends upon the distance between all the edges of the element and is not just the sum. This is a crucial feature that is present in FEM formulations that provides them with a richer structure than truss or peridynamics formulations.
e
ð38Þ
5.3. Infinitesimal deformations
Next using the fact that e
rei
rej
¼
J rei J 2 þ J rej J 2 ðdij Þ2 2
;
e dij
¼ Jρ
e ij J
¼
J rei rej J
ð39Þ
Here ρeij denotes the vector connecting element node i to element node j (see Fig. 3). We can further write Eq. (38) in the equivalent form
Ψ e ¼ Ψ e ð J rei J 2 ; ðdeij Þ2 Þ ði; j ¼ 1; …; nÞ
ð40Þ
Thus invariance under pure rotation about the origin implies that the strain energy function can be dependent only upon the distance of the nodes from the origin and the distance of the nodes from each other. Next, we consider the translation vector c and note that
Ψ e ¼ Ψ e ð J rei J ; ðdeij Þ2 Þ ¼ Ψ e ð J rei þ ce J 2 ; ðdeij Þ2 Þ ði; j ¼ 1; …; nÞ
ð41Þ
we now choose ce ¼ re1 and arrive at
Ψ ¼ Ψ ðJρ e
e
e 2 e 2 i1 J ; ðdij Þ Þ ¼
Ψ
e
e ððdij Þ2 Þ
ði; j ¼ 1; 2; …; nÞ
The results established above do not assume that either the displacements or the strain energy are in any way restricted. For the case of infinitesimal deformations where the strain energy can be at most quadratic, we can obtain the result in a direct way as follows: For the case of infinitesimal deformations, we see that if we write
ρeij ¼ Reij þ ueij
where R eij ¼ R ei R ej is the relative position of node i with respect to node j in the undeformed configuration and ueij ¼ uei uej is the relative displacement of node i with respect to node j. Now a routine calculation of ρeij ρeij , keeping terms that are at most linear in ueij (and assuming that the undeformed configuration is one of minimum energy) reveals that
ψ e ¼ ψ e ðξeij Þ;
ð42Þ
e
5.2. Nodal forces derived from the strain energy Now that we have shown that the requirement that the strain energy of a body that is discretized into finite elements, should be invariant to translations and rotations of the elements implies that the strain energy has to have the simplified form (35), we go on to show that this implies that the nodal forces are always directed along the vectors connecting the nodes of the element. The magnitude of these forces can depend upon the distances from all the nodes that share the element but its direction is along the edges that connect the nodes. To see this, we note that for elastic bodies, the nodal forces for each element are given by an expression similar to Castigliano's Theorem I [3], and results in the following form: e e n X ∂Ψ ∂Ψ e ¼2 e 2 ρij e ∂ri j ¼ 1 ð∂dij Þ
ð46Þ
e
fi ¼
e
ξeij ¼ ueij R eij no sum on i and j
where ξij is the axial extension1 along the edge that connects nodes i and j in element e. This in turn implies that the nodal force is given by
and thus our proof is completed.□
fi ¼
ð45Þ
e e n X ∂Ψ ∂Ψ e ¼ e R ij e ∂ri j ¼ 1 ξij
ð47Þ
6. Conclusions In this paper we established a very simple but revealing result concerning the forces at the nodes of finite element representations of hyperelasticity. In particular, we established that forces at a node are only functions of the extensional strains along the edges and thus the finite element formulation is “truss-like” but with the forces along the truss members being dependent upon the strains of the neighboring truss members. This result shows that these FEM formulations can be treated in a manner similar (but not identical) to truss systems.
ð43Þ
where we have also used the chain rule for differentiation to simplify the expression. Upon taking the magnitude of the above expression and noting e e e that J ρeij J ¼ dij we see that forces of magnitude ð∂ψ e =∂ðdij Þ2 Þdij are acting along each edge in the element connecting the node i with the other nodes j. The similarity of this expression with that for a truss system and for peridynamics [7] is noteworthy. However there is one key difference, for a simple truss system and for peridynamics, it is
Acknowledgments The authors gratefully acknowledge the support of this work by their respective named endowed professorships. The authors are also happy to acknowledge the constructive comments by Ms. Parisa Khodabkhshi on the contents of the paper. 1 It can easily be shown that the axial extension ξij defined in Eq. (46) remains invariant under infinitesimal rotations as befits a small deformation theory.
e
40
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
References [1] J.N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill, New York, 2006. [2] J.N. Reddy, Introduction to Continuum Mechanics, 2nd ed., Cambridge University Press, New York, 2013. [3] J.N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, 2nd ed., John Wiley & Sons, New York, 2002.
[4] J.N. Reddy, Applied Functional Analysis and Variational Methods in Engineering, McGraw-Hill, New York, 1986. [5] J.N. Reddy, An Introduction to Nonlinear Finite Element Analysis, Oxford University Press, Cambridge, UK, 2015. [6] A.J.M. Spencer, R.S. Rivlin, Isotropic integrity bases for vectors and second-order tensors, Arch. Ration. Mech. Anal. 9 (1) (1962) 45–63. [7] S.A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (2000) 175–209.
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
On the force–displacement characteristics of finite elements for elasticity and related problems J.N. Reddy n, Arun Srinivasa Department of Mechanical Engineering, Texas A & M University, College Station, TX 77843-3123, USA
art ic l e i nf o
a b s t r a c t
Article history: Received 23 October 2014 Received in revised form 25 February 2015 Accepted 28 April 2015
The displacement based finite element models for hyperelastic materials are considered, and it is shown that —based on the invariance of the discretized strain energy to translations and rigid body rotations of the elements—the magnitude of the nodal forces can be written purely in terms of the strains along the edges composing the elements, while their directions are along the unit vectors corresponding to the edges. In other words, the nodal forces are “truss like” with one major difference: the magnitude of the forces depends upon the extensional strains of neighboring edges. In order to motivate the ideas, we first explicitly demonstrate this result for triangular finite element formulation of plane elasticity problems and then show that the result is valid for any finite element formulation of hyperelasticity. & 2015 Elsevier B.V. All rights reserved.
Keywords: Plane elasticity Triangular finite elements Extensional strains along edges Mathematical proof
1. Introduction The aim of this paper is to show a relatively simple result that has hitherto not been documented in the FEM literature, despite a long history of finite element analysis of structural problems (which were originally motivated by matrix structural analysis of trusses). The idea is stated simply thus: for hyperelasticity problems, the magnitude of the nodal forces can be written purely in terms of the strains along the edges composing the elements, while their directions are along the unit vectors corresponding to the edges, i.e., no matter what the element type, we can show that the problem can always be reduced to a truss-like problem (with the “links” of the truss being represented by straight lines connecting the nodes of each element), but with one difference—namely that the magnitude of each of the forces at a node (which are along each of the links connecting the node to other nodes, like in a truss) depends not only upon the strain in the link but also on the strains along the links connecting the other nodes. This provides a surprisingly simple picture of the finite element method but provides additional insight and new possibilities for generalization and solution of problems. In order to prove this result, we proceed in two steps: in the first step, we start with the traditional FEM formulation of the isotropic 2-D elasticity problem using the linear triangular element and show explicitly that the magnitude of the forces at the nodes of a constant strain triangular (CST) element depends only upon the extensional strains of the edges of the elements that share the node, and that the directions of the nodal forces are along the edges that have the node in common. In other words, the problem is like that of a generalization of n
Corresponding author. Tel.: þ 979 862 2417. E-mail address: [email protected] (J.N. Reddy).
http://dx.doi.org/10.1016/j.finel.2015.04.011 0168-874X/& 2015 Elsevier B.V. All rights reserved.
a “truss”, i.e. the forces are along the truss elements but their magnitudes dependent upon the strains in neighboring truss elements. We then show that further that the result is not specific to a constant strain triangle, but is true of all finite element formulations of elasticity problems (both linear and non-linear). The result follows from the fact that the discretized strain energy of the finite element representation has to (1) depend only upon nodal positions (once the shape functions are chosen) and (2) must be invariant to translations and rigid body rotations of the elements, together with well-established representation theorems using the integrity bases for isotropic functions. 2. Governing equations The governing equations for plane elasticity (plane stress or plane strain) problems are summarized here (see Reddy [1,2]). The presentation here is limited to linearized elasticity, where no distinction is made between the Cauchy stress tensor and the second Piola–Kirchhoff stress tensor and between the current coordinates x and the reference coordinates X. The domain of the elastic solid is denoted by Ω, which is an open bounded subset of R2 . The boundary of Ω is denoted by Γ ¼ ∂Ω ¼ Ω Ω, where Ω represents the closure of Ω. A typical point belonging to Ω is denoted as x. We employ the customary designations for the Sobolev spaces H s ðΩÞ and H s ðΓ Þ, where s Z0 [4,5]: Z Z H s ðΩÞ ¼ u : j Ds uj 2 dΩ o1 ; H s ðΓ Þ ¼ u : j Ds uj 2 dΓ o 1 Ω
Γ
ð1Þ
36
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
The corresponding norms are denoted as J J Ω;s and J J Γ ;s . Likewise the inner products associated with these spaces are denoted as ð; ÞΩ;s and ð; ÞΓ ;s respectively. The product spaces Hs ðΩÞ ¼ ½H s ðΩÞ2 are constructed in the usual way. The governing equations of a plane elastic continuum consist of equations of equilibrium, strain–displacement relations, and stress– strain relations. These equations are expressed in matrix form. Equations of equilibrium: DT σ þ f ¼ 0
ð2Þ
where (σ xx ; σ yy ; σ xy ) are the rectangular Cartesian components of the stress tensor σ , fx and fy are the components of the body force vector f (measured per unit volume) along the x- and y-directions, respectively, and 8 9 " # ( ) σ = > < xx > fx ∂=∂x 0 ∂=∂y T σ yy ; f¼ ; σ¼ D ¼ ð3Þ f 0 ∂=∂y ∂=∂x > > y :σ ; xy
Strain–displacement relations:
εxx ¼
∂ux ; ∂x
or
ε ¼ Du;
εyy ¼
∂uy ; ∂y
2εxy ¼
8 9 ε = > < xx > ε ¼ εyy ; > : 2ε > ; xy
∂ux ∂uy þ ∂y ∂x
( u¼
ux
ð4Þ
)
uy
ð5Þ
where (εxx ; εyy ; εxy ) are the rectangular Cartesian components of the strain tensor ε and ux and uy are the rectangular Cartesian components of the displacement vector u. Stress–strain (or constitutive) relations: 8 9 2 9 38 σ = ε = c11 c12 0 > > < xx > < xx > 7 σ yy ¼ 6 ð6Þ 4 c12 c22 0 5 εyy > > > :σ > ; 0 0 c33 : 2εxy ; xy or
σ ¼ Cε;
2
c11 6 C ¼ 4 c12 0
c12
0
3
0 7 5 c33
c22 0
ð7Þ
where cij (cji ¼ cij Þ are the elasticity (material) constants for an orthotropic medium with the material principal directions (x1 ; x2 ; x3 ) coinciding with the coordinate axes ðx; y; zÞ used to describe the problem. The cij can be expressed in terms of the engineering constants ðE1 ; E2 ; ν12 ; G12 Þ for an orthotropic material, say for plane stress problems, by c11 ¼
E1
;
c22 ¼
ð1 ν12 ν21 Þ c12 ¼ ν12 c22 ¼ ν21 c11 ;
E2
ð1 ν12 ν21 Þ c33 ¼ G12
ð8Þ
Eqs. (2)–(8) can be combined into a single equation in terms of the displacement vector u as DT CDu þ f ¼ 0
ð9Þ
Eq. (9) contains, for the two-dimensional case, two equations of equilibrium in terms of the two displacement components ux and uy.
where Bðw; uÞ is a bilinear form, F ðwÞ is a linear form, Z Bðw; uÞ ¼ h ðDwÞT CðDuÞ dx e I ZΩ wT f dx þ h wT t ds F ðwÞ ¼ h Ωe
Γe
Here h denotes the thickness of the elastic solid, ðÞT denotes the transpose of the enclosed quantity, V and W are vector spaces (e.g., appropriate subsets of the product (Sobolev) space H1 ðΩÞ ¼ H 1 ðΩÞ H 1 ðΩÞ). The quantity w δu represents the weight or test function. Unlike classical solutions that are defined unambiguously point-wise, weak solutions exist with respect to test functions and are therefore understood in the context of distributions. As a result, the weak solution (and its derivatives) is typically defined unambiguously in Ω up to a set of measure zero. 4. Finite element model The finite element model associated with Eq. (11) is obtained by restricting the solution space to a finite-dimensional subspace V h of the infinite-dimensional function space V, and the weight function to a finite-dimensional sub-space W h W. As a result, in the discrete problem we seek uh A V h such that Bðwh ; uh Þ ¼ F ðwh Þ
8 wh A W h
The strains can be computed as
ε ¼ Du ¼ DΨΔ BΔ;
σ ¼ CBΔ
To obtain the finite element model, we substitute Eq. (13) into the variational statement in Eq. (10) Z Z I e e e T e T 0 ¼ he ðδΔ ÞT BT CBΔ dx he ðδΔ ÞT Ψ f dx he ðδΔ ÞT Ψ t ds Ωe
or
The principle of virtual displacements [3] is equivalent to the following variational problem: find u A V such that
where
ð10Þ
ð14Þ
where D is defined earlier in Eq. (3) and B is a 3 2n matrix 2 ∂ψ 3 ψ2 ψn 1 0 ∂∂x 0 ⋯ ∂∂x 0 6 ∂x ∂ψ ψ2 ψn 7 1 6 7 0 ∂∂y ⋯ 0 ∂∂y B ¼ DΨ ¼ 6 0 ð15Þ 7 ∂y 4 ∂ψ ∂ψ ∂ψ ∂ψ 5 ∂ψ n ∂ψ n 1 1 2 2 ⋯ ∂y ∂x ∂y ∂x ∂y ∂x
e T
¼ ðδΔ Þ K Δ F Q
8wAW
ð12Þ
The quantity h in the definition of the subspaces V h and W h represents the average size of all the elements in a given finite element discretization. We assume that the domain Ω R2 is e discretized into a set of NE non-overlapping sub-domains Ω , h e NE called finite elements, such that Ω Ω ¼ ⋃e ¼ 1 Ω . The geometry of each element is characterized using the standard isoparametric ^ e to the physical bijective mapping from the master element Ω e element Ω . e Let u be approximated over Ω by the finite element interpolations (the element label e on the variables is omitted in the interest of brevity) ( ) ( ) ux w 1 ¼ δ ux u¼ ¼ ΨΔ; w ¼ δu ¼ ¼ ΨδΔ uy w 2 ¼ δ uy " # ψ1 0 ψ2 0 … ψn 0 Ψ¼ 0 ψ 0 ψ2 … 0 ψn 1 n oT Δ ¼ u1x u1y u2x u2y … unx uny ð13Þ
3. The principle of virtual displacements and the variational problem
Bðw; uÞ ¼ F ðwÞ
ð11Þ
e
e
e
Ωe
e
Γe
ð16Þ
Ke Δ ¼ Fe þ Q e e
Ke ¼ he
ð17Þ
Z
Z Ωe
BT CB dx;
Fe ¼ he
Ω
Ψ T f dx; e
Q e ¼ he
I Γe
ΨT t ds
ð18Þ
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
The element stiffness matrix Ke is of the order 2n 2n and the element load vector Fe and the vector of internal forces Q e is of order 2n 1, where n is the number of nodes in a Lagrange finite element. For the linear right-angled triangular (master) element (n ¼3) with horizontal side of a and vertical side of b, matrix B is given by 2 1 3 1 a 0 0 0 0 a 6 7 0 1b 0 0 0 1b 7 B ¼ DΨ ¼ 6 ð19Þ 4 5 1 1 1 1 b a 0 a b 0 Fig. 1 shows the nodal displacements and forces of the linear master triangular element. The stiffness matrix takes the form (area of the triangle, A ¼ ab=2 and h is the thickness of the body) 2c 3 þ c33 Þ 11 33 12 þ c332 ðc12 ab ca112 cab c332 cab a2 b b 6 7 6 ðc12 þ c33 Þ c33 þ c22 c12 c33 c33 c22 7 2 2 6 ab ab ab a2 a2 b b 7 6 7 c11 c11 c12 7 12 cab 0 0 abh6 6 a2 e ab 7 a2 K ¼ ð20Þ 6 7 c33 c33 33 2 6 cab ca332 0 0 7 ab a2 6 7 c33 c33 c33 6 c33 7 ab 0 0 7 2 2 6 ab b b 4 5 c12 c22 c12 c22 ab 2 0 0 2 ab b
The nodal forces are given by 2c þ c33 Þ 11 þ c332 ðc12 ab a2 b 6 6 ðc12 þ c33 Þ c33 þ c22 2 6 ab a2 b 6 c11 12 cab abh6 6 a2 ¼ 6 33 > > F 2y > 2 6 cab ca332 > > > 6 > > > > c 6 33 33 F > > 3x > cab > 6 b2 > > > 4 ; : F 3y > c12 c22 ab 2 b
ca112
33 cab
12 cab
ca332
c11 a2
c332 b 33 cab
0
0
0
c33 a2 c33 ab
c33 ab c33 2 b
c12 ab
0
0
0
38 9 12 cab > u1x > > 7> > > > 1> > > > c222 7 > > uy > > b 7 > 7> 2 > c12 7> = x 7 ab 7 2 > uy > 0 7 > > > 7> > > > u3 > > 0 7 > > 7> x > > 5> > > 3> ; : c22 u y 2 b
ð21Þ where Fix and Fiy are the horizontal (x-direction) and vertical (ydirection) forces, respectively, at node i, i ¼ 1; 2; 3 (see Fig. 2). For example, at node 1 we have abh c11 c33 1 ðc12 þ c33 Þ 1 F 1x ¼ uy þ 2 ux þ 2 2 ab a b c11 c33 c33 c12 2 u2x u2y 2 u3x u3y ab ab a b abh c11 1 c33 ¼ u u2x þ 2 u1x u3x 2 a2 x b c
i c12 1 33 uy u3y þ u1y u2y þ ab ab abh c11 c c ¼ εxx þ 12 εyy þ 33 γ xy ð22Þ 2 a a b abh c12 þ c33 1 c33 c22 1 þ 2 uy ux þ F 1y ¼ 2 2 ab a b c12 2 c33 2 c33 3 c22 3 ux 2 uy ux 2 uy ab ab a b c
c22
abh c12 1 c33 1 33 2 ux ux þ ux u3x þ 2 u1y u3y þ 2 u1y u2y ¼ 2 ab ab a b abh c12 c c εxx þ 22 εyy þ 33 γ xy 2 b b a
ð23Þ
where
εyy ¼ εð2Þ ¼
u3y u1y
;
εxx ¼ εð3Þ ¼
b u3x u1x u2y u1y þ γ xy ¼ γ 23 ¼ b a
u x3 u 3y
Side 1 Undeformed element Side 2
3
Deformed element
3
θ
b
2
ux2
1 1
u1x
u1y a
2
u 2y
x
Side 3
Fig. 1. The linear triangular element with the 6 displacement degrees of freedom.
b
9 8 F 1x > > > > > > >F > > > > 1y > > > > > > > < F 2x =
¼
y
37
u2x u1x ; a
From the strain transformation equations, we have
εn ¼ εxx sin 2 θ γ xy cos θ sin θ þ εyy cos 2 θ
ð24Þ
Fig. 2. The linear triangular element with the 6 displacement degrees of freedom.
with a sin θ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 2 a2 þ b
b cos θ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 þ b
Thus, defining the axial strain along the edges 2–3 as εð1Þ and using the above strain transformation rules, we can obtain the shear strains in terms of the three axial strains εðiÞ (i¼1, 2, 3) as a b
b a
γ xy ¼ εð3Þ þ εð2Þ
2
a2 þ b ð1Þ ε ab
ð25Þ
Therefore, we can express F1x and F1y solely in terms of the extensional (normal) strains in sides 1, 2, and 3: " !# 2 abh c11 ð3Þ c12 ð2Þ c33 a ð3Þ b ð2Þ a2 þ b ð1Þ ε þ ε þ ε þ ε ε F 1x ¼ 2 a a a b b ab " # 2 abh a2 þ b 1 c11 c33 ð3Þ ð1Þ ð2Þ c33 ε ðc12 þ c33 Þε a 2 þ 2 ε ¼ ð26Þ 2 2 a a ab b " !# 2 abh c12 ð3Þ c22 ð2Þ c33 a ð3Þ b ð2Þ a2 þb ð1Þ ε þ ε þ ε þ ε ε 2 a b b a b ab " # 2 2 abh a þ b c22 c33 ð2Þ 1 ð1Þ ð3Þ c c ε b þ ε ð þc Þ ε ¼ 33 12 33 2 2 b a2 b a2 b
F 1y ¼
Similarly, abhhc12 ð2Þ c11 ð3Þ i ε þ ε ; F 2x ¼ 2 a a
F 3x ¼
F 2y ¼
ð27Þ
i c33 hh 2 ð3Þ 2 2 a ε þ b εð2Þ ða2 þ b Þεð1Þ 2a ð28Þ
i c33 hh 2 ð3Þ 2 2 a ε þ b εð2Þ ða2 þb Þεð1Þ ; 2b
F 3y ¼
abhhc22 ð2Þ c12 ð3Þ i ε þ ε 2 b b ð29Þ
One can verify that the nodal forces are in equilibrium (i.e., F 1x þ F 2x þ F 3x ¼ 0 and F 1y þ F 2y þF 3y ¼ 0). If we define the unit vectorspalong ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the three edges as 2 e1 2 ¼ i; e1 3 ¼ j, e2 3 ¼ ð ai þ bjÞ= a2 þ b Þ (where i and j are the unit vectors along the coordinate axes, it is easy to see that the nodal forces at node 1 are trivially along the edges 1–2 and 1–3. Next, if we consider the nodal forces at node 2, we see that (using
38
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
the geometry of the triangular element) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2
a a2 þ b F 2y e2 3 F2 ¼ F 2x iþ F 2y j ¼ F 2x þ F 2y e1 2 þ b b
ð30Þ
and
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 b a2 þ b F3 ¼ F 3x iþ F 3y j ¼ F 3y þ F 3x e1 3 F 3x e2 3 a a
ð31Þ
In other words, the nodal forces at each node can be written in terms of the unit vectors along the edges. For the case of b ¼a, we have the result ah 2c33 εð1Þ ðc12 þc33 Þεð2Þ ðc11 þ c33 Þεð3Þ 2 ah 2c33 εð1Þ ðc22 þ c33 Þεð2Þ ðc12 þ c33 Þεð3Þ F 1y ¼ 2 ah c33 ah ð1Þ c12 εð2Þ þc11 εð3Þ ; F 2y ¼ 2ε εð2Þ εð3Þ F 2x ¼ 2 2 ah ah c22 εð2Þ þ c12 εð3Þ F 3x ¼ c33 2εð1Þ εð2Þ εð3Þ ; F 3y ¼ 2 2 F 1x ¼
Ψ e ðrei Þ ¼ Ψ ððrei Þn Þ;
ð32Þ
5. Strain energy function and nodal forces associated with a finite element discretization of hyperelasticity 5.1. Form of the strain energy function Having explicitly demonstrated that the nodal forces for right angled triangular elements for 2-D isotropic plane elasticity can be recast in the form of a “truss-like” problem, we now turn to the general problem of FEM formulations of hyperelasticity and prove that the result is not a special quirk of triangular elements but it follows from fundamental invariance properties of the strain energy function. Here, we need to generalize the notion of an “edge” as being “the straight line joining any two nodes of a finite element”. Thus, a four-node quadrilateral has six edges, four along the boundary and two along the diagonals. Similarly an eight-node cube will have 28 edges. Thus, we suppose that the elastic medium Ω has been h discretized into a finite element mesh Ω consisting of NN nodes and NE elements. It is evident that, for a Green elastic material, given the interpolation functions and the location of the nodes, the strain energy Ψ of the discretized domain can be determined entirely by the position vectors ri of the NN nodes, since these are the generalized coordinates, that (along with the fixed interpolation functions) uniquely determine the configuration of the deformed body. Thus, for a finite element model of a Green elastic body, we must have NE X
Ψ e ðrei Þ
ðrei Þn ¼ Q e rei þ ce
ð34Þ
ðrei Þn e
A routine calculation further reveals that the direction of the nodal force vectors is along the edges of the triangle. In other words, the FEM formulation of the plane problem in elasticity using right angled triangular elements resembles a truss formed by the nodes and edges of the elements. The nodal forces are along the truss elements (edges of the FEM formulation) while their magnitudes depend upon the extensional strains of all the neighboring truss elements. This is not a surprising result for this specific case, especially considering the special geometry chosen. However, we now show that for any finite element discretization of hyperelasticity, the above result holds.
Ψ¼
elements, i runs from 1 to 3. If we choose some other higher order element that i runs from 1 to the number of nodes per element. Thus the notation Ψe stands for the strain energy of the eth element, while the notation rei stands for the current position of the ith node of the eth element (see Fig. 3). Irrespective of how the discretized strain energy is computed, it has to meet a fundamental criterion: translations and rigid body rotations (in the case of infinitesimal deformations, it is the infinitesimal rotation) superimposed on a given configuration of any element (i.e. on any given state of strain) should not change the strain energy function . Thus the elemental strain energy Ψe must satisfy (for each element)
ð33Þ
is the position vector with superposed rigid-body where motion, Q is a proper orthogonal tensor, and ce is a constant vector. The above requirement is the counterpart of the notion of invariance under superposed rigid body motions [2] for all classical continua. We now state and prove one of the important results of this study. Theorem 1. A necessary and sufficient condition for the strain energy to meet the condition in Eq. (34) is that
Ψ e ¼ Ψ e ððdeij Þ2 Þ e dij
ð35Þ
J rei rej J
¼ is the distance between the ith and jth nodes of where the eth element (see Fig. 3). Proof. The fact that form in Eq. (35) is sufficient is a consequence of the fact that the distance between any two nodes of an element is unaffected by translations and rigid body rotations of the element. In other words, it is obvious from Eq. (34), the properties of orthogonal tensors (i.e., Q T Q ¼ I, the identity tensor and Q ¼ 7 1), and the fact that ce is a constant that 2 n n n n n n n e ¼ J rei rej J 2 ¼ rei rej rei rej dij
¼ Q rei rej Q rei rej
2
e ¼ rei rej Q T Q rei rej ¼ dij ð36Þ
y NN
26
•3
27 Ωe
2
•
1
•15
24 13 12
3•
2
e
•
•
15
1
1
2
r15 = r1e
3
4
r 12e
7
10
23 22 11
Global node numbers
Element node numbers
e¼1
where superscript “e” represents the element number, a lower case subscript denotes the node number within an element, that is, ðÞei stands for the enclosed quantity at the ith node of the eth e element Ω . Thus if we choose to use mesh of triangular
x Fig. 3. A finite element mesh of linear triangular elements showing the position of the 15th global node and its node number relative to the eth element Ωe. Thus the node, which has a global number of 15, is the 1st node of the eth element.
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
To prove necessity, we note that since Eq. (34) is valid for all c we first choose c ¼ 0 and obtain
assumed that
Ψ e ððrei Þn Þ ¼ Ψ e ðQ rei Þ
Ψ e ððdeij Þ2 Þ ¼
ð37Þ
In other words, Ψ should be an isotropic scalar function of the vectors rei . Now by virtue of well known representation theorems for isotropic functions of vectors (see [6]), the requirement in Eq. (37) implies that Ψe is a function of the integrity basis for its vector arguments. In the current case, this implies that the strain energy can depend only upon all possible inner products between the vectors, that is ðΨ Þn Ψ ððrei Þn Þ ¼ Ψ ðQ rei Þ ) Ψ ¼ Ψ ðrei rej Þ ði; j ¼ 1; …; nÞ e
e
e
e
NNE X NNE X i
39
Ψ eij ððdeij Þ2 Þ
ð44Þ
j4i
That is, for the peridynamics approaches, the strain energy is written as just the sum of the strain energies along each edge with no cross dependence. On the other hand for FEM formulations, the strain energy depends upon the distance between all the edges of the element and is not just the sum. This is a crucial feature that is present in FEM formulations that provides them with a richer structure than truss or peridynamics formulations.
e
ð38Þ
5.3. Infinitesimal deformations
Next using the fact that e
rei
rej
¼
J rei J 2 þ J rej J 2 ðdij Þ2 2
;
e dij
¼ Jρ
e ij J
¼
J rei rej J
ð39Þ
Here ρeij denotes the vector connecting element node i to element node j (see Fig. 3). We can further write Eq. (38) in the equivalent form
Ψ e ¼ Ψ e ð J rei J 2 ; ðdeij Þ2 Þ ði; j ¼ 1; …; nÞ
ð40Þ
Thus invariance under pure rotation about the origin implies that the strain energy function can be dependent only upon the distance of the nodes from the origin and the distance of the nodes from each other. Next, we consider the translation vector c and note that
Ψ e ¼ Ψ e ð J rei J ; ðdeij Þ2 Þ ¼ Ψ e ð J rei þ ce J 2 ; ðdeij Þ2 Þ ði; j ¼ 1; …; nÞ
ð41Þ
we now choose ce ¼ re1 and arrive at
Ψ ¼ Ψ ðJρ e
e
e 2 e 2 i1 J ; ðdij Þ Þ ¼
Ψ
e
e ððdij Þ2 Þ
ði; j ¼ 1; 2; …; nÞ
The results established above do not assume that either the displacements or the strain energy are in any way restricted. For the case of infinitesimal deformations where the strain energy can be at most quadratic, we can obtain the result in a direct way as follows: For the case of infinitesimal deformations, we see that if we write
ρeij ¼ Reij þ ueij
where R eij ¼ R ei R ej is the relative position of node i with respect to node j in the undeformed configuration and ueij ¼ uei uej is the relative displacement of node i with respect to node j. Now a routine calculation of ρeij ρeij , keeping terms that are at most linear in ueij (and assuming that the undeformed configuration is one of minimum energy) reveals that
ψ e ¼ ψ e ðξeij Þ;
ð42Þ
e
5.2. Nodal forces derived from the strain energy Now that we have shown that the requirement that the strain energy of a body that is discretized into finite elements, should be invariant to translations and rotations of the elements implies that the strain energy has to have the simplified form (35), we go on to show that this implies that the nodal forces are always directed along the vectors connecting the nodes of the element. The magnitude of these forces can depend upon the distances from all the nodes that share the element but its direction is along the edges that connect the nodes. To see this, we note that for elastic bodies, the nodal forces for each element are given by an expression similar to Castigliano's Theorem I [3], and results in the following form: e e n X ∂Ψ ∂Ψ e ¼2 e 2 ρij e ∂ri j ¼ 1 ð∂dij Þ
ð46Þ
e
fi ¼
e
ξeij ¼ ueij R eij no sum on i and j
where ξij is the axial extension1 along the edge that connects nodes i and j in element e. This in turn implies that the nodal force is given by
and thus our proof is completed.□
fi ¼
ð45Þ
e e n X ∂Ψ ∂Ψ e ¼ e R ij e ∂ri j ¼ 1 ξij
ð47Þ
6. Conclusions In this paper we established a very simple but revealing result concerning the forces at the nodes of finite element representations of hyperelasticity. In particular, we established that forces at a node are only functions of the extensional strains along the edges and thus the finite element formulation is “truss-like” but with the forces along the truss members being dependent upon the strains of the neighboring truss members. This result shows that these FEM formulations can be treated in a manner similar (but not identical) to truss systems.
ð43Þ
where we have also used the chain rule for differentiation to simplify the expression. Upon taking the magnitude of the above expression and noting e e e that J ρeij J ¼ dij we see that forces of magnitude ð∂ψ e =∂ðdij Þ2 Þdij are acting along each edge in the element connecting the node i with the other nodes j. The similarity of this expression with that for a truss system and for peridynamics [7] is noteworthy. However there is one key difference, for a simple truss system and for peridynamics, it is
Acknowledgments The authors gratefully acknowledge the support of this work by their respective named endowed professorships. The authors are also happy to acknowledge the constructive comments by Ms. Parisa Khodabkhshi on the contents of the paper. 1 It can easily be shown that the axial extension ξij defined in Eq. (46) remains invariant under infinitesimal rotations as befits a small deformation theory.
e
40
J.N. Reddy, A. Srinivasa / Finite Elements in Analysis and Design 104 (2015) 35–40
References [1] J.N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill, New York, 2006. [2] J.N. Reddy, Introduction to Continuum Mechanics, 2nd ed., Cambridge University Press, New York, 2013. [3] J.N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, 2nd ed., John Wiley & Sons, New York, 2002.
[4] J.N. Reddy, Applied Functional Analysis and Variational Methods in Engineering, McGraw-Hill, New York, 1986. [5] J.N. Reddy, An Introduction to Nonlinear Finite Element Analysis, Oxford University Press, Cambridge, UK, 2015. [6] A.J.M. Spencer, R.S. Rivlin, Isotropic integrity bases for vectors and second-order tensors, Arch. Ration. Mech. Anal. 9 (1) (1962) 45–63. [7] S.A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (2000) 175–209.