- Email: [email protected]
High order anisotropic finite elements for three-dimensional isotropic hyperelastic continua
238
High order anisotropic finite elements for three-dimensional isotropic hyperelastic continua A. Duster ^'*, S. Hartmann^ ^Technische Universitdt Miinchen, Fakultdt fur Bauingenieur- und Vermessungswesen, Lehrstuhl fiir Bauinformatik, Arcisstn 21, 80290 Miinchen, Germany ^ Universitdt Kassel, Institut fiir Mechanik, Mdnchebergstr 7, 34109 Kassel, Germany
Abstract In this presentation we apply the p-version finite element method to three-dimensional hyperelastic bodies. On the basis of a pure displacement formulation we take hyperelastic constitutive equations for near incompressibility into account. In order to discretize efficiently three-dimensional continua, an anisotropic element formulation of high order is applied, yielding locking-free approximations for thick as well as thin-walled solids at finite strains. A numerical example demonstrates the basic properties of our nonUnear p-version approach for hyperelastic problems. Keywords: p-Version finite element method; Anisotropic elements; Hyperelasticity; Finite strains
1. Introduction Usually, the finite element method using low order shape functions is applied to finite strain problems. This leads in the case of nearly incompressible material behaviour to volumetric locking phenomena. In order to overcome these problems, a series of mixed finite elements have been proposed (see, for example, Belytschko et al. [1], Reese and Wriggers [7] and the literature cited therein). On the other hand high order finite elements do not show, beyond a certain order of the shape functions, locking effects, which has been investigated by several authors, see for example Suri [8] and Szabo and Babuska [9]. Here, we combine anisotropic p-version finite elements and the finite strain theory to compute three-dimensional continua with nearly incompressible, isotropic hyperelastic material.
deformation gradient. The unimodular left Cauchy-Green tensor B = FF is used to formulate a strain-energy function of the form ^(/,lB,IlB) = f/(/) + u;(%,IlB),
where % = trB and Ilf = trB denote the first and second invariant of B. Due to Hartmann and Neff [6] we assume a strain-energy function V^(/,IB,IIB) with
^ ( / ) = ^ (/' + / ' ' - 2 ) ,
which is polyconvex and coercive, i.e. there exists a solution of the underlying boundary-value problem. The Cauchy stress decomposes in a natural manner into a spherical and a deviatorical part resulting from the specified decomposition of the deformation gradient. T = UXJ)l +
* Tel.: +49 (89) 2892-5060; Fax: +49 (89) 2892-5051; E-mail: [email protected] © 2003 Elsevier Science Ltd. All rights reserved. Computational Fluid and Solid Mechanics 2003 K.J. Bathe (Editor)
(2)
U;(IB,IIB) = Of(II - 27) + ciode - 3) + Coi (II^/' - 3x/3), (3)
2. Constitutive model For the case of the aforementioned investigations, we restrict ourselves to isotropic hyperelasticity. We decompose the deformation gradient F = FF into a volume changing part F = /^/^I and a volume preserving part F = /~^/^F, detF = 1. / = detF symbolizes the determinant of the
(1)
^fe)^ J\6R
(4)
)
where T defines the Cauchy stress tensor and I symbolizes the second order identity tensor. The superscript D denotes the deviator operator, A^ = A — l/3(tr A)I, while trA = a,' defines the trace operator. Fig. 1 shows the basic uniaxial behaviour of rubber materials, i.e. it represents an S-shaped curve in a force-stretch diagram. For more details we refer to [6].
A. Duster, S. Hartmann /Second MIT Conference on Computational Fluid and Solid Mechanics
239
10 c3
<°
2 ^
-10 -15 0.5
1
1.5 2 stretch L/LQ
2.5
Fig. 1. Uniaxal stress-stretch behaviour (material parameters, see Section 4).
3. p-Version FEM at finite deformations The implementation of the /7-version is based on a hexahedral element, allowing for an anisotropic Ansatz of the displacement field (see Fig. 2). The polynomial degree of each separate component of the displacement field can be chosen individually and may also be varied in the three local directions of the element. This anisotropic Ansatz allows the efficient computation of three-dimensional plate and shelllike structures. For a detailed presentation of the element formulation and its application to different linear and nonlinear problems of solid mechanics we refer to [2-5].
Ql = [(-1,1) X (-1,1) X (-1,1)]
The overall solution process is realized by applying an incremental/iterative procedure. Therefore, the nonlinear system of equations arising from the spatial discretization based on the p-version is linearized by means of a Newton-Raphson scheme. The load is applied incrementally, including an adaptation of the load increment which is controlled by the number of iterations of the NewtonRaphson scheme. Since the presented numerical example is highly nonlinear, the adaptation of the load steps was found to be crucial for computing a numerical solution.
S^r'^'^iQllS^^''^'^(QllSP^P^HQl)
Fig. 2. Standard hexahedral Q^^: definition of nodes, edges, faces and polynomial degree.
240
A. Duster, S. Hartmann / Second MIT Conference on Computational Fluid and Solid Mechanics
t t t t
!-V " " "^
point A t
t = 1 [mm]
/W;c=0
/
n/
\
/
\
.iL^.^i
5 [mm],
G
10 [mm]
I
5 [mm] I
Fig. 3. Inhomogeneous compression test (geometry and element mesh).
4. Numerical examples
in such a way that the local coordinate f of all elements coincides with the z-axis. The convergence of the displacement component M^ of point A is plotted in Fig. 4 for both an isotropic and an anisotropic increase of the polynomial degree pk,k = ^,r],^ based on the trunk space 5tf'^"''^(^^t) (for the definition of the trunk space see [2,4]). In the case of the isotropic discretization the polynomial degree of the elements is increased in all three local directions (|,y;,f) simultaneously, whereas in the case of the anisotropic discretization the polynomial degree is only increased in in-plane direction (§,/;) and is kept fixed at p^ = 1 in f-direction in order to account for the plain strain assumptions. From this it is obvious that the isotropic p-extension provides a fast convergence of the finite element approximation. The efficiency can be further improved if an anisotropic discretization is applied which takes into account the special situation due to the plane strain assumptions. Fig. 5 shows
In order to investigate the anisotropic high order finite element formulation for hyperelastic bodies, we consider a numerical example which was presented by Reese and Wriggers [7]. In contrast to the strain-energy function chosen in [7] we apply the one defined by Eq. (1) using the material parameters K = 5000.0 [MPa], a = 0.00367 [MPa], cio =0.1788 [MPa] and CQI =0.1958 [MPa]. The geometry, boundary conditions and finite element mesh of the example are sketched in Fig. 3. The horizontal displacements Ux of the upper surface of the plate are suppressed and the pressure r„ is increased to its maximum value of In =5 [MPa]. Plain strain assumptions are introduced by suppressing the displacements in z-direction at the forefront and backside. Due to synmietry, only one half of the system is meshed with 9 hexahedral elements. One layer of elements in z-direction is applied and the mesh is constructed
-0.5
1".,-1.5 '"' s
<
•
\ \ \
-2
o -2.5 -3
r i ~ D ~ p ~ l 2 ^ 6 -4 ^}; ^^ \ ' ' ^ ' / ^ = p ^ = 1,2,3, ...,6,/7 = 1 - - X - - -
.
\\
n\ \
\\ \ \ \ \ \ \ \ \^
-3.5 -4 ^_ 0
>< 1
•
~~~--^ 1
_
I
. 1
1
1
200 400 600 800 1000 1200 1400 degrees of freedom A^
Fig. 4. Convergence of displacement component Uy of point A.
A. Duster, S. Hartmann /Second MIT Conference on Computational Fluid and Solid Mechanics
241
Y-DISPLACEMENTS _ - 0.72476
I ^-^'^^ P
-0.29463 - -0,80432 • -1.314 * -1.8237 -2.3334 • -2,8431 • -3.3528 • -3.8624
iip^ )eformation(>cl): DISPLACEMENTS of Load_Analysis; step 0.5.
Fig. 5. Deformation of the structure (no scaling) and contour plot of M^ [mm].
the deformation of the structure combined with a contour plot of the vertical displacement. Although the elements are strongly distorted, a locking-free solution is obtained, if the polynomial degree chosen is sufficiently high.
5. Conclusions A high order anisotropic three-dimensional element formulation is applied to hyperelastic continua. Concerning a numerical example with plain strain assumptions, it is numerically demonstrated that the presented anisotropic element formulation leads to a locking-free approximation if a /7-extension is performed. When extending this approach to real three-dimensional problems further advantages are expected. Thin-walled structures like plates and shells can be efficiently discretized, since the polynomial degree in thickness direction may be chosen differently from the one in in-plane direction. This decreases the number of degrees of freedom and therefore reduces the computational time. A pure displacement formulation is sufficient since locking effects arising in thin-walled situations are avoided if the polynomial degree is raised. Future work will include the investigation of this approach for a series of different benchmark problems, focussing especially on three-dimensional thin-walled continua for which the presented element formulation is developed for.
References [1] Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. John Wiley and Sons, Chichester, 2000. [2] Dtister A, Broker H, Rank E. The p-version of the finite element method for three-dimensional curved thin walled structures. Int J Numer Methods Eng 2001;52:673-703. [3] Duster A, Rank E. A p-version finite element approach for two- and three-dimensional problems of the J2flowtheory with non-linear isotropic hardening. Int J Numer Methods Eng 2002;53:49-63. [4] Duster A. High order finite elements for three-dimensional, thin-walled nonlinear continua. PhD thesis, Lehrstuhl ftir Bauinformatik, Technische Universitat Munchen, 2001. Pubfished in: Shaker Verlag, ISBN: 3-8322-0189-0, 2002. [5] Duster A, Niggl A, Rank E. Thermo-elastic computations of geometrically non-linear three-dimensional thin-walled continua based on high order finite elements. Proceedings of the Fifth World Congress on Computational Mechanics, Vienna, Austria 2002. [6] Hartmann S, Neff P. Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for nearincompressibility. Int J Solids Struct 2003; in press. [7] Reese S, Wriggers P. A stabilization technique to avoid hourglassing in finite elasticity. Int J Numer Methods Eng 2000;48:79-109. [8] Suri M. Analytical and computational assessment of locking in the hp finite element method. Comput Methods Appl MechEng 1996;133:347-371. [9] Szabo B, Babuska I. Finite Element Analysis. John Wiley and Sons, Inc., New York, 1991.
High order anisotropic finite elements for three-dimensional isotropic hyperelastic continua A. Duster ^'*, S. Hartmann^ ^Technische Universitdt Miinchen, Fakultdt fur Bauingenieur- und Vermessungswesen, Lehrstuhl fiir Bauinformatik, Arcisstn 21, 80290 Miinchen, Germany ^ Universitdt Kassel, Institut fiir Mechanik, Mdnchebergstr 7, 34109 Kassel, Germany
Abstract In this presentation we apply the p-version finite element method to three-dimensional hyperelastic bodies. On the basis of a pure displacement formulation we take hyperelastic constitutive equations for near incompressibility into account. In order to discretize efficiently three-dimensional continua, an anisotropic element formulation of high order is applied, yielding locking-free approximations for thick as well as thin-walled solids at finite strains. A numerical example demonstrates the basic properties of our nonUnear p-version approach for hyperelastic problems. Keywords: p-Version finite element method; Anisotropic elements; Hyperelasticity; Finite strains
1. Introduction Usually, the finite element method using low order shape functions is applied to finite strain problems. This leads in the case of nearly incompressible material behaviour to volumetric locking phenomena. In order to overcome these problems, a series of mixed finite elements have been proposed (see, for example, Belytschko et al. [1], Reese and Wriggers [7] and the literature cited therein). On the other hand high order finite elements do not show, beyond a certain order of the shape functions, locking effects, which has been investigated by several authors, see for example Suri [8] and Szabo and Babuska [9]. Here, we combine anisotropic p-version finite elements and the finite strain theory to compute three-dimensional continua with nearly incompressible, isotropic hyperelastic material.
deformation gradient. The unimodular left Cauchy-Green tensor B = FF is used to formulate a strain-energy function of the form ^(/,lB,IlB) = f/(/) + u;(%,IlB),
where % = trB and Ilf = trB denote the first and second invariant of B. Due to Hartmann and Neff [6] we assume a strain-energy function V^(/,IB,IIB) with
^ ( / ) = ^ (/' + / ' ' - 2 ) ,
which is polyconvex and coercive, i.e. there exists a solution of the underlying boundary-value problem. The Cauchy stress decomposes in a natural manner into a spherical and a deviatorical part resulting from the specified decomposition of the deformation gradient. T = UXJ)l +
* Tel.: +49 (89) 2892-5060; Fax: +49 (89) 2892-5051; E-mail: [email protected] © 2003 Elsevier Science Ltd. All rights reserved. Computational Fluid and Solid Mechanics 2003 K.J. Bathe (Editor)
(2)
U;(IB,IIB) = Of(II - 27) + ciode - 3) + Coi (II^/' - 3x/3), (3)
2. Constitutive model For the case of the aforementioned investigations, we restrict ourselves to isotropic hyperelasticity. We decompose the deformation gradient F = FF into a volume changing part F = /^/^I and a volume preserving part F = /~^/^F, detF = 1. / = detF symbolizes the determinant of the
(1)
^fe)^ J\6R
(4)
)
where T defines the Cauchy stress tensor and I symbolizes the second order identity tensor. The superscript D denotes the deviator operator, A^ = A — l/3(tr A)I, while trA = a,' defines the trace operator. Fig. 1 shows the basic uniaxial behaviour of rubber materials, i.e. it represents an S-shaped curve in a force-stretch diagram. For more details we refer to [6].
A. Duster, S. Hartmann /Second MIT Conference on Computational Fluid and Solid Mechanics
239
10 c3
<°
2 ^
-10 -15 0.5
1
1.5 2 stretch L/LQ
2.5
Fig. 1. Uniaxal stress-stretch behaviour (material parameters, see Section 4).
3. p-Version FEM at finite deformations The implementation of the /7-version is based on a hexahedral element, allowing for an anisotropic Ansatz of the displacement field (see Fig. 2). The polynomial degree of each separate component of the displacement field can be chosen individually and may also be varied in the three local directions of the element. This anisotropic Ansatz allows the efficient computation of three-dimensional plate and shelllike structures. For a detailed presentation of the element formulation and its application to different linear and nonlinear problems of solid mechanics we refer to [2-5].
Ql = [(-1,1) X (-1,1) X (-1,1)]
The overall solution process is realized by applying an incremental/iterative procedure. Therefore, the nonlinear system of equations arising from the spatial discretization based on the p-version is linearized by means of a Newton-Raphson scheme. The load is applied incrementally, including an adaptation of the load increment which is controlled by the number of iterations of the NewtonRaphson scheme. Since the presented numerical example is highly nonlinear, the adaptation of the load steps was found to be crucial for computing a numerical solution.
S^r'^'^iQllS^^''^'^(QllSP^P^HQl)
Fig. 2. Standard hexahedral Q^^: definition of nodes, edges, faces and polynomial degree.
240
A. Duster, S. Hartmann / Second MIT Conference on Computational Fluid and Solid Mechanics
t t t t
!-V " " "^
point A t
t = 1 [mm]
/W;c=0
/
n/
\
/
\
.iL^.^i
5 [mm],
G
10 [mm]
I
5 [mm] I
Fig. 3. Inhomogeneous compression test (geometry and element mesh).
4. Numerical examples
in such a way that the local coordinate f of all elements coincides with the z-axis. The convergence of the displacement component M^ of point A is plotted in Fig. 4 for both an isotropic and an anisotropic increase of the polynomial degree pk,k = ^,r],^ based on the trunk space 5tf'^"''^(^^t) (for the definition of the trunk space see [2,4]). In the case of the isotropic discretization the polynomial degree of the elements is increased in all three local directions (|,y;,f) simultaneously, whereas in the case of the anisotropic discretization the polynomial degree is only increased in in-plane direction (§,/;) and is kept fixed at p^ = 1 in f-direction in order to account for the plain strain assumptions. From this it is obvious that the isotropic p-extension provides a fast convergence of the finite element approximation. The efficiency can be further improved if an anisotropic discretization is applied which takes into account the special situation due to the plane strain assumptions. Fig. 5 shows
In order to investigate the anisotropic high order finite element formulation for hyperelastic bodies, we consider a numerical example which was presented by Reese and Wriggers [7]. In contrast to the strain-energy function chosen in [7] we apply the one defined by Eq. (1) using the material parameters K = 5000.0 [MPa], a = 0.00367 [MPa], cio =0.1788 [MPa] and CQI =0.1958 [MPa]. The geometry, boundary conditions and finite element mesh of the example are sketched in Fig. 3. The horizontal displacements Ux of the upper surface of the plate are suppressed and the pressure r„ is increased to its maximum value of In =5 [MPa]. Plain strain assumptions are introduced by suppressing the displacements in z-direction at the forefront and backside. Due to synmietry, only one half of the system is meshed with 9 hexahedral elements. One layer of elements in z-direction is applied and the mesh is constructed
-0.5
1".,-1.5 '"' s
<
•
\ \ \
-2
o -2.5 -3
r i ~ D ~ p ~ l 2 ^ 6 -4 ^}; ^^ \ ' ' ^ ' / ^ = p ^ = 1,2,3, ...,6,/7 = 1 - - X - - -
.
\\
n\ \
\\ \ \ \ \ \ \ \ \^
-3.5 -4 ^_ 0
>< 1
•
~~~--^ 1
_
I
. 1
1
1
200 400 600 800 1000 1200 1400 degrees of freedom A^
Fig. 4. Convergence of displacement component Uy of point A.
A. Duster, S. Hartmann /Second MIT Conference on Computational Fluid and Solid Mechanics
241
Y-DISPLACEMENTS _ - 0.72476
I ^-^'^^ P
-0.29463 - -0,80432 • -1.314 * -1.8237 -2.3334 • -2,8431 • -3.3528 • -3.8624
iip^ )eformation(>cl): DISPLACEMENTS of Load_Analysis; step 0.5.
Fig. 5. Deformation of the structure (no scaling) and contour plot of M^ [mm].
the deformation of the structure combined with a contour plot of the vertical displacement. Although the elements are strongly distorted, a locking-free solution is obtained, if the polynomial degree chosen is sufficiently high.
5. Conclusions A high order anisotropic three-dimensional element formulation is applied to hyperelastic continua. Concerning a numerical example with plain strain assumptions, it is numerically demonstrated that the presented anisotropic element formulation leads to a locking-free approximation if a /7-extension is performed. When extending this approach to real three-dimensional problems further advantages are expected. Thin-walled structures like plates and shells can be efficiently discretized, since the polynomial degree in thickness direction may be chosen differently from the one in in-plane direction. This decreases the number of degrees of freedom and therefore reduces the computational time. A pure displacement formulation is sufficient since locking effects arising in thin-walled situations are avoided if the polynomial degree is raised. Future work will include the investigation of this approach for a series of different benchmark problems, focussing especially on three-dimensional thin-walled continua for which the presented element formulation is developed for.
References [1] Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. John Wiley and Sons, Chichester, 2000. [2] Dtister A, Broker H, Rank E. The p-version of the finite element method for three-dimensional curved thin walled structures. Int J Numer Methods Eng 2001;52:673-703. [3] Duster A, Rank E. A p-version finite element approach for two- and three-dimensional problems of the J2flowtheory with non-linear isotropic hardening. Int J Numer Methods Eng 2002;53:49-63. [4] Duster A. High order finite elements for three-dimensional, thin-walled nonlinear continua. PhD thesis, Lehrstuhl ftir Bauinformatik, Technische Universitat Munchen, 2001. Pubfished in: Shaker Verlag, ISBN: 3-8322-0189-0, 2002. [5] Duster A, Niggl A, Rank E. Thermo-elastic computations of geometrically non-linear three-dimensional thin-walled continua based on high order finite elements. Proceedings of the Fifth World Congress on Computational Mechanics, Vienna, Austria 2002. [6] Hartmann S, Neff P. Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for nearincompressibility. Int J Solids Struct 2003; in press. [7] Reese S, Wriggers P. A stabilization technique to avoid hourglassing in finite elasticity. Int J Numer Methods Eng 2000;48:79-109. [8] Suri M. Analytical and computational assessment of locking in the hp finite element method. Comput Methods Appl MechEng 1996;133:347-371. [9] Szabo B, Babuska I. Finite Element Analysis. John Wiley and Sons, Inc., New York, 1991.