A computational framework for polyconvex large strain elasticity

TL;DR: In this paper, the deformation gradient (the fibre map), its adjoint (the area map) and its determinant (the volume map) are introduced as independent kinematic variables of a convex strain energy function.

About: This article is published in Computer Methods in Applied Mechanics and Engineering.The article was published on 2015-01-01 and is currently open access. It has received 115 citations till now. The article focuses on the topics: Interpolation & Piecewise.

Summary

1. Introduction

• Large strain elastic and inelastic analysis by finite elements or other computational techniques is now well established for many engineering applications [1–8].
• The resulting formulation opens up new interesting possibilities in terms of using various interpolation spaces for different variables [28–31], leading to enhanced type of formulations [24].
• This new notation not only greatly simplifies the formulation but it also provides new expressions for operators such as the tangent elasticity tensor which are useful from a computational and theoretical point of view.
• This definition makes it possible to introduce HellingerReissner [34, 39] type of mixed variational principles in the context of large strain elasticity, with a significantly reduced set of variables over a more traditional Hu-Washizu type of functional.
• Section 3 reviews the definition of polyconvex elastic strain energy functions and defines a new set of stresses conjugate to the main kinematic variables.

2.1. Motion and deformation

• The standard notation and definitions for the deformation gradient and its determinant are used: F = ∂x ∂X = ∇0x; J = detF = dv dV (1) where x represents the current position of a particle originally at X and ∇0 denotes the gradient with respect to material coordinates.
• (2) Clearly, the components of this tensor are the order 2 minors of the deformation gradient and it is often referred to as the co-factor or adjoint tensor.
• This tensor and its derivatives will feature heavily in the formulation that follows as it is a key variable for polyconvex elastic models.
• Its evaluation and, more importantly, the evaluation of its derivatives using equation (2) is not ideal, and a more convenient formula can be derived for 3-dimensional applications.
• This relies on the definition of the tensor cross product given in the next section.

2.2. Tensor cross product

• One of the key elements of the framework proposed is the extension of the standard vector cross product to define the cross product between second order tensors and between tensors and vectors.
• For this purpose the above definition can be readily particularised to second order two-point tensors or material tensors as, (A B)iI = EijkEIJKAjJBkK ; (A B)IJ = EIKLEJMNAKMBLN .
• (6) Box 1 shows the practical evaluation of these products.

2.3. Alternative expressions for the area and volume maps

• The above formulas simplify the manipulation of the derivatives ofH and J by avoiding differentiating the inverse of the deformation gradient.
• They will be key to the development of the framework presented below.

3.1. The strain energy

• Polyconvexity is now well accepted as a fundamental mathematical requirement that must be satisfied by admissible strain energy functions used to describe elastic materials in the large strain regime.
• Doing this, however, has no practical effect on the resulting formulation as this will be driven by derivatives of the strain energy.
• Appropriate values for α and β and suitable functions f will be found in the sections below.

3.4. Tangent elasticity operator

• This leaves only the contribution from the first positive definite term in equation (48).
• It is therefore easy to note that polyconvexity implies ellipticity [13].

3.6. The Second Piola-Kirchhoff, Kirchhoff and Cauchy stress tensors

• Note, however, that the function W̃ need not to be strictly convex with respect to its variables.
• This derivation will not be pursued here as it will not be used in the computational framework proposed later.
• In addition to the first and second Piola-Kirchhoff stresses, it is necessary to derive expressions for the Cauchy and Kirchhoff stresses as often these tensors are needed in order to express plasticity models or simply to display solution results.

4. Variational formulations

• This section presents several possible variational formulations.
• The section starts reviewing the standard displacement based variational principle.
• For this case the present framework does not provide any practical advantages, other than the evaluation of the first Piola-Kirchhoff stress and the tangent modulus via equations (39) and (48), respectively.
• It is simply presented here in order to provide a useful background for comparison with mixed and complementary energy variational principles presented later in the section.
• The aim of the section is to present the concepts in as simple a manner as possible, rather than to provide precise mathematical statements given that, in practice, the concepts presented here will be implemented in well defined discrete finite element spaces described in the following section.

4.1. Standard displacement based variational principle

• The solution of large strain elastic problems is often expressed by means of the total energy minimisation variational principle as: Π(x∗) = min x (83) where x∗ denotes the exact solution.
• The first Piola-Kirchhoff tensor P x is evaluated in the standard fashion using equation (39) in terms of the gradient of the deformation ∇0x.

