1 Introduction

Among different strain measures in the finite-deformation theory, the logarithmic (Hencky) strain possesses special properties that have led to a variety of applications in solid mechanics, e.g. [1,2,3,4,5]. These, for instance, include the fully uncoupled additive decomposition of deformation into the volumetric and isochoric parts, or straightforward characterization of the incompressibility constraint. Moreover, the elastic strain energy defined as a quadratic function of the Hencky strain (in short, Hencky strain energy) behaves reasonably well for a substantially wider range of elastic strains compared to the simple St. Venant-Kirchhoff strain energy and provides a physically more relevant material response at finite strains [6, 7]. For an overview of the special properties of the Hencky strain energy, see [8].

The main difficulty related to the use of the Hencky strain energy is that, in addition to computation of the logarithm of a \(3\times 3\) matrix, the first and second derivatives of the matrix logarithm are required in order to derive the stress and the tangent operator, respectively. The most common approaches to represent the matrix logarithm along with its derivatives are based on the spectral decomposition, e.g. [9], or series expansion, e.g. [10]. The former requires evaluation of the eigenvalues and eigenvectors, which may pose numerical challenges in the vicinity of repeated eigenvalues, while the latter is convergent only when the matrix is close to identity, and lacks computational efficiency when high-order approximations are used, see [11] for an overview regarding the common approaches to compute general matrix functions, see also [12].

The primary motivation for the present study comes from our work on finite-strain phase-field modelling of martensitic phase transformation in shape memory alloys (SMAs), e.g. [13]. In the related models, elastic strain energy functions are needed for arbitrary anisotropy, e.g. within the diffuse austenite–martensite interfaces. The St. Venant-Kirchhoff strain energy is a common choice, as e.g. in [14], but it is known to behave poorly in some situations, for instance, under compression (this is also illustrated in the present paper). Hence, the Hencky strain energy is an attractive alternative, as it behaves well for a much wider range of strains [7, 15]. The Hencky strain energy may also be a reasonable replacement for the anisotropic St. Venant-Kirchhoff strain energy in other application areas, for instance, in crystal plasticity [16, 17], and in micromechanics of SMA single crystals [18, 19].

In the present study, the focus is on the use of the rational approximation method, known as Padé approximants, for evaluation of the tensor logarithm (or matrix logarithm, depending on the context). Padé approximants have proved to be more accurate with respect to the classical Taylor series, possess a wider range of applicability and converge to the solution more rapidly as the order of approximation increases, see [20] for a general description of Padé approximation method. The associated derivatives of the matrix logarithm, needed in practical applications in computational mechanics, are then computed straightforwardly by differentiating the corresponding Padé approximant which is available in an explicit form.

The Padé approximation method is actually a general approach that can be used to approximate arbitrary tensor functions. It seems that the only related applications in computational mechanics concern the tensor logarithm [21, 22] and the tensor exponential [21,22,23,24]. In particular, the approximated Hencky strain was used by Brünig [21, 22] to express the Hencky strain energy in isotropic elasticity.

In this paper, we examine the applicability of Padé approximants in finite-element modelling of Hencky-type hyperelasticity problems with an emphasis on the anisotropic elasticity. For this purpose, the performance of Padé approximants is evaluated against the truncated Taylor series, and the closed-from representation of the matrix logarithm [11] is considered as a reference for the comparison.

A general description of the Padé approximation method is presented in Sect. 2. In Sect. 3, we discuss the Hencky-type constitutive models of hyperelasticity. A side remark on the stress conjugate to the Hencky strain is provided in Sect. 4. Finally, performance of Padé approximants in finite-element computations is examined in Sect. 5.

2 Padé approximants

2.1 Padé approximants of scalar and tensor functions

The Padé approximant \(R_{(m,n)}(x)\) of order (m, n) of a scalar function f(x) is defined as a rational polynomial,

$$\begin{aligned} f(x) \approx R_{(m,n)}(x) = \frac{N_{(m,n)}(x)}{D_{(m,n)}(x)} , \end{aligned}$$
(1)

where

$$\begin{aligned} N_{(m,n)}(x) = \sum _{j=0}^m a_j x^j , \quad D_{(m,n)}(x) = 1+\sum _{k=1}^n b_k x^k , \end{aligned}$$
(2)

and the coefficients \(a_j\) and \(b_k\) are uniquely determined such that the Maclaurin series expansion of the Padé approximant \(R_{(m,n)}(x)\) agrees with that of the function f(x) up to the order of \(m+n\). Note that \(D_{(m,n)}\) has been normalized such that the zeroth-order term in \(D_{(m,n)}\) is equal to unity, \(b_0=1\), and \(N_{(m,n)}\) and \(D_{(m,n)}\) have no common factors [20]. Both \(N_{(m,n)}\) and \(D_{(m,n)}\) depend on the approximation order (m, n), hence the dual index is used to indicate this dependence.

Consider now an isotropic tensor-valued function \(f(\mathrm{\mathbf{X}})\) of a tensor argument \(\mathrm{\mathbf{X}}\), where \(\mathrm{\mathbf{X}}\) is a second-order tensor in a three-dimensional vector space. For a symmetric tensor \(\mathrm{\mathbf{X}}\), function \(f(\mathrm{\mathbf{X}})\) can be defined using the spectral representation, viz.

$$\begin{aligned} \mathrm{\mathbf{X}} = \sum _{i=1}^3 \lambda _i {{\varvec{\omega }}_i \otimes {\varvec{\omega }}_i} , \quad f(\mathrm{\mathbf{X}}) = \sum _{i=1}^3 f(\lambda _i) {{\varvec{\omega }}_i \otimes {\varvec{\omega }}_i} , \end{aligned}$$
(3)

where \(\lambda _i\) and \({\varvec{\omega }}_i\) are, respectively, the eigenvalues and eigenvectors of \(\mathrm{\mathbf{X}}\). Alternatively, the Taylor series expansion can be used to define \(f(\mathrm{\mathbf{X}})\), which is applicable also to non-symmetric \(\mathrm{\mathbf{X}}\). Here and below, we refer to a tensor function, but the same argument applies to the respective matrix function, where the matrix can be treated as a representation of the tensor in a given orthonormal basis.

In analogy to the scalar function case, the Padé approximant \(R_{(m,n)}(\mathrm{\mathbf{X}})\) of order (m, n) of function \(f(\mathrm{\mathbf{X}})\) has the form

$$\begin{aligned} f(\mathrm{\mathbf{X}}) \approx R_{(m,n)}(\mathrm{\mathbf{X}}) = N_{(m,n)}(\mathrm{\mathbf{X}}) (D_{(m,n)}(\mathrm{\mathbf{X}}))^{-1} , \end{aligned}$$
(4)

where \(N_{(m,n)}\) and \(D_{(m,n)}\) are the polynomials of the Padé approximant \(R_{(m,n)}(x)\) of the corresponding scalar function f(x). Clearly, \(N_{(m,n)}(\mathrm{\mathbf{X}})\) and \(D_{(m,n)}(\mathrm{\mathbf{X}})\) are here tensor polynomials of the form (2) with the scalar argument x replaced by the tensor \(\mathrm{\mathbf{X}}\), with no restriction on the symmetry of X. As in the scalar case, the approximation is here performed at \(\mathrm{\mathbf{X}}=\mathrm{\mathbf{0}}\). Considering that \(N_{(m,n)}(\mathrm{\mathbf{X}})\) and \(D_{(m,n)}(\mathrm{\mathbf{X}})\) commute, the order in Eq. (4) can be changed, so that \(R_{(m,n)}(\mathrm{\mathbf{X}}) = (D_{(m,n)}(\mathrm{\mathbf{X}}))^{-1} N_{(m,n)}(\mathrm{\mathbf{X}})\).

As in the case of the original tensor function \(f(\mathrm{\mathbf{X}})\) with a symmetric argument \(\mathrm{\mathbf{X}}\), the Padé approximation of \(f(\mathrm{\mathbf{X}})\) amounts to applying the Padé approximant to the eigenvalues in the spectral decomposition, thus

$$\begin{aligned} f(x) \approx R_{(m,n)}(x) , \quad f(\mathrm{\mathbf{X}}) \approx \sum _{i=1}^3 R_{(m,n)}(\lambda _i) {{\varvec{\omega }}_i} \otimes {{\varvec{\omega }}_i} . \end{aligned}$$
(5)

Note that the above representation is here provided for illustration purposes only, and the rational polynomial representation (4) would always be preferable as it does not require spectral decomposition. Actually, the spectral decomposition allows one to compute the exact tensor function according to Eq. (3) so that the approximation can be avoided.

It has been proved [25] that for a fixed order of Padé approximant (\(m+n=\text {constant}\)), the approximation error is the smallest for diagonal Padé approximants, i.e. when \(m=n\). Accordingly, the focus of the present study is restricted to diagonal Padé approximants.

2.2 Padé approximation of the logarithm function