4.2. Mixed Variational Principle

• It is important to note that for, certain choices of interpolation spaces insufficiently rich, the above variational formulation may not ensure that the kinematic compatibility constraints are enforced at every point in the domain.
• There are, however, a number of cross derivative terms that do not vanish.

4.3. Complementary Energy Principle

• . (104) Comparing the term in square brackets in the first integral with the definition of the complementary energy given by equation (42), enables a com- plementary variational principle to be established as3: ΠM(x ∗,ΣF ∗,Σ∗ H ,Σ∗J) = min x∈H max ΣF ,ΣH ,ΣJ .
• (105) This represents a Helinger-Reissner type of variational principle [1].
• D21;1ΠC [δv;u] =D 2 1;1ΠM [δv;u] (108) A new term, however, emerges when taking the second derivatives with respect to stresses, leading to a constitutive expression similar to equation (100) but now involving the second derivative of the complementary energy function, also known as lowing terms are identical.
• This would also change the sign of D21;1ΠC [δv;u] but this term is neither positive nor negative since it contains tensor cross products of the gradient of δv and u.
• For practical purposes, however, this will not be necessary, as the stress variables will typically be eliminated locally in each finite element.

4.4. Variational principles for incompressible and nearly incompressible models

• Many applications of practical importance rely on the decomposition of the strain energy into isochoric and volumetric components.
• For such cases, it is possible to modify the variational formulations above in such a way that different approaches are used for the isochoric and volumetric components.
• The third integral term above enforces the compatibility between these two measures. .
• (111) The stationary conditions of these hybrid functionals are evaluated in the same fashion as above.
• The evaluation of second derivatives required for a Newton-Raphson process proceeds along the same lines as in the previous sections.

5. Finite Element implementation

• It is not the purpose of this paper to provide an exhaustive analysis of different finite element interpolation strategies to discretise the equations derived in previous section.
• Two particular examples will be provided in order to demonstrate the validity of the formulation presented.

5.1. General remarks

• The implementation of the various variational principles described in the previous section is based on a finite element partition of the domain into a set of elements.
• In general, different interpolations can (and are often) used to describe different variables.
• In order to complete the finite element formulation it is necessary to derive equations for the components of the tangent matrix by discretising the tangent operators defined in the previous section.
• (128) Some of the terms in the above matrices are very straightforward to obtain.
• For instance the term KDΣ relating stresses and strains in equation (127) follows from the discretisation of the corresponding tangent operator component given in equation (101) as: Kab DΣ = − ∫ V dV (129) where I denotes the components of the fourth order identity tensor.

5.2. Complementary energy application case

• The equations provided above can be implemented using a variety of finite element spaces.
• Of course, not all choices will lead to effective or valid finite element formulations.
• Stress and strain variables can actually be discretised independently on each element of the mesh.
• This enables a static condensation process to be carried out before assembly of the global tangent matrix.
• The resulting element is very similar to that proposed in reference [24].

5.3. Stabilised linear tetrahedron for incompressible elasticity

• As a final case study of the application of the present framework, a linear tetrahedron for modelling full incompressible elasticity is described.
• (143) Note that in this case the pressure is continuous across elements and therefore cannot be eliminated locally.
• Unfortunately, it is well know that an incompressible mixed formulation with linear displacements and linear pressures on tetrahedral element spaces does not satisfy the necessary LBB condition and leads to unstable solutions [3].
• The classical solution to this problem is to introduce Petrov-Galerkin stabilisation [46–49].

6. Numerical examples

• The objective of this section is to present a series of numerical examples in order to prove the robustness, accuracy and applicability of the computational framework presented above.
• All of the numerical results presented correspond to the following selection of functional spaces: continuous quadratic interpolation of the displacement field x, piecewise linear interpolation of the strain and stress fields F , H , ΣF and ΣH and piecewise constant interpolation of the Jacobian J and its associated stress conjugate ΣJ .
• With these functional spaces, the three mixed formulationsM7F,MCF andM5F will render identical results.
• For the incompressible case, a two field {x, p} mixed formulation is employed, as described in Section 4.4.

6.1. Patch test

• The first numerical example includes a standard three dimensional patch test in order to assess the correctness of the computational implementation.
• This problem was already presented in [24].
• (167) For both constitutive models Ws and Wq, the linear elasticity constitutive operator in the reference configuration renders the same material parameters, namely shear modulus µ = 756kPa and Poisson’s ratio ν = 0.49546.
• To achieve this deformation, non-zero normal Dirichlet boundary conditions are applied on the boundary faces perpendicular to the OX axis and zero normal Dirichlet boundary conditions are defined on two adjacent faces perpendicular to the OY and OZ axes.
• As expected, hence passing the patch test, for the three mixed formulations defined above, namely M7F, MCF and M5F, the results are identical for both meshes (the only numerical difference due to machine accuracy).

6.2. Cook type cantilever problem.

• The geometry of the problem is shown in Figure 4(a) where, as it can be observed, the cantilever is clamped on its left end and subjected to an upwards parabolic shear force distribution applied on its right end of maximum value τmax = 16kPa.
• Figures 5 to 7 display the contour plot of different stress magnitudes (σxx, σyy and pressure) using the constitutive model Wp and the fine discretisa- tion.
• Results are presented comparing the mixed formulations, which render identical results , 6(a) and 7(a)) versus the DF formulation , 6(b) and 7(b)).
• The results obtained with the constitutive model Wq match very well those presented in reference [24].

6.3. Compressible short column subjected to transverse shear force.

• Note that this force is not considered to be a follower-load during the deformation process.
• As a closed form solution is not available for this problem, the finest mesh is used to generate numerically the so-called “benchmark” solution (for each mixed formulation) for comparison purposes.
• Figure 9 displays the contour plot distribution of the stress σxx for different stages of the deformation.
• The results presented are identical for all of the mixed formulations.
• Figure 10 shows the order of accuracy of the different unknown variables for the mixed formulations (all yielding identical convergence pattern).

6.4. Incompressible long column subjected to transverse shear force.

• As in the previous section, the size of the square defining the cross section remains a = 1m.
• Note that this force is not considered to be a follower-load during the deformation process.
• The selection of these material parameters yield a linear elasticity constitutive operator in the origin defined by shear modulus µ = 756kPa and Poisson’s ratio ν = 0.5.
• For this purpose, two interpolation techniques will be considered: first, a P1-P1 linear continuous interpolation for both displacement and pressure fields (e.g. with the help of stabilisation) and, second, a P2-P0 continuous quadratic interpolation for the displacement field and piecewise constant interpolation for the pressure field (e.g. without stabilisation).
• The results show an excellent agreement between the non-stabilised P2-P0 formulation and 13(d)) and the stabilised P1-P1 formulation when using α = 0.2 and 13(b)).

6.5. Twisting cantilever beam

• The beam is clamped at its left end and subjected to a torsion on its right end.
• As can be observed, the section is not restricted to in-plane torsion and zero Neumann boundary conditions are imposed normal to the cross sectional area.
• A similar example has been presented by the authors in previous references [47, 49].
• The analysis is carried out with a finite element discretisation comprised of 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).
• Results displayed in Figures 15(b) and 15(d) are obtained by using the DF implementation whereas results in Figures 15(a) and 15(c) are obtained by using the alternative mixed formulations, namely M7F, MCF and M5F, all yielding identical results.

7. Concluding remarks

• This paper has provided a novel approach to formulate large polyconvex large strain elasticity in the computational context.
• The key novel contributions of the work presented here are: .
• The definition of stresses {ΣF ,ΣH ,ΣJ} conjugate to the main extended kinematic variable set {F ,H , J} and the development of relationships between these stresses and more classical stress tensors such as the Piola-Kirchhoff stress.
• The formulation of several mixed variational principles including novel Hellinger-Reissner two field versions which reduce the number of problem variables over more traditional three field Hu-Washizu type of principles.
• Moreover, the use of discretisation spaces in which deformation or stress fields are discontinuous across element faces enables these variables to be resolved locally leading to computational cost very comparable to those of displacement based approaches.

Figures

Figure 8: Compressible short column of height h = 5m and squared cross section defined by its size a = 1m. (a) Boundary conditions: clamped bottom end and parabolic stress distribution at the top end. Axes OX1 and OX2 coincide with oy and oz in (b), respectively. (b) Example of finite element discretisation: coarsest mesh with (2 × 2 × 10) × 6 tetrahedral elements.