Consider now the logarithm function and apply the Padé approximation technique to this function. With reference to the target application in continuum mechanics, the approximation of \(\log (x)\) is sought in the vicinity of \(x=1\). Upon setting \(x'=x-1\), the approximation of \(\log (x)\) at \(x=1\) is

$$\begin{aligned} \log (1+x') \approx R_{(m,n)}(x') = \frac{N_{(m,n)}(x')}{D_{(m,n)}(x')}. \end{aligned}$$
(6)

The corresponding approximation of the function \(\log (x)\) can then be readily obtained in the following form,

$$\begin{aligned} \log (x) \approx R^*_{(m,n)}(x) = \frac{N^*_{(m,n)}(x)}{D^*_{(m,n)}(x)} , \end{aligned}$$
(7)

where

$$\begin{aligned} N^*_{(m,n)}(x) = N_{(m,n)}(x-1) , \quad D^*_{(m,n)}(x) = D_{(m,n)}(x-1) . \end{aligned}$$
(8)

The respective polynomials of the diagonal Padé approximants of orders (1, 1) through (5, 5) are provided in Table 1.

Table 1 Polynomials \(N^*_{(m,n)}(x)\) and \(D^*_{(m,n)}(x)\) of selected diagonal Padé approximants of \(\log (x)\)
Full size table

Clearly, following Sect.  2.1, the Padé approximants of the tensor logarithm can be obtained by replacing the scalar polynomials \(N^*_{(m,n)}\) and \(D^*_{(m,n)}\) in Eq. (7) by the corresponding tensor polynomials in terms of \(\mathrm{\mathbf{X}}\). For instance, the approximants of order (1, 1) and (2, 2) are given by

$$\begin{aligned}&{\log (\mathrm{\mathbf{X}}) \approx R^*_{(1,1)}(\mathrm{\mathbf{X}}) = 2(\mathrm{\mathbf{X}}-\mathrm{\mathbf{1}})(\mathrm{\mathbf{X}}+\mathrm{\mathbf{1}})^{-1} ,} \end{aligned}$$
(9)
$$\begin{aligned}&{\log (\mathrm{\mathbf{X}}) \approx R^*_{(2,2)}(\mathrm{\mathbf{X}}) = 3(\mathrm{\mathbf{X}}^2-\mathrm{\mathbf{1}})(\mathrm{\mathbf{X}}^2+4\mathrm{\mathbf{X}}+\mathrm{\mathbf{1}})^{-1},} \end{aligned}$$
(10)

where \(\mathrm{\mathbf{1}}\) denotes the second-order identity tensor.

An alternative method of approximating the tensor logarithm employs the truncated Taylor series. The Taylor series expansion of the scalar logarithm function \(f(x)=\log (x)\) evaluated at \(x=1\) and truncated to order m reads

$$\begin{aligned} \log (x) \approx T_m(x) = \sum _{k=1}^m (-1)^{k+1}\frac{(x-1)^k}{k} , \end{aligned}$$
(11)

which is convergent when \(m\rightarrow \infty \) for \(0<x\le 2\). The truncated Taylor series approximation of the tensor logarithm is thus given by

$$\begin{aligned} \log (\mathrm{\mathbf{X}})\approx & {} T_m(\mathrm{\mathbf{X}}) {}= \mathrm{\mathbf{X}}-\mathrm{\mathbf{1}} - \frac{1}{2} (\mathrm{\mathbf{X}}-\mathrm{\mathbf{1}})^2 + \cdots \nonumber \\&+ \frac{(-1)^{m+1}}{m} (\mathrm{\mathbf{X}}-\mathrm{\mathbf{1}})^m . \end{aligned}$$
(12)

Figure 1 illustrates the performance of the diagonal Padé approximants (\(m=n\)) and of the truncated Taylor series in approximating the scalar logarithm function \(\log (x)\). It can be seen that the Padé approximants behave well in a very wide range of values of x (\(0.01< x< 100\)), and even the Padé approximant of low-order (2, 2) provides a good approximation of the logarithm. On the contrary, the truncated Taylor series rapidly diverge for \(x>2\) regardless of the truncation order.

Fig. 1
figure 1

Approximation of the scalar logarithm function \(\log (x)\) at \(x=1\) by a Padé approximants and b truncated Taylor series of selected orders

Full size image

To assess the approximation error quantitatively, the absolute relative error of selected diagonal Padé approximants and truncated Taylor series is shown in Fig. 2. It follows that for \(0.2< x < 5\), the relative error for the Padé approximant of order (2, 2) lies mostly below \(10^{-2}\) and that for order (5, 5) below \(10^{-5}\). For x close to unity, the accuracy of Padé approximants is remarkably high, for instance, the Padé approximant of order (5, 5) is exact to the machine precision for \(0.9< x < 1.1\). On the other hand, the relative error for the truncated Taylor series is significantly higher than that for the Padé approximants. It can be seen that, as the order of approximation increases, the accuracy of the truncated Taylor series does not grow as rapidly as in the case of the Padé approximants.

Fig. 2
figure 2

Relative error of the approximation of the scalar logarithm function \(\log (x)\) by a Padé approximants and b truncated Taylor series of selected orders. The dashed line in panel b shows the relative error of the Padé approximant of order (3, 3), as in panel a. The vertical line at \(x=2\) in panel b indicates the upper limit of convergence of the Taylor series

Full size image

As shown by Kenney and Laub [25], under the condition that \(\Vert \mathrm{\mathbf{X}}\Vert \le c\), where \(\Vert \cdot \Vert \) denotes a matrix norm, the error in the approximation of the matrix logarithm, \(\log (\mathrm{\mathbf{X}})\), by a Padé approximant of order (m, n) is bounded by the error in the approximation (of the same order) of the scalar logarithm function \(\log (c)\). Accordingly, in view of the excellent performance of Padé approximants and their superiority over truncated Taylor series in the approximation of the scalar logarithm, see Figs. 1 and 2, Padé approximants can be used for approximation of the matrix logarithm with a higher degree of confidence.

3 Hencky-type hyperelasticity

Let us first recall the basic notions of the finite-strain framework. The primary kinematic quantity is the deformation gradient \(\mathrm{\mathbf{F}}\), defined as the reference gradient of the deformation mapping. The deformation gradient admits the polar decomposition, \(\mathrm{\mathbf{F}}=\mathrm{\mathbf{R}}\mathrm{\mathbf{U}}\), where \(\mathrm{\mathbf{U}}\) is the symmetric stretch tensor and \(\mathrm{\mathbf{R}}\) is a rotation tensor. The right Cauchy–Green deformation tensor \(\mathrm{\mathbf{C}}\),

$$\begin{aligned} \mathrm{\mathbf{C}} = \mathrm{\mathbf{F}}^\mathrm{T}\mathrm{\mathbf{F}} = \mathrm{\mathbf{U}}^2 , \end{aligned}$$
(13)

is used to define the Green strain tensor \(\mathrm{\mathbf{E}}\),

$$\begin{aligned} \mathrm{\mathbf{E}} = \mathrm{\mathbf{E}}^{(2)} = \frac{1}{2} (\mathrm{\mathbf{C}}-\mathrm{\mathbf{1}}) , \end{aligned}$$
(14)

and the Hencky (or logarithmic) strain tensor \(\mathrm{\mathbf{H}}\),

$$\begin{aligned} \mathrm{\mathbf{H}} = \mathrm{\mathbf{E}}^{(0)} = \log (\mathrm{\mathbf{U}}) = \frac{1}{2}\log (\mathrm{\mathbf{C}}) , \end{aligned}$$
(15)

both belonging to the Seth–Hill family [26]. The well-known benefit of the logarithmic strain measure is that its trace, \({\text {tr}}\mathrm{\mathbf{H}}=\log (J)\), \(J=\det \mathrm{\mathbf{F}}\), describes the volumetric deformation, and the corresponding split into volumetric and isochoric parts is additive, as in the small-strain framework.

A hyperelastic material model is specified by the elastic strain energy \(W=W(\mathrm{\mathbf{F}})\) such that

$$\begin{aligned} \mathrm{\mathbf{P}} = \frac{\partial W}{\partial \mathrm{\mathbf{F}}} , \end{aligned}$$
(16)

where \(\mathrm{\mathbf{P}}\) is the first Piola-Kirchhoff stress tensor. The standard transformation rules relate \(\mathrm{\mathbf{P}}\) to the Cauchy stress, \({\varvec{\sigma }}=J^{-1}\mathrm{\mathbf{P}}\mathrm{\mathbf{F}}^\mathrm{T}\), and to the second Piola-Kirchhoff stress tensor \(\mathrm{\mathbf{S}}=\mathrm{\mathbf{F}}^{-1}\mathrm{\mathbf{P}}\).

In the St. Venant-Kirchhoff model, the elastic strain energy is quadratic in the Green strain tensor \(\mathrm{\mathbf{E}}\),

$$\begin{aligned} W_\mathrm{SVK} = \frac{1}{2} \mathrm{\mathbf{E}} \cdot \mathbb {L} \mathrm{\mathbf{E}} , \quad W_\mathrm{SVK}^\mathrm{(iso)} = \frac{\kappa }{2} ( {\text {tr}} \mathrm{\mathbf{E}} )^2 + \mu \, \mathrm{\mathbf{E}}' \cdot \mathrm{\mathbf{E}}', \end{aligned}$$
(17)

where \(W_\mathrm{SVK}^\mathrm{(iso)}\) corresponds to the special case of isotropy. Here, \(\mathbb {L}\) is the usual fourth-order elastic stiffness tensor, cf. [27], \(\kappa \) is the bulk modulus, \(\mu \) the shear modulus, and \(\mathrm{\mathbf{E}}'\) the deviatoric part of \(\mathrm{\mathbf{E}}\), \(\mathrm{\mathbf{E}}'=\mathrm{\mathbf{E}}-\frac{1}{3}({\text {tr}}\mathrm{\mathbf{E}})\mathrm{\mathbf{1}}\). The St. Venant-Kirchhoff model is the simplest generalization of the linear small-strain elasticity to the finite-deformation framework. At the same time, as it is well known, the St. Venant-Kirchhoff strain energy lacks important properties, notably rank-one convexity and polyconvexity [28, 29].

The Hencky elastic strain energy is quadratic in the Hencky strain \(\mathrm{\mathbf{H}}\),

$$\begin{aligned} W_\mathrm{H} = \frac{1}{2} \mathrm{\mathbf{H}} \cdot \mathbb {L} \mathrm{\mathbf{H}} , \quad W_\mathrm{H}^\mathrm{(iso)} = \frac{\kappa }{2} ( {\text {tr}} \mathrm{\mathbf{H}} )^2 + \mu \, \mathrm{\mathbf{H}}' \cdot \mathrm{\mathbf{H}}' , \end{aligned}$$
(18)

and, as above, \(W_\mathrm{H}^\mathrm{(iso)}\) corresponds to the isotropic case. The Hencky strain energy is not rank-one convex everywhere [8, 30]. However, compared to the St. Venant-Kirchhoff strain energy, it possesses a much larger domain of rank-one convexity. Therefore, the range of strains in which the model behaves well is much wider [7, 15]. For instance, it has been shown that, for positive Lamé constants in the isotropic case, the respective elastic strain energy is elliptic whenever every principal stretch is between 0.212 and 1.396 [31]. The Hencky elastic strain energy is thus favoured over the St. Venant-Kirchhoff energy when the analyses are not restricted to small elastic strains, e.g. [13].

Recently, the exponentiated Hencky model has been introduced [8] as a modification of the isotropic Hencky model in the following form,

$$\begin{aligned} W_\mathrm{eH}^\mathrm{(iso)} = \frac{\kappa }{2\hat{k}} \exp (\hat{k}( {\text {tr}} \mathrm{\mathbf{H}} )^2) + \frac{\mu }{k} \exp ( k \mathrm{\mathbf{H}}' \cdot \mathrm{\mathbf{H}}' ) , \end{aligned}$$
(19)

where k and \(\hat{k}\) are two additional dimensionless material parameters. This energy is polyconvex in the two-dimensional case for positive shear and bulk moduli and when \(k>1/3\) and \(\hat{k}>1/8\) [32], and has a larger domain of rank-one convexity with respect to the Hencky elastic strain in the three-dimensional case. The exponentiated Hencky model has also been extended to the anisotropic cases of transverse isotropy and orthotropy [33].

The common drawback associated with practical application of the Hencky-type models is that the tensor logarithm must be computed along with its first and second derivatives that are needed to determine the stress and the tangent operator, respectively. Here, it is appealing to employ the Padé approximants that have a simple, explicit form and deliver a good approximation of the tensor logarithm, as illustrated in Sect. 2.2. For instance, application of the Padé approximant of order (2, 2) leads to the following approximation of the Hencky strain, cf. Eq. (10),

$$\begin{aligned} \mathrm{\mathbf{H}} \approx \mathrm{\mathbf{H}}_{(2,2)} = \frac{3}{2} (\mathrm{\mathbf{C}}^2 - \mathrm{\mathbf{1}}) (\mathrm{\mathbf{C}}^2 + 4\mathrm{\mathbf{C}} + \mathrm{\mathbf{1}})^{-1} , \end{aligned}$$
(20)

where \(\mathrm{\mathbf{H}}_{(m,n)}\) is introduced as a notation for the approximate Hencky strain based on the Padé approximant of order (m, n).

Importantly, since \(\mathrm{\mathbf{H}}_{(m,n)}\) is given as an explicit formula, its derivatives can be derived directly, also by employing general-purpose automatic differentiation (AD) tools.Footnote 1 Note that a closed-form representation of the matrix logarithm has been recently developed by Hudobivnik and Korelc [11], which relies on the AceGen system [35] and its implementation of the AD technique. However, the closed-form matrix logarithm is not directly available in other computing environments, thus the Padé approximants may be an attractive alternative, particularly when a small approximation error is not a major drawback. For instance, in Hencky-type hyperelasticity, the specific form of the elastic strain energy is anyway a model of a more complex reality, so that a reasonable approximation of the logarithmic strain seems acceptable.

In Sect. 5, performance of the Padé approximants of the matrix logarithm is examined in sample finite-element computations. The computations are carried out using the AceFEM system, which is integrated with the AceGen system, hence the closed-form representation of the matrix logarithm, available in AceGen, is used as a reference.

4 Remark on the stress conjugate to the logarithmic strain

The logarithmic strain \(\mathrm{\mathbf{H}}\) has definite advantages and is frequently used in various contexts. However, its work-conjugate stress, denote it by \(\mathrm{\mathbf{T}}^{(0)}\), is not available in closed form. In the case of isotropic elasticity, \(\mathrm{\mathbf{T}}^{(0)}\) is simply the rotated Kirchhoff stress \({\varvec{\tau }}=J{\varvec{\sigma }}\), namely \(\mathrm{\mathbf{T}}^{(0)}=\mathrm{\mathbf{R}}^\mathrm{T}{\varvec{\tau }}\mathrm{\mathbf{R}}\). In the general case, the form of \(\mathrm{\mathbf{T}}^{(0)}\) depends on the multiplicity of the principal stretches [36]. In view of the ill-posedness of the eigenproblem in the vicinity of repeated eigenvalues, the mentioned closed-form expressions may be inconvenient in practice. The stress conjugate to the Eulerian logarithmic strain \(\log (\mathrm{\mathbf{V}})\) is discussed in [37].

The stress power can be expressed in terms of any pair of conjugate strain and stress measures. Consider thus the logarithmic strain \(\mathrm{\mathbf{H}}\) with its conjugate stress \(\mathrm{\mathbf{T}}^{(0)}\) and, for instance, the Green strain \(\mathrm{\mathbf{E}}=\mathrm{\mathbf{E}}^{(2)}\) with its conjugate second Piola-Kirchhoff stress \(\mathrm{\mathbf{S}}=\mathrm{\mathbf{T}}^{(2)}\). The stress power then reads

$$\begin{aligned} \mathrm{\mathbf{S}} \cdot \dot{\mathrm{\mathbf{E}}} = \mathrm{\mathbf{T}}^{(0)} \cdot \dot{\mathrm{\mathbf{H}}} = \mathrm{\mathbf{T}}^{(0)} \cdot \mathbb {H} \dot{\mathrm{\mathbf{E}}} = (\mathbb {H}^\mathrm{T} \mathrm{\mathbf{T}}^{(0)}) \cdot \dot{\mathrm{\mathbf{E}}} , \end{aligned}$$
(21)

where \(\mathbb {H}=\partial \mathrm{\mathbf{H}}/\partial \mathrm{\mathbf{E}}\) is a fourth-order tensor such that \(\dot{\mathrm{\mathbf{H}}}=\mathbb {H}\dot{\mathrm{\mathbf{E}}}\). Since \(\dot{\mathrm{\mathbf{E}}}\) is arbitrary, we have

$$\begin{aligned} \mathrm{\mathbf{S}} = \mathbb {H}^\mathrm{T} \mathrm{\mathbf{T}}^{(0)} , \quad \mathrm{\mathbf{T}}^{(0)} = \mathbb {H}^\mathrm{-T} \mathrm{\mathbf{S}} . \end{aligned}$$
(22)

Note that tensors \(\mathrm{\mathbf{H}}\) and \(\mathrm{\mathbf{E}}\) are uniquely defined one in terms of the other, i.e. \(\mathrm{\mathbf{H}}=\frac{1}{2}\log (2\mathrm{\mathbf{E}}+\mathrm{\mathbf{1}})\) and \(\mathrm{\mathbf{E}}=\frac{1}{2}(\exp (2\mathrm{\mathbf{H}})-\mathrm{\mathbf{1}})\). Hence, if \(\dot{\mathrm{\mathbf{H}}}=\mathrm{\mathbf{0}}\) then \(\dot{\mathrm{\mathbf{E}}}=\mathrm{\mathbf{0}}\) is the only solution of the linear equation \(\mathbb {H}\dot{\mathrm{\mathbf{E}}}=\mathrm{\mathbf{0}}\), and thus \(\mathbb {H}\) is invertible. Tensor \(\mathbb {H}\) possesses the minor symmetries, \((\mathbb {H})_{ijkl}=(\mathbb {H})_{jikl}=(\mathbb {H})_{ijlk}\), but not the major symmetry, \((\mathbb {H})_{ijkl}\ne (\mathbb {H})_{klij}\), in general.

It follows that the closed-form representation of the matrix logarithm, as developed by Hudobivnik and Korelc [11], allows one to evaluate \(\mathbb {H}\) as the derivative of \(\mathrm{\mathbf{H}}\) with respect to \(\mathrm{\mathbf{E}}\) (or another suitable strain measure) and thus to compute a numerically-exact representation of \(\mathrm{\mathbf{T}}^{(0)}\). Note that the closed-form representation of the matrix logarithm [11] delivers a computer code (algorithm) that is generated using the AD technique. The explicit formulae for \(\mathbb {H}\) and \(\mathrm{\mathbf{T}}^{(0)}\) would thus not be available, rather, the respective computer code would be derived using the AD technique.

Assume now that the Hencky strain \(\mathrm{\mathbf{H}}\) is approximated by \(\mathrm{\mathbf{H}}_{(m,n)}\), \(\mathrm{\mathbf{H}}\approx \mathrm{\mathbf{H}}_{(m,n)}\), where (m, n) is the order of the corresponding Padé approximant of the tensor logarithm, see, for instance, \(\mathrm{\mathbf{H}}_{(2,2)}\) given by Eq. (20). Note that \(\mathrm{\mathbf{H}}_{(m,n)}\) itself is a strain measure, and thus one can define the corresponding conjugate stress \(\mathrm{\mathbf{T}}_{(m,n)}^{(0)}\). Starting from the stress power \(\mathrm{\mathbf{T}}_{(m,n)}^{(0)}\cdot \dot{\mathrm{\mathbf{H}}}_{(m,n)}\), in analogy to Eq. (22), we have

$$\begin{aligned} \mathrm{\mathbf{T}}^{(0)} \approx \mathrm{\mathbf{T}}_{(m,n)}^{(0)} = \mathbb {H}_{(m,n)}^\mathrm{-T} \mathrm{\mathbf{S}} , \end{aligned}$$
(23)

where \(\mathbb {H}_{(m,n)}=\partial \mathrm{\mathbf{H}}_{(m,n)}/\partial \mathrm{\mathbf{E}}\). Since \(\mathrm{\mathbf{H}}_{(m,n)}\) is an approximation of \(\mathrm{\mathbf{H}}\), \(\mathrm{\mathbf{T}}_{(m,n)}^{(0)}\) provides an approximation of \(\mathrm{\mathbf{T}}^{(0)}\). The accuracy of approximation is controlled by the approximant order (m, n) and can be very high for moderate strains, as shown in Sect. 2.2. Note that \(\mathrm{\mathbf{H}}_{(m,n)}\) is given by a single explicit formula, see e.g. Eq. (20), and thus \(\mathbb {H}_{(m,n)}\) can also be derived in an explicit form. Alternatively, a general-purpose AD tool can be used to generate the respective computer code.

To summarize, we have shown that \(\mathrm{\mathbf{T}}^{(0)}\), the stress conjugate to the logarithmic strain \(\mathrm{\mathbf{H}}=\mathrm{\mathbf{E}}^{(0)}\), can be computed, i.e. evaluated numerically, using the AD technique and the related closed-form representation of the matrix logarithm [11]. On the other hand, \(\mathrm{\mathbf{T}}^{(0)}\) can be approximated by an explicit formula, with controlled accuracy, using the Padé approximants of the tensor logarithm.

5 Performance of Padé approximants in finite-element computations

5.1 Preliminaries

In this section, a finite-element study is performed with the aim to examine the applicability and computational efficiency of Padé approximants in Hencky-type hyperelasticity. Finite-element simulations are performed for Hencky and exponentiated Hencky models, Eqs. (18) and (19), respectively, employing Padé approximants of the Hencky strain \(\mathrm{\mathbf{H}}\) and the results are compared to those obtained for the respective Hencky-type models employing truncated Taylor series. Hereinafter, to avoid confusion, the Hencky models in Eq. (18) will be referred to as ‘classical Hencky models’. Diagonal Padé approximants of order (2, 2) through (5, 5) and Taylor series truncated to order 2 through 8 are studied. At the same time, to provide a baseline for this study, simulations are also performed for Hencky-type models with the exact closed-form representation of the matrix logarithm [11]. In addition, some comparisons are made with the results obtained for the St. Venant-Kirchhoff model, Eq. (17), and for the neo-Hookean model.Footnote 2

Indentation of a hyperelastic block is examined as a model problem and two types of materials are considered, namely isotropic and anisotropic of cubic symmetry. In the latter case, simulations are carried out for [001] and [011] orientations and the material is assumed to possess the elastic constants of single-crystal austenitic CuAlNi, which is characterized by high anisotropy (\(c_{11}=142\) GPa, \(c_{44}=96\) GPa, \(c_{12}=126\) GPa [39], with the Zener ratio \(A=2c_{44}/(c_{11}-c_{12})=12\), for the interpretation of the elastic constants \(c_{ij}\), see [27]). Note that we use the elastic constants of CuAlNi for illustration purposes only, and the range of elastic strains considered in the present finite-element study largely exceeds that of metals. Elastic properties of the isotropic material (\(\mu =39\) GPa and \(\kappa =128\) GPa, so that \(E=106\) GPa and \(\nu =0.36\)) are determined via averaging the anisotropic elastic constants of cubic CuAlNi using the Voigt–Reuss–Hill averaging scheme [40]. The isotropic exponentiated Hencky model, Eq. (19), involves two dimensionless material parameters in addition to the bulk and shear moduli. Following [41], \(k=2\) and \(\hat{k}=3\) are adopted here (the effect of parameters k and \(\hat{k}\) on the material response is discussed in [8]).

5.2 Finite-element implementation

In the present study, the external load (indentation) is applied through a frictionless contact interaction between the hyperelastic body and a rigid spherical indenter. The contact problem is solved by enforcing the unilateral contact constraint using the augmented Lagrangian technique [42]. The related formulation is quite standard and thus is not provided here.

The finite-element implementation involves the derivation of the weak form of equilibrium (virtual work principle) followed by the discretization of the weak form using the standard Galerkin method. The resulting discretized governing equations are then solved by using the Newton method simultaneously for the nodal displacements and contact Lagrange multipliers. In the present implementation, the standard isoparametric 8-noded hexahedral element with trilinear interpolation functions is used for the displacement field.

The finite-element code generation is efficiently performed via the symbolic programming tool, AceGen, which exploits the automatic differentiation (AD) and expression optimization techniques [35, 43]. As discussed in Sect. 3, AceGen provides a closed-form representation of the matrix logarithm [11] that allows a straightforward and exact implementation of the Hencky-type hyperelasticity models. In fact, models based on this closed-form representation play an important role in the present study, as they serve as a reference and allow a clear evaluation of the performance of Padé approximants in hyperelasticity problems.

The finite-element computations are performed in AceFEM, a finite-element environment that is fully interfaced with AceGen. A direct linear solver (MKL PARDISO) has been used.

5.3 Indentation of a hyperelastic block

An elastic block of the size \(L \times L \times H =40 \times 40 \times 20\) mm\(^3\) is compressed by a rigid spherical indenter of the radius \(R=15\) mm, Fig. 3. The vertical displacements are constrained at the bottom surface, otherwise the block is free to expand laterally during indentation. By exploiting the symmetry of the problem, only one quarter of the block is modelled, with adequate boundary conditions imposed along the planes of symmetry. The actual computational domain is discretized in \(20 \times 20 \times 20\) hexahedral elements leading to 8000 elements and 26,901 degrees of freedom, comprising of displacements and contact Lagrange multipliers.

Fig. 3
figure 3

Indentation of a hyperelastic block: geometry and finite-element mesh

Full size image
Fig. 4
figure 4

Normalized load–normalized indentation depth responses (\(P/(E R^2)\) vs. h/R) obtained for a isotropic models, b cubic anisotropic models with [001] orientation, and c cubic anisotropic models with [011] orientation. The markers represent the last converged load steps. The curly arrow in panel c indicates the load step at which folding instability occurs

Full size image

First, a brief discussion of the simulation results for different hyperelasticity models is provided, and whenever applicable, the closed-form representation of the matrix logarithm is used. The normalized load–normalized indentation depth responses (for convenience, referred to as ‘P–h response’ in the sequel) obtained for isotropic and cubic anisotropic models are shown in Fig. 4. To make the quantities dimensionless, the indentation depth is normalized by the radius R, and the load is normalized by \(E R^2\), where the Young’s modulus E is determined as described in Sect. 5.1. In each case, the simulation, which is carried out using an adaptive load-step control, is continued up to the failure of the Newton scheme and the last converged solution is indicated by a marker at the end of the P–h curve.

It can be seen in Fig. 4 that the P–h responses have the same slope close to zero indentation depth. This observation stems from the fact that all isotropic and all cubic anisotropic models coincide in the small strain regime. In both isotropic and cubic anisotropic cases, the St. Venant-Kirchhoff model exhibits a considerably poorer performance with respect to the other models, as the corresponding Newton schemes crash at early stages of indentation. As discussed in Sect. 3, this is probably due to the loss of rank-one convexity that implies the poor performance of the St. Venant-Kirchhoff elastic strain energy under compression [15, 29, 44]. Note also that, even for small indentation depths, the responses yielded by the St. Venant-Kirchhoff models deviate significantly from the other responses, see the enlarged views in Fig. 4. Among isotropic models, the neo-Hookean model has achieved the largest indentation depth of \(h/R=1.28\). The indentation depth of \(h/R=1.1\) has been achieved by the exponentiated Hencky model. However, in view of the strong stiffening characteristic for this model at large elastic strains, the corresponding maximum normalized load is far greater than in the other models.

For the purpose of illustration, the deformed configurations obtained for the Hencky-type models at the final load step of the respective simulations are shown in Fig. 5. The corresponding maximum and minimum principal stretches achieved within the hyperelastic block are presented in Fig. 6. These are provided here as a reference for the study reported below, in which the Padé approximants of the Hencky strain will be employed, and the accuracy of the approximation of the tensor logarithm depends on the eigenvalues of its argument, cf. Sect. 2. For instance, the maximum and minimum principal stretches of, respectively, 2.89 and 0.16 have been obtained for the isotropic classical Hencky model. Accordingly, with reference to Figs. 1 and 2, a general idea regarding the corresponding maximum error of the approximation of the Hencky strain can be gained.

Fig. 5
figure 5

The deformed configuration for the isotropic and cubic anisotropic Hencky-type models at the last converged load step of the respective computations (top). For [011] cubic classical Hencky model, the deformed configuration at the load step just before the occurrence of the folding instability is also shown (bottom), see the curly arrow in Fig. 4c. The color contours indicate the vertical displacement

Full size image

The final deformed configuration obtained for the classical Hencky model with [011] cubic anisotropy reveals an oscillatory deformation mode, which is characterized by formation of a folding pattern in the mesh layers beneath the indenter within the y–z plane. The curly arrow in Fig. 4c indicates the load step at which the folding instability is encountered. Such instability has also been observed when the improved enhanced-strain element TSCG12 [45] has been used instead of the standard isoparametric element, which suggests that the instability is related to the material model itself (possibly due to the loss of rank-one convexity of the strain energy) rather than to the finite-element formulation. The corresponding deformed configuration at the load step just before the occurrence of the instability is depicted in Fig. 5. A similar deformation pattern was also noticed for the isotropic St. Venant-Kirchhoff model under compression [44], where such behaviour was attributed to the lack of compressive resistance of the model as a consequence of exceeding the domain of rank-one convexity. Although the folding instability affects the reliability of the simulation results, the respective case has not been excluded from the study of the approximation methods, which is presented below.

Fig. 6
figure 6

Maximum and minimum principal stretches achieved for the Hencky-type models at the last converged load step, see Fig. 5

Full size image

Now, the simulation results are discussed regarding the performance of the Padé and Taylor-series approximations in Hencky-type hyperelasticity. Figs. 7, 8, 9 and 10 illustrate the P–h responses obtained for the Hencky-type models with different representation of the Hencky strain \(\mathrm{\mathbf{H}}\). Here, we consider the isotropic classical Hencky model (Fig. 7), the isotropic exponentiated Hencky model (Fig. 8), and the anisotropic classical Hencky model with [001] and [011] orientations (Figs. 9 and 10, respectively). It can be seen that in both isotropic and cubic anisotropic cases, the models with Padé approximants outperform those with truncated Taylor series in terms of the maximum attainable indentation depth, and the corresponding P–h responses are in a good agreement with that obtained for the respective models with the closed-form representation. The superiority of the models with Padé approximants is pronounced even for the Padé approximant of low order (2, 2). Interestingly, models with odd-order Taylor series show a better performance with respect to those with even-order, note the difference in the respective approximations of the scalar logarithm function in Fig. 1b. For instance, in isotropic cases, the former have reasonably approximated the P–h responses in a wide range of indentation depth, i.e. for \(h/R < 0.75\), but have visibly overpredicted the response outside this range. The latter, however, exhibit a narrower range of functionality, e.g. \(h/R < 0.6\) for the classical Hencky model with the Taylor series of order 8, and terminate at earlier stages of the simulation due to the failure of the Newton scheme.

Fig. 7
figure 7

Performance of isotropic classical Hencky model with different representations of the Hencky strain \(\mathrm{\mathbf{H}}\). The P–h responses for the models with a Padé approximants, b even-order Taylor series, and c odd-order Taylor series. The markers indicate the last converged load steps

Full size image
Fig. 8
figure 8

Performance of isotropic exponentiated Hencky model with different representations of the Hencky strain \(\mathrm{\mathbf{H}}\). The P–h responses for the models with a Padé approximants, b even-order Taylor series, and c odd-order Taylor series. The markers indicate the last converged load steps

Full size image
Fig. 9
figure 9

Performance of [001]-oriented cubic anisotropic classical Hencky model with different representations of the Hencky strain \(\mathrm{\mathbf{H}}\). The P–h responses for the models with a Padé approximants, b even-order Taylor series, and c odd-order Taylor series. The markers indicate the last converged load steps

Full size image
Fig. 10
figure 10

Performance of [011]-oriented cubic anisotropic classical Hencky model with different representations of the Hencky strain \(\mathrm{\mathbf{H}}\). The P–h responses for the models with a Padé approximants, b even-order Taylor series, and c odd-order Taylor series. The markers and curly arrows indicate, respectively, the last converged load steps and the load steps at which the instability occurs

Full size image

It is noteworthy to point out that, despite the occurrence of the folding instability in the [011] cubic anisotropic case, the models with Padé approximants have yielded a reasonable approximation of the corresponding P–h response. Moreover, all the Padé approximant models have correctly captured the load at which the instability is encountered. On the other hand, among different Taylor-series models, only the Taylor-series models of orders 5 and 7 have been able to capture the instability, however, with an inconsistent load level with respect to that of the closed-form representation, see the curly arrows in Fig. 10a,c.

Next, a detailed study is carried out with the aim to investigate the computational efficiency of the Hencky-type models with Padé and Taylor-series approximations. The results reported in Figs. 7, 8, 9 and 10 indicate an overall poor performance of the low-order Taylor-series models. Moreover, a relatively narrow range of functionality was observed for the even-order Taylor-series models. As a consequence, only the Taylor-series model of order 7 is included in the further analyses.

Two representative Hencky-type models are considered in this study, namely the isotropic and [001]-oriented cubic anisotropic classical Hencky models. Recall that a folding instability appeared in the case of [011]-orientation, thus the corresponding model has been excluded from this study. A twice finer mesh has been utilized, which consists of 64, 000 elements and leads to approximately 203, 400 degrees of freedom. For the sake of illustration, the effect of mesh density on the final deformed configuration and on the P–h response for the [001]-oriented Hencky model is shown in Fig. 11.

Fig. 11
figure 11

The effect of mesh density on a the deformed configuration for [001]-oriented cubic anisotropic classical Hencky model at the last converged load step of computation, and b the P–h response. The markers in panel b indicate the last converged load steps

Full size image

To evaluate the computational efficiency of different approximation methods, the total number of load steps required to complete the simulations and the average assembly time per iteration are compared in the bar charts in Fig. 12. The average assembly time per iteration is defined as the average time required to assemble the global residual vector and the global tangent matrix and is influenced by the computational cost of evaluation of the Hencky strain and its first and second derivatives. It has been observed that the total simulation wall-clock time follows exactly the same trend as that of the total number of load steps. To make the comparisons relevant, all the computations are performed up to the normalized indentation depth \(h_\text {max}/R=0.75\). Recall that an adaptive load-stepping procedure is employed in the simulations, where the size of the load step is automatically tuned (increased or decreased) based on the number of global Newton iterations needed to converge to the solution at the previous load step.

Fig. 12
figure 12

Computational performance of the classical Hencky models with different representations of the Hencky strain \(\mathrm{\mathbf{H}}\) in terms of the average assembly time per iteration and the total number of load steps: a isotropic classical Hencky and b [001]-oriented cubic classical Hencky

Full size image

It can be seen from Fig. 12 that, as expected, the average assembly time increases with increasing the order of Padé approximant, but only to a small extent, such that in both isotropic and cubic anisotropic cases, the assembly time obtained for the model with the Padé approximant of order (5 ,5) is only about 1.25 times greater than that for the model with order (2, 2). The computational cost of the model with the Taylor series of order 7 is similar to that with the Padé approximant of order (4, 4). However, as illustrated in Figs. 7, 8, 9 and 10, the former performs visibly worse in terms of overall accuracy. Comparison with the assembly time for the model with the closed-form representation shows a high computational efficiency of the implementation based on the Padé approximants. Recall that code generation is here performed using AceGen, which combines automatic differentiation and expression optimization techniques. The conclusions concerning the computational efficiency are thus specific for the present implementation approach.

Concerning the total number of load steps, in the isotropic case, the model with the Taylor series of order 7 exhibits the largest value, thus the longest simulation wall-clock time. At the same time, the total number of load steps of the models with Padé approximants is not significantly influenced by the order of approximation. In the cubic case, the best performance in terms of the number of load steps has been achieved by the model with Padé approximant of order (2, 2).

6 Conclusion

We have employed the Padé approximants of the tensor logarithm as a simple, accurate and computationally efficient method to approximate the logarithmic (Hencky) strain tensor. For moderate strains, accuracy close to the machine precision can be obtained depending on the approximation order. The resulting approximate Hencky strain tensor is obtained in an explicit form, which admits a straightforward computer implementation, including evaluation of the first and second derivatives that are needed to compute the stress and the tangent operator. In this context, we have also discussed the possibility of computing, in an exact or approximate manner, the stress tensor conjugate to the Hencky strain.

As an application, we have considered Hencky-type hyperelastic models in which the elastic strain energy is expressed in terms of the Hencky strain, in particular, the Hencky strain energy that is quadratic in the Hencky strain. Practical implementation of the Hencky-type models is significantly simplified when the Hencky strain is approximated, e.g., using a Padé approximant. At the same time, the corresponding error introduced to the elastic strain energy function is justified considering that the specific elastic strain energy is anyway a model, i.e. an approximation of a more complex phenomenon.

The results of the finite-element study illustrate a superior performance of the Hencky-type models employing the Padé approximants with respect to those employing the truncated Taylor series in both isotropic and cubic anisotropic hyperelasticity problems. More specifically, the superiority is reflected in a larger attainable strain range, in a better agreement of the mechanical response with that obtained for the model with numerically-exact closed-form representation of the matrix logarithm, and in a competitive computational efficiency. The results also demonstrate that the models with odd-order Taylor series possess a wider range of functionality compared to those with even order. A poor performance of the St. Venant-Kirchhoff model has also been illustrated in both isotropic and anisotropic hyperelasticity.