Figure 4: Cook type cantilever problem. (a) Geometry of the problem and boundary conditions, where a parabolic upwards shear stress distribution with maximum value τmax = 300kPa is applied on its right end. Axes OX1 and OX2 coincide with ox and oy in (b), respectively. (b) Representation of the finest discretisation employed with (14× 14× 5) × 6 tetrahedral elements. (c) Contour plot of displacement field uy(m) for the polyconvex model Wq defined in (165) with material parameters in (167). (d) Contour plot of stress field σxx(Pa) for the polyconvex model Wq defined in (165) with material parameters in (167).

Figure 14: Twisting cantilever beam problem of length h = 6m and squared cross section defined by its size a = 1m. (a) Boundary conditions: clamped left end and twisting rotation applied at the right end. Axes OX1, OX2 and OX3 coincide with ox, oy and oz in (b), respectively. (b) Example of finite element discretisation: 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).

Figure 18: Contour plot of the distribution of the components (a) H22, (b) H21, (c) H32, and (d) J , for the twisting cantilever beam example. Results obtained with the alternative mixed formulations. Mooney-Rivlin model Wu defined in (179) with material parameters given in (180). Results shown for a discretisation of 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).

Figure 9: Contour plot of the stress σxx(Pa) distribution for the compressible short column example using a mixed formulation. Results after application of an incremental loading of (a) 5% (b) 25% (c) 50% and (d) 100% of the total external shear stress τ = τmax. MooneyRivlin model Wp defined in (169) with material parameters given in (174). Results shown for a discretisation of (4× 4× 20)× 6 tetrahedral elements (3, 321× 3 degrees of freedom associated to the spatial coordinates x).

Figure 3: Three dimensional patch test. Quadratic convergence of the Newton-Raphson linearisation procedure.

Figure 2: Three dimensional patch test. (a) View of half undistorted mesh in the reference configuration. (b) View of half distorted mesh in the reference configuration. (c) Example of deformed geometry after stretching of ∆L/L = 0.5 in the OX direction for a MooneyRivlin model Ws defined in (163) with parameters defined in (164).

Table 2: Cook type cantilever problem. Stress components (kPa) and displacements (m) at points A, B and C as displayed in Figure 4(b). Results are obtained using the mixed formulations. Coarse, medium and fine discretisations of (7× 7× 5)×6, (10× 10× 6)×6 and (14× 14× 5)×6 tetrahedral elements, respectively. Wp model (columns 2 to 4) defined in (169) with parameters given in (170). Wq model (columns 5 to 7) defined in (165) with material parameters in (167).

Table 3: Cook type cantilever problem. Stress components (kPa) and displacements (m) at points A, B and C as displayed in Figure 4(b). Results obtained using the DF formulation. Coarse, medium and fine discretisations of (7× 7× 5)× 6, (10× 10× 6)× 6 and (14× 14× 5)×6 tetrahedral elements, respectively. Wp model (columns 2 to 4) defined in (169) with parameters given in (170). Wq model (columns 5 to 7) defined in (165) with material parameters in (167).

Figure 17: Contour plot of the distribution of the components (a) F11, (b) F12, (c) F32, (d) H12, for the twisting cantilever beam example. Results obtained with the alternative mixed formulations. Mooney-Rivlin model Wu defined in (179) with material parameters given in (180). Results shown for a discretisation of 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).

Figure 13: Contour plot of the displacement ux(m) distribution for the incompressible long column example. (a) Stabilised P1-P1 formulation with τp = 0 and α = 0.1. (b) Stabilised P1-P1 formulation with τp = 0 and α = 0.2. (c) Stabilised P1-P1 formulation with τp = 0 and α = 0.5. (d) Non-stabilised P2-P0 formulation. Incompresible Mooney-Rivlin model Ŵi defined in (176) with material parameters given in (177). Results shown in (a) to (c) are for a discretisation of 17,575 tetrahedral elements (3, 859 × 3 degrees of freedom associated to the spatial coordinates x). Results shown in (d) are for a discretisation of 3,962 tetrahedral elements (6, 675 × 3 degrees of freedom associated to the spatial coordinates x).

Figure 10: Compressible short column example: order of accuracy of different strain and stress magnitudes for the mixed formulations. (a) Order of accuracy of the kinematic variables x, F , H and J . (b) Order of accuracy of the kinetic variables ΣF , ΣH and ΣJ . Mooney-Rivlin model Wp defined in (169) with material parameters given in (174).

Figure 7: Cook type cantilever problem. Contour plot of the distribution of the hydrostatic pressure p(Pa) for (a) mixed formulations and (b) DF formulation. Mooney-Rivlin model Wp defined in (169) with material parameters given in (170). Discretisation of (14× 14× 5)× 6 tetrahedral elements.

Figure 6: Cook type cantilever problem. Contour plot of the distribution of stress coefficient σyy(Pa) for (a) mixed formulations and (b) DF formulation. Mooney-Rivlin model Wp defined in (169) with material parameters given in (170). Discretisation of (14× 14× 5)× 6 tetrahedral elements.

Figure 16: Lateral and plan view of the contour plot of the hydrostatic pressure distribution p(Pa) for the twisting cantilever beam example. Results obtained with the alternative mixed formulations. Mooney-Rivlin model Wu defined in (179) with material parameters given in (180). Results shown for a discretisation of 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).

Table 1: Cook type cantilever problem. Mesh discretisation details. Column 2: number of tetrahedral elements (Elems.). Column 3: number of degrees of freedom (Dofs.) associated to the spatial coordinates x. Column 4: number of degrees of freedom (Dofs.) associated to the strain/stress fields F ,H,ΣF ,ΣH . Column 5: number of degrees of freedom (Dofs.) associated to the strain/stress fields J,ΣJ .

Figure 11: Compressible short column example: quadratic convergence of the NewtonRaphson linearisation procedure. Mooney-Rivlin model Wp defined in (169) with material parameters given in (174).

Figure 12: Contour plot of the hydrostatic pressure distribution p(kPa) for the incompressible long column example. (a) Stabilised P1-P1 formulation with τp = 0 and α = 0.1. (b) Stabilised P1-P1 formulation with τp = 0 and α = 0.2. (c) Stabilised P1-P1 formulation with τp = 0 and α = 0.5. (d) Non-stabilised P2-P0 formulation. Incompresible Mooney-Rivlin model Ŵi defined in (176) with material parameters given in (177). Results shown in (a) to (c) are for a discretisation of 17,575 tetrahedral elements (3, 859× 3 degrees of freedom associated to the spatial coordinates x). Results shown in (d) are for a discretisation of 3,962 tetrahedral elements (6, 675× 3 degrees of freedom associated to the spatial coordinates x).

Figure 1: Deformation mapping of a continuum and associated kinematics magnitudes: F ,H, J .

Figure 5: Cook type cantilever problem. Contour plot of the distribution of stress coefficient σxx(Pa) for (a) mixed formulations and (b) DF formulation. Mooney-Rivlin model Wp defined in (169) with material parameters given in (170). Discretisation of (14× 14× 5)× 6 tetrahedral elements.

Figure 15: Contour plot of the stress σxy(Pa) ((a) and (b)) and σyy(Pa) ((c) and (d)) distribution for the twisting cantilever beam example. Results obtained with the DF formulation ((b) and (d)) and alternative mixed formulations ((a) and (c)). MooneyRivlin model Wu defined in (179) with material parameters given in (180). Results shown for a discretisation of 2,304 tetrahedral elements (4, 009× 3 degrees of freedom associated to the spatial coordinates x).

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "A computational framework for polyconvex large strain elasticity" ?

This paper presents a new computational formulation for large strain polyconvex elasticity. The paper also introduces conjugate stresses to these kinematic variables which can be used to define a generalised convex complementary energy function and a corresponding complementary energy principle of the Hellinger-Reissner type, where the new conjugate stresses are primary variables together with the deformed geometry. A key element to the developments presented in the paper is the definition of a new tensor cross product which facilitates the algebra associated with the adjoint of the deformation gradient. 

Q2. What are the future works mentioned in the paper "A computational framework for polyconvex large strain elasticity" ?

The definition of stresses { ΣF, ΣH, ΣJ } conjugate to the main extended kinematic variable set { F, H, J } and the development of relationships between these stresses and more classical stress tensors such as the Piola-Kirchhoff stress. Future work will consider the extension of the present framework to solid dynamics using appropriate conservation laws for the extended set of kinematic variables in the manner proposed in references [ 47, 49–54 ].