Crossing Scales: Data-Driven Determination of the Micro-scale Behavior of Polymers From Non-homogeneous Tests at the Continuum-Scale

We propose an efficient method to determine the micro-structural entropic behavior of polymer chains directly from a sufficiently rich non-homogeneous experiment at the continuum scale. The procedure is developed in 2 stages: First, a Macro-Micro-Macro approach; second, a finite element method. Thus, we no longer require the typical stress-strain curves from standard homogeneous tests, but we use instead the applied/reaction forces and the displacement field obtained, for example, from Digital Image Correlation. The approach is based on the P-spline local approximation of the constituents behavior at the micro-scale (a priori unknown). The sought spline vertices determining the polymer behavior are first pushed up from the micro-scale to the integration point of the finite element, and then from the integration point to the element forces. The polymer chain behavior is then obtained immediately by solving a linear system of equations which results from a least squares minimization error, resulting in an inverse problem which crosses material scales. The result is physically interpretable and directly linked to the micro-structure of the material, and the resulting polymer behavior may be employed in any other finite element simulation. We give some demonstrative examples (academic and from actual polymers) in which we demonstrate that we are capable of recovering “unknown” analytical models and spline-based constitutive behavior previously obtained from homogeneous tests.

1 Introduction

Modern applications and the easiness of 3D printing of polymers even at the micro-scale (e.g., via dual-photon polymerization), have renewed the interest in large deformation modeling of these entropic materials. Polymeric materials can now be found in a wide range of biomedical applications (stents, sutures, spinal cages, soft tissue implants, and tissue engineering scaffolds, … ), see Bergström and Hayman (2016). Even most human soft biological tissues, which are made of a matrix (elastine, proteoglycans) plus fibers (collagen), withstand large reversible deformations within the physiological range, and therefore, use hyperelasticity as ground for more complex aspects; Chagnon et al. (2015), Chagnon et al. (2017). As the simplest procedure to guarantee true elasticity (reversible, non-dissipative processes) the cornerstone in hyperelasticity is the free energy function, the state function, from which the stresses are uniquely derived from the strains (or vice-versa), regardless of path. Since the 3D strain energy function cannot be measured directly, the classical approach in constitutive modeling establishes a predefined form for the free energy. This function typically contains some parameters that are adjusted according to the experimentally (stress-strain) observed material behavior. Although it is relatively simple to tune model parameters to predict (up to a desirable precision) a single experimental curve, determining the parameters that produce accurate results for different modes of deformation is not trivial, as it is apparent from the unaccountable number of hyperelastic models available; Volokh (2016). Theoretically, if the proposed model is correct, this set of parameters should exist and although determined for specific tests they should predict well other modes of deformation. In practice, when parameters are obtained from a single experimental curve, they fail to generalize to other deformation states; this is the reason why in practice multiple tests are recommended to determine the parameters of the free energy function (Marckmann and Verron, 2006, 3, p.12). Using multiple tests to calibrate the parameters alleviates the deviations of the model for other modes of deformation (at the cost of accuracy for a given test), but at the same time it raises the question of whether the proposed form for the free energy really captures the physical phenomena behind experimental data or this assumed form is just a complex interpolation scheme that adapts its parameters to fit the curves used during calibration. We remark that if the physics behind were accurately represented, a single curve should be sufficient to capture the general multiaxial behavior of isotropic, incompressible polymers under reversible deformations.

This kind of problems has encouraged many researchers to pursue different approaches. One of them is the model-free data-driven computing paradigm. In this approach, basic conservation laws and essential constraints are satisfied but the constitutive laws are eliminated in the benefit of data; Kirchdoerfer and Ortiz (2016), Kirchdoerfer and Ortiz (2018), Eggersmann et al. (2019), Ibañez et al. (2017), Ibañez et al. (2018). Regarding the leading role of data for some of these references, works that address the efficient handling of data have also been published, see Zheng et al. (2020), Korzeniowski and Weinberg (2021). On the other hand, other approaches attempt to surrogate the constitutive law with an input-output relation through Artificial neural networks (ANNs), Nguyen-Thanh et al. (2020), Liu et al. (2020). Both approaches (model-free data-driven and surrogate-like ones) show promising results, however, the predictive capability of models that just rely on data is strongly dependent on the amount and quality of data being employed. In addition, since there is no expression for the free energy function most of those models are very difficult to interpret from a physical standpoint. It seems clear that an approach solely based on data might not be the best option for this kind of problems (the more the model needs to learn, the more data is required). This need for introducing physics information in full data-driven models has led to other works based on Physics-informed neural networks (PINNs) Liu et al. (2020) and on thermodynamically consistent data-driven approaches; González et al. (2018). Still, the increased generality of those approaches increase the amount of data required when compared to classical constitutive modeling techniques and their interpretability is much less direct. To summarize, an optimal approach for the constitutive modeling of polymers should: 1) include information about the physical equations without assuming a fixed given form for the free energy function; 2) use data to complement what we know about the physical phenomena and fill in the gaps in our knowledge, but without using more data than actually needed; 3) interpretability is also very important because understanding the solution and being able to identify its physical meaning avoids many pitfalls allowing us to search for the answer within a smaller solution space and identify spurious solutions. Interpretability also facilitates the imposition of desired (physics-based) requirements to the sought solution, for example, smoothness, monotonic increase or decrease, isotropy, etc.

With all those requirements in mind we developed the WYPiWYG (What-You-Prescribe is What-You -Get) approach to constitutive modelling, Latorre and Montáns (2013), based on some seminal ideas from the Sussman-Bathe model for isotropic, incompressible materials, Sussman and Bathe (2009). The WYPiWYG approach determines the free energy function or its contributions, but in contrast to classical phenomenological models, which presume a form for the energy function and fit the model parameters to the experimental data, our approach starts with some basic fundamental assumptions about the material behavior (isotropy/anisotropy, Valanis-Landel decomposition, invariant-based contributions to the energy function) and then obtains numerically the constitutive equation from equilibrium using a local approximation scheme based on splines. This local approximation philosophy is similar to the way shape functions in finite elements interpolate the displacement field, instead of computing coefficients of predefined analytical functions as in the Navier and Rayleigh methods. The generality of this approach is demonstrated on the models elaborated for anisotropic materials, Latorre and Montáns (2014), auxetic materials Crespo and Montáns (2018) and models for the active and passive response of skeletal muscle, Moreno et al. (2020). So far phenomenological WYPiWYG hyperelasticity circumvents the need to prescribe the shape of the energy function while maintaining the same model interpretability that the phenomenological models have. However, the amount of data required to characterize the behavior of polymeric materials is similar to the amount of data required by other parametric phenomenological models like the Ogden Model, Ogden (1972).

Other alternative to phenomenological models are those based on the micro-structure which employ additional information about the structure of polymers to get better predictive capabilities with fewer data. Most micro-structural models assume that the polymer is fully entropic and thus all the work employed in its deformation directly translates into a variation of its entropy. This physical insight has been exploited by researches, leading to well known models as the Neo-Hookean model, Flory and Rehner (1943); Gent (1989); Treloar (1975), and the 8-chain model, Arruda and Boyce (1993). The expressions of those models depend just on the first invariant of the Green-Cauchy tensor, $I_1^C$ and some material parameters, but in contrast to the phenomenological ones, the material parameters are linked to the micro-structure resulting in some additional physical insight. Although they were conceived to be characterized from a simple extension test, the results on other modes of deformation are not satisfactory, even with the additional information about the microstructural behavior. This fall from expectations added up to the conclusions of Mooney and its Mooney plots Mooney (1940), which showed that $I_1^C$ was not the right (or at least the unique) variable to describe the polymer behavior due to the controversial slope C₂ that consistently appeared on simple extension experimental data.

The extension of the WYPiWYG approach to microstructural modelling with the aim of overcoming those dificulties resulted in a Macro-Micro-Macro (MMM) approach to obtain the polymer constituents behavior directly from experimental data with no assumptions about its analytical form or parameters to calibrate; Amores et al. (2020). Just some basic assumptions were made: 1) homogenization of the chains free energy to obtain the free energy of the continuum $\mathrm{\Psi }\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)={\int }_S{\psi }_{ch}\left({\lambda }_{ch}\right)dS/S$, and 2) the computation of the micro-stretch variable is obtained from the continuum stretch tensor, λ_ch = r ⋅U ⋅r (a non-affine measure of deformation in agreement with the lack of relevant contribution of chains orientation change in the entropy reduction). The MMM approach predicts well any general deformation mode requiring just a single experimental curve to characterize the chain behavior, see (Amores et al., 2020, Figure 3). Since a similar Data-Driven MMM framework using the affine deformation measure ${\left({\lambda }_{ch}^C\right)}^2=\mathit{r}\cdot \mathit{C}\cdot \mathit{r}$ was not able to offer the same results, in Amores et al. (2021) we questioned the affine micro stretch assumption from theoretical grounds, which seems to be the most popular in micro-structural models; Treloar (1975), Arruda and Boyce (1993), Alastrué et al. (2009), Sáez et al. (2011), Khiêm and Itskov (2016). In that same work the non-affine measure of deformation λ_ch = r ⋅U ⋅r did show to be in accordance with the “controversial” C₂ slope observed by Mooney that up to the date could not be successfully explained from the classical statistical theory.

The framework presented in Amores et al. (2020) seems to be in accordance with both the chain statistical theory and experimental results, but needs homogeneous tests to characterize the chain behavior. Hence, our work here is to pursue a more general approach by employing arbitrary continuum non-homogeneous tests and using Digital Image Correlation (DIC), crossing scales from the continuum to the polymer constituent macromolecules.

The procedure consists of linking two stages. One is the previously introduced MMM method, and the other one is to link that method to a finite element analysis of non-homogeneous continuum problems continuum problems, see (Cite to Figure 1) outline. We assume in the latter that the non-homogeneous field of displacements (via DIC), plus the test loads (via load cell) are known, the input data could be either 1D, 2D or 3D depending on the case, but it is important to note that in 2D and 1D, it should be possible to employ reasonable assumptions to determine the principal stretches and stresses in the eliminated directions (incompressibility plus plane stress allow to determine both the stretch and the stress out of the plane just from the plane information). Then, the polymer chain behavior is modeled by P-splines, which vertices are to be determined—P-splines are penalyzed interpolating B-splines to guarantee smoothness; see Eilers and Marx (2021). That structure (the unknown vertices) are transferred to the continuum scale via integration in all the material directions and the result attached to the finite element integration point (the continuum constitutive behavior). Hence, the nodal forces of the finite element are set as a direct, explicit function of the unknown P-spline vertices of the polymer chains and the prescribed deformation gradient. By a least squares formulation, a linear system of equations is established, which allows for the immediate determination (i.e., simply solving a linear system of equations) of the P-spline vertices of the chain behavior from the macroscopic loads in the specimen and the macroscopic field of deformation. The physics equations present on the procedure include at the FEM level the compatibility equations (computing the strain quantities from the displacement field) and the equilibrium equations (null force residual), incompressible hyperelasticity with volumetric-deviatoric decoupling at the integration point level, properly including the micro-macro connections (energy homogenization and affine micro-stretch) at the chain level.

FIGURE 1

FIGURE 1. On the right side the FEM part (green contour), from bottom to top there are different levels, first at the solid level, the nodal displacement vector, u, the nodal internal forces, f_int, and the external nodal forces, f_ext, are encountered. Second, at the element level, the elemental nodal displacement vector, u^e, the elemental nodal internal forces, ${\mathit{f}}_{\text{int}}^e$, and the elemental external nodal forces, ${\mathit{f}}_{\text{ext}}^e$. Finally, the integration point (at the macro-scale) characterized by the Right Green Cauchy tensor C, and the PK2 stress tensor, S, at this point. On the left hand side the Macro-Micro-Macro approach (magenta contour), the general idea would be “zooming” (hence the magnifying glass) on the integration point and consider that for polymeric materials, the chains oriented in a certain direction (representative direction) suffer a deformation described by the continuum stretches ellipsoid (upper left). The a priori unknown chain mechanical behavior, $P_{ch}\left({\lambda }_{ch}\right)$ (bottom left), will be evaluated for each direction and the effect in each representative direction will be integrated for all the representative directions on the reference unit micro-sphere (center left) to obtain the continuum principal stress components. Finally, the Macro-Micro transition (center of the figure) dictates how the continuum energy is obtained from the constituents free energy and how the micro-structural stretches (non-affine chain stretch) are obtained from the continuum deformation measures. From a more practical perspective: 1) The Macro-Micro-Macro approach expresses the continuum principal stresses in terms of the unknown P-Spline vertices for the constituents mechanical behavior, $\mathit{S}\left(\widehat{\mathit{v}}\right)$, and 2) A classical FEM assembly process is employed to write nodal equilibrium equations from the information at the integration points, ${\mathrm{F}}_{\text{int}}\cdot \widehat{\mathit{v}}={\mathit{f}}_{\text{ext}}$. The substeps involved are: 1.1) The continuum strain, C, tensor defines the micro-structural deformations, λ_ch on the representative directions; 1.2) The micro-mechanical behavior (unknown a priori) is evaluated in all the representative directions in terms of the unknown vertices $\widehat{\mathit{v}}$, $P_{ch}\left({\lambda }_{ch}\right)={\mathit{B}}_{row}\left(x\right)\cdot \widehat{\mathit{v}}$; 1.3) The constituents mechanical behavior in all directions is integrated over the unit micro-sphere to obtain the continuum principal stresses at the integration point, $\mathit{S}\left(\widehat{\mathit{v}}\right)$; 2.1) The integration point stress $\mathit{S}\left(\widehat{\mathit{v}}\right)$ is integrated over the element volume to obtain the elemental nodal internal forces, ${\mathit{f}}_{\text{int}}^e\left(\widehat{\mathit{v}}\right)$; 2.2) The elemental nodal internal forces are assembled to obtain the nodal internal forces of the solid, ${\mathit{f}}_{\text{int}}\left(\widehat{\mathit{v}}\right)$; 2.3) The unknown vertices are determined from equilibrium, ${\mathit{f}}_{\text{int}}\left(\widehat{\mathit{v}}\right)={\mathrm{F}}_{\text{int}}\cdot \widehat{\mathit{v}}={\mathit{f}}_{\text{ext}}$.

In the following sections we introduce the procedure, first using a continuum hyperelastic formulation and then the micromechanical one. We also demonstrate the applicability through an academic example and an example using the well-known Treloar’s rubber.

2 Methodology

2.1 Detailed Procedure Description

One way to obtain the displacement-based finite element formulation is through the principle of virtual work, which for the quasi-static case reads δW = δW_int − δW_ext = 0. Considering a conventional FEM discretization and interpolation of the displacement field, both virtual works (internal and external) could also be expressed in terms of the internal and external nodal forces, $\delta W=\delta \mathit{u}\cdot \left({\mathit{f}}_{\text{int}}-{\mathit{f}}_{\text{ext}}\right)=0$ or equivalently f_int = f_ext. To conclude, a weighted integration of the equilibrium equation leads to an alternative expression for the virtual work, see (Eq. 1), that if compared with $\delta W=\delta \mathit{u}\cdot \left({\mathit{f}}_{\text{int}}-{\mathit{f}}_{\text{ext}}\right)=0$ provides an expression for the internal and external nodal forces. The expression for the weighted integration of the equilibrium equation in the reference configuration can be expressed as:

${\int }_{\mathrm{\Omega }}\left[\delta \mathit{u}\otimes \nabla \right]:\mathit{P}\left(\mathit{u}\right)d\mathrm{\Omega }={\int }_{\mathrm{\Omega }}\delta \mathit{u}\cdot \mathit{b}d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\delta \mathit{u}\cdot \mathit{t}d\mathrm{\Gamma }(1)$

where δu and u are the virtual displacement field, and the displacement field, respectively, P, the 1st PK (First Piola-Kirchhoff) stress tensor, b the volumetric forces and t the surface forces. The symbol ⊗ represents the dyadic production, so $\left[\delta \mathit{u}\otimes \nabla \right]:\mathit{P}\left(\mathit{u}\right)=\mathit{P}:\nabla \delta \mathit{u}$ is the internal virtual work density.

2.1.1 Internal Force Term

As it has already been mentioned, a FEM discretization for the internal force term, ∫_ΩP: ∇δudΩ, and the interpolation of the displacement field in the reference unit element, $\delta \mathit{u}={\sum }_{a=1}^{n_n}{h}^a\left({}^0\mathit{\xi }\right)\delta {\mathit{u}}^a$ leads to the expression for the nodal internal forces, see (Eq. 13). Note that $h^a\left({}^0\mathit{\xi }\right)$ are the shape functions, δu^a the virtual displacement vector for the node a, and that in $\delta u_i^a$, a = 1, … , n_n and i indicates the spatial dimension (i = 1, 2 in 2D). The expressions for the internal forces are described with depth underneath, (Eq. 2) and (Eq. 3):

$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}\mathit{P}:\nabla \delta \mathit{u}d\mathrm{\Omega }& =\sum _{e=1}^{n_{el}}{\int }_{\mathrm{\Omega }e}\mathit{P}:\nabla \delta \mathit{u}d{\mathrm{\Omega }}_e\hfill \\ \hfill & =\sum _{e=1}^{n_{el}}{\int }_{\square }\mathit{P}:\nabla \delta \mathit{u}{J}^{e}d\square \hfill \\ \hfill & =\sum _{e=1}^{n_{el}}\sum _{j=1}^{n_{qp}}{\mathit{P}}_j:\nabla \delta {\mathit{u}}_j{J}_j^e{w}_j\hfill \\ \hfill & =\sum _{e=1}^{n_{el}}\sum _{j=1}^{n_{qp}}\sum _{a=1}^{n_n}{\mathit{P}}_j:\left(\delta {\mathit{u}}^a\otimes \nabla h^a\left({}^0\mathit{\xi }_j\right)\right)J_j^e{w}_j\hfill \end{array}(2)$

where all the variables with subscript j are computed in the integration point j. The previous equations can be rewritten doing the sum over the DOF of the element (n_dofs = n_n × 2 in the 2D solid elements) instead of doing the sum over the nodes:

${\int }_{\mathrm{\Omega }}\mathit{P}:\nabla \delta \mathit{u}d\mathrm{\Omega }=\sum _{e=1}^{n_{el}}\sum _{j=1}^{n_{qp}}\sum _{i=1}^{n_{\text{dofs}}}\delta u_i{\mathit{P}}_j:\nabla {\mathit{h}}_j^i{J}_j^e{w}_j,(3)$

In Eq. 3, i is the index that runs through the local degrees of freedom in the element, for the cases studied here (2D plane stress problems) i = 1, … , 2n_n, δu_i is the virtual displacement at the local degree of freedom i and $\nabla {\mathit{h}}_j^i$ is a second order tensor that projects the contribution of P_j to the ith local degree of freedom. Regarding $\nabla {\mathit{h}}_j^i$, for i = 2, $\nabla {\mathit{h}}_j^2$ is the projector for the second local degree of freedom of the element, this DOF corresponds to the first node of the element a = 1 and the second dimension 2, therefore, $\nabla {\mathit{h}}_j^2={\mathit{e}}_2\otimes \nabla h^1\left({}^0\mathit{\xi }_j\right)$. If i = 3 instead, $\nabla {\mathit{h}}_j^3$ is the projector for the third local degree of freedom of the element, which corresponds to the second node of the element a = 2 and the first dimension 1, therefore, $\nabla {\mathit{h}}_j^3={\mathit{e}}_1\otimes \nabla h^2\left({}^0\mathit{\xi }_j\right)$. Now Eq. 3 is rewritten in matrix form to identify the components of the internal force vector:

$\begin{array}{cc}\hfill \sum _{e=1}^{n_{el}}\sum _{j=1}^{n_{qp}}\sum _{i=1}^{n_{\text{dofs}}}\delta u_i{\mathit{P}}_j:\nabla {\mathit{h}}_j^i{J}_j^e{w}_j& =\sum _{e=1}^{n_{el}}\left[\begin{array}{ccc}\hfill {\left(f_{\text{int}}^e\right)}_1\hfill & \hfill \dots \hfill & \hfill {\left(f_{\text{int}}^e\right)}_{n_{\text{dofs}}}\hfill \end{array}\right]\left[\begin{array}{c}\hfill \delta u_1^e\hfill \\ \hfill ⋮\hfill \\ \hfill \delta u_{n_{\text{dofs}}}^e\hfill \end{array}\right]\hfill \\ \hfill & =\sum _{e=1}^{n_{el}}{\mathit{f}}_{\text{int}}^e\cdot \delta {\mathit{u}}^e\hfill \end{array}(4)$

where

${\left(f_{\text{int}}^e\right)}_i=\sum _{j=1}^{n_{qp}}{\mathit{P}}_j:\nabla {\mathit{h}}_j^i{J}_j^e{w}_j(5)$

Since the PK1 tensor is a two-leg tensor placed in 2 configurations at the same time (material in the right and spatial in the left), it might be more suitable to rewrite the term ${\mathit{P}}_j:\nabla {\mathit{h}}_j^i$ in terms of the PK2 (Second Piola–Kirchhoff) stress tensor, S, which lies completely in material configuration, P = XS, being X = ∂^tx/∂⁰x the deformation gradient:

$\begin{array}{c}{\mathit{P}}_j:\nabla {\mathit{h}}_j^i=\left({\mathit{X}}_j{\mathit{S}}_j\right):\nabla {\mathit{h}}_j^i={\mathit{S}}_j:{\mathit{X}}_j^T\nabla {\mathit{h}}_j^i={\mathit{S}}_j:\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{X}}_j^T\nabla {\mathit{h}}_j^i\right)\\ =\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{X}}_j^T\nabla {\mathit{h}}_j^i\right):{\mathit{S}}_j\end{array}(6)$

In order to simplify the notation we define ${\mathit{l}}_j^i≔\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{X}}_j^T\nabla {\mathit{h}}_j^i\right)$, on the other hand the double contraction will be computed using the Voigt notation. Note that S₃ = 0 since we are considering the case of plane stress and that the components in the Voigt notation are with respect to the principal directions of deformation/stress. Then, we write in principal directions (denoted as the X_ppal system of representation)

${\mathit{P}}_j:\nabla {\mathit{h}}_j^i={\left[\begin{array}{ccc}\hfill l_{11}\hfill & \hfill l_{22}\hfill & \hfill 2l_{12}\hfill \end{array}\right]}_{j,X_{\text{ppal}}}^i{\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill 0\hfill \end{array}\right]}_{j,X_{\text{ppal}}}(7)$

To compute the principal stress components of a polymeric and quasi-incompressible material, we assume that the volumetric and deviatoric contributions of the energy can be separated (a typical assumption in quasi-incompressibility), $\mathrm{\Psi }\left(\mathit{A}\right)=\mathcal{U}\left(J\right)+\mathcal{W}\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)$, where the volumetric term $\mathcal{U}\left(J\right)$ is just a penalization function, being $J=\det \left(\mathit{X}\right)=d^{t}V/d^{0}V$ the volume ratio, and $\mathit{A}=1/2\left({\mathit{X}}^T\mathit{X}-\mathit{I}\right)$ the Green-Lagrange strain tensor. A typical choice for the penalization could be $\mathcal{U}\left(J\right)=\kappa /2{\left(J-1\right)}^2$ with the bulk modulus, κ, selected such that the stress produced by a volumetric deformation grows rapidly. If an estimation for the shear modulus of the material, μ, is known, κ/μ∼10⁴, might suffice to ensure the satisfaction of the quasi-incompressibility condition J = λ₁ ⋅ λ₂ ⋅ λ₃ ≈ 1. The decoupling of the energy in a volumetric and deviatoric energy results in a similar decoupling for the PK2 stress tensor, S = S^v + S^d. We note that S^v and S^d are not the volumetric and deviatoric part, respectively, of a second order tensor S (as obtained from the respective mathematical operators), but S^v is the contribution to S that comes from the volumetric energy $\mathcal{U}(J)$, and S^d is the contribution to S that comes from the deviatoric energy, $\mathcal{W}({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d)$:

$\mathit{S}=\frac{d\mathrm{\Psi }}{d\mathit{A}}=\frac{d\mathcal{U}}{dJ}\frac{dJ}{d\mathit{A}}+\sum _{i=1}^3\frac{\partial \mathcal{W}}{\partial {\lambda }_i^d}\frac{d{\lambda }_i^d}{d\mathit{A}}(8)$

$={\mathit{S}}^v+{\mathit{S}}^d=pJ{\mathit{C}}^{-1}+\sum _i\sum _k\frac{{\mathcal{W}}_k{J}^{-1/3}}{{\lambda }_i^2}\left({\lambda }_i{\delta }_{ik}-\frac{1}{3}{\lambda }_k\right){\mathit{N}}_i\otimes {\mathit{N}}_i(9)$

where N_i are the principal referential directions of deformation, ${\mathcal{W}}_k=\partial \mathcal{W}/\partial {\lambda }_k^d$ and ${\lambda }_i^d={\lambda }_i{J}^{-1/3}$, the isochoric stretches. As we have already mentioned the test cases are plane stress and therefore, $S_3=0=\mathit{S}:\left({\mathit{N}}_3\otimes {\mathit{N}}_3\right)$:

$S_3=\frac{pJ}{{\lambda }_3^2}+\frac{2}{3}\frac{{\mathcal{W}}_3{J}^{-1/3}}{{\lambda }_3}-\frac{1}{3}\left(\frac{{\mathcal{W}}_1{J}^{-1/3}}{{\lambda }_3^2}{\lambda }_1+\frac{{\mathcal{W}}_2{J}^{-1/3}}{{\lambda }_3^2}{\lambda }_2\right)=0(10)$

From the previous equation, the pressure term could be isolated and introduced in the equations for S₁ and S₂:

$S_j=J^{-2/3}\frac{1}{{\lambda }_j^d}\left(\frac{\partial \mathcal{W}}{\partial {\lambda }_j^d}-\frac{{\lambda }_3^d}{{\lambda }_j^d}\frac{\partial \mathcal{W}}{\partial {\lambda }_3^d}\right)\mathrm{f}\mathrm{o}\mathrm{r}j=\mathrm{1,2}(11)$

Once the principal stretches are known for a particular integration point, the principal PK2 stress components in that same point could be determined through evaluation of the deviatoric contribution derivatives. When the stress tensor has been computed in all the integration points of the element, the internal nodal forces for that element are computed using (Eq. 4). On the other hand, if the functions for the deviatoric contribution derivatives are unknown, a cubic P-Spline local approximation can be employed, see Amores et al. (2019), Eilers and Marx (2021). The B-splines (or its penalized version, P-splines) are one of the approaches for expressing any general function y(x) as the product of a set of known basis functions (cubic in this case), $B_i\left(x\right)$ and it’s a priori unknown corresponding weights (also called vertices), ${\widehat{v}}_i$, if the number of vertices of the B-Spline is n_vert, the expression for the unknown function is a scalar product, $y\left(x\right)={\sum }_{i=1}^{n_{\text{vert}}}{B}_i\left(x\right){\widehat{v}}_i=\mathit{B}\left(x\right)\cdot \widehat{\mathit{v}}$. In contrast to the basis of features proposed in Flaschel et al. (2021), where a basis of preassumed functions is used, the local B-Spline basis is general with no assumptions about the possible function expressions (the only assumption would be that locally the function is at most a cubic function when cubic B-Splines are used). Additionally, the local P-Spline approximation in this work is typically performed directly on the derivatives so there is no need to compute derivatives of the approximated function. With the mentioned P-Spline for the derivative of the energy function, the principal stresses could then be represented as $S_i\left({\lambda }_1^d,{\lambda }_2^d\right)={\mathit{S}}_{\text{row}}^i\left({\lambda }_1^d,{\lambda }_2^d\right)\cdot \widehat{\mathit{v}}$ where ${\mathit{S}}_{\text{row}}^i$ is a vector defined with 1D P-Splines basis vector ($\mathit{B}\left({\lambda }_1^d\right)$, $\mathit{B}\left({\lambda }_2^d\right)$), see Section 3.1 and Section 3.2 for specific expressions of this vector. Therefore, depending on the approach that is used to compute the continuum principal stresses, the unknown vertices, $\widehat{\mathit{v}}$, could represent the derivative of the energy either on the macro-structure or in the micro-structure:

$\begin{array}{ll}\hfill {\mathit{P}}_j:\nabla {\mathit{h}}_j^i& =\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{X}}_j^T\nabla {\mathit{h}}_j^i\right):\hfill \\ {\mathit{S}}_j\hfill & ={\left[\begin{array}{ccc}\hfill l_{11}\hfill & \hfill l_{22}\hfill & \hfill 2l_{12}\hfill \end{array}\right]}_{j,X_{\text{ppal}}}^i{\left[\begin{array}{c}\hfill {\left[{\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_2^d\right)\right]}^T\hfill \\ \hfill {\left[{\mathit{S}}_{\text{row}}^2\left({\lambda }_1^d,{\lambda }_2^d\right)\right]}^T\hfill \\ \hfill {\left[\mathbf{0}\right]}^T\hfill \end{array}\right]}_j\left[\widehat{\mathit{v}}\right]\hfill \\ \hfill & ={\left(l_{11}{\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_2^d\right)+l_{22}{\mathit{S}}_{\text{row}}^2\left({\lambda }_1^d,{\lambda }_2^d\right)\right)}_j^i\cdot \widehat{\mathit{v}}\hfill \\ \hfill & ={\left({}^{*}\mathit{S}_{\text{row}}\right)}_j^i\left({\lambda }_1^d,{\lambda }_2^d\right)\cdot \widehat{\mathit{v}}\hfill \end{array}(12)$

and for the element local degree of freedom i

${\left(f_{\text{int}}^e\right)}_i=\left[\sum _{j=1}^{n_{qp}}{J}_j^e{w}_j{\left({}^{*}\mathit{S}_{\text{row}}\right)}_j^i\left({\lambda }_1^d,{\lambda }_2^d\right)\right]\cdot \widehat{\mathit{v}}(13)$

which means that every component of the elemental nodal force vector is obtained as the product of a known row multiplied by the unknown vertices, or expressed in a different manner:

${\mathit{f}}_{\text{int}}^e={\mathrm{F}}_{\text{int}}^e\cdot \widehat{\mathit{v}}(14)$

where ${\mathrm{F}}_{\text{int}}^e$ is a known matrix, given in square brackets in Eq. 13.

2.1.2 External Forces Term

In case the body forces and the tractions are not zero, the external force vector for each element, ${\mathit{f}}_{\text{ext}}^e$ has also to be computed and assembled into f_ext, looking at the right hand side of (Eq. 1).

${\left(f_{\text{ext}}^e\right)}_i=\sum _{j=1}^{n_{qp}}{\mathit{h}}_j^i\cdot \mathit{b}{J}_j^e{w}_j+\sum _{s=1}^{n_{sbound}^e}\sum _{j=1}^{n_{qps}}{\bar{\mathit{h}}}_j^i\cdot \mathit{t}{J}_j^s{w}_j^s(15)$

${\mathit{f}}_{\text{ext}}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{f}}_{\text{ext}}^e(16)$

On the previous equations, the reaction forces on the boundary where displacements are imposed are accounted in the traction term, for the sake of simplicity we are going to consider that those nodal reaction forces are known, while this is not typically the case in a experimental setting, instead, the total reaction forces are known rather than the nodal forces. To deal with this fact, 2 different approaches can be followed: 1) take another artificial boundary far enough from the original one in which we can suppose that the reaction force is evenly distributed according to the Saint-Venant’s principle or 2) consider two independent set of equations one for the free dofs and other for the DOF with fixed displacements, in the fixed displacement DOF the resultant of the internal forces equals the reaction force at each of the boundaries.

2.1.3 System of Linear Equations for the P-Spline Vertices

The problem to solve is an overdetermined linear system of equations ${\mathrm{F}}_{\text{int}}\cdot \widehat{\mathit{v}}={\mathit{f}}_{\text{ext}}$, where $\widehat{\mathit{v}}$ are the unknown vertices of the P-Spline that approximates the derivative of the energy function at the macro-scale (phenomenological energy function) or the derivative of the energy function of the constituents at the micro-scale. The solution process consists of finding the vertices, $\widehat{\mathit{v}}$, that minimizes the mean square error, MSE, in the force residual:

$MSE\left(\widehat{\mathit{v}}\right)=\left({\mathrm{F}}_{\text{int}}\cdot \widehat{\mathit{v}}-{\mathit{f}}_{\text{ext}}\right)\cdot \left({\mathrm{F}}_{\text{int}}\cdot \widehat{\mathit{v}}-{\mathit{f}}_{\text{ext}}\right)(17)$

The solution for the previous minimization problem is analytical and result in the mentioned linear system of equations:

$\widehat{\mathit{v}}={\mathit{A}}_{\text{sys}}^{-1}{\mathit{b}}_{\text{sys}};{\mathit{A}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T{\mathrm{F}}_{\text{int}},{\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T{\mathit{f}}_{\text{ext}}(18)$

As we have already mentioned, one of the advantages of using the P-Spline-based local approximation is that although we do not assume the form of the energy function, additional requirements can be added to the solution, a typical one is smoothness, which can be translated into a penalization on the second order finite differences of the P-spline vertices, if the solution obtained is not monotonically increasing, this property could also be imposed by iterative penalization on the vertices that do not satisfy the condition see Amores et al. (2019), Eilers and Marx (2021). With all the penalizations, the general system of equations to solve would be:

$\begin{array}{c}{\widehat{\mathit{v}}}^k={\left({\mathit{A}}_{\text{sys}}^k\right)}^{-1}{\mathit{b}}_{\text{sys}};\\ {\mathit{A}}_{\text{sys}}^k={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathrm{F}}_{\text{int}}+{\mathit{D}}_2^T{\mathbf{\Omega }}_2{\mathit{D}}_2+{\mathit{D}}_1^T{\mathbf{\Omega }}_1^k{\mathit{D}}_1,\\ {\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathit{f}}_{\text{ext}},\end{array}(19)$

where W is a diagonal matrix that weights the relative importance of the equations in ${\mathrm{F}}_{\text{int}}$, D₁ and D₂ are the first and second differences matrices respectively, Ω₂ and Ω₁, are also diagonal weight matrices for the penalizations in the second and first differences of the vertices, note that k indicates the step of the iterative penalization for monotonic smoothing. When all the intervals are monotonically increasing the weighting matrix ${\mathbf{\Omega }}_1^k=\mathbf{0}$ and that penalization automatically disappears. In contrast to the overdetermined system of linear equations produced by piecewise approximation presented in this work, the methodology in Flaschel et al. (2021) employs a basis of functions which extend over the whole function domain, although this procedure seems similar to the local one, the global support typically leads to dense solutions which are less interpretable. For the sake of interpretability Flaschel et al. (2021) proposes a sparsity promotion with ℓ_p regularization, but at the price of requiring a fixed point iteration procedure to obtain the solution and introducing p as an hyperparameter that has also to be determined in the solution process. On the other hand, the additional matrices added to the system matrix in our work ensure smoothness of the solution (making the solution more robust against noise) and allows penalization in monotonicity if required (this penalization has not been required for the examples presented in this article).

3 Results

3.1 Toy Example: Recovering an Analytical Neo-Hookean With a P-Spline for the Continuum Free Energy

In this section we are going to demonstrate just the FEM part of the methodology using a P-Spline based approximation for the functions that form the deviatoric continuum free energy. The experimental data is a virtual test (FEM results) of an analytical quasi-incompressible Neo-Hookean model with parameters μ = 3.5 MPa and κ = 1 × 10⁵ MPa. The FEM model employed for the virtual tests is a 2D plate with a hole on plane stress conditions. The solicitation is an imposed displacement in the upper border with u = 4. The dimensions of the plate and the mesh employed in the simulation is shown in Figure 2, the elements are quadratic quadrilaterals with nine integration points. The formulation employed for plane stress in large deformations is detailed in Supplementary Appendix SA, the reader can also find the details about the Neo-Hookean material model in Supplementary Appendix SB. Regarding the software employed, Julia programming language, Bezanson et al. (2017) has been used, in particular the package FerriteFem.jl, see Carlsson and Ekre (2021), for the FEM simulations and Amores (2022) for the P-Splines functionality.

FIGURE 2

FIGURE 2. FEM model used during virtual tests. The elements employed are quadratic quadrilateral with nine integration points.

A typical assumption made on the phenomenological approach for isotropic incompressible solids is the Valanis Landel decomposition, $\mathcal{W}\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)=\omega \left({\lambda }_1^d\right)+\omega \left({\lambda }_2^d\right)+\omega \left({\lambda }_3^d\right)$. With that assumption, the expressions for the principal PK2 stresses can be obtained in terms of an unknown function ${\omega }^{\prime }\left(x\right)$:

$S_j=J^{-2/3}\frac{1}{{\lambda }_j^d}\left({\omega }^{\prime }\left({\lambda }_j^d\right)-\frac{{\lambda }_3^d}{{\lambda }_j^d}{\omega }^{\prime }\left({\lambda }_3^d\right)\right)\mathrm{f}\mathrm{o}\mathrm{r}j=\mathrm{1,2}(20)$

Once ${\omega }^{\prime }\left(x\right)$ is obtained, the principal stresses are determined for any deformation state under the plane stress hypothesis or the contribution to the principal stresses coming from the deviatoric energy for any general case (not in plane stress), see (Eq. 8), this is all that can be determined from a constitutive point of view because for an incompressible solid, the pressure does not come from the constitutive equation but from equilibrium conditions. Using P-splines, the unknown function can be written in terms of some known basis functions and some unknown vertices in a similar manner that it is done in FEM formulation with the displacement field:

$\begin{array}{c}{\omega }^{\prime }\left(x\right)=\sum _{i=1}^{n_{\text{vert}}}{B}_i\left(x\right){\widehat{\omega }}_i=\left[\begin{array}{ccc}\hfill B_1\left(x\right)\hfill & \hfill \dots \hfill & \hfill B_{n_{\text{vert}}}\left(x\right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\widehat{\omega }}_1\hfill \\ \hfill ⋮\hfill \\ \hfill {\widehat{\omega }}_{n_{\text{vert}}}\hfill \end{array}\right]\\ ={\mathit{B}}_{row}\left(x\right)\cdot {\widehat{\mathit{\omega }}}^{\prime }\end{array}(21)$

With the previous expression, the principal stresses could be written as:

$\begin{array}{cc}& S_j=J^{-2/3}\frac{1}{{\lambda }_j^d}\left({\mathit{B}}_{row}\left({\lambda }_j^d\right)-\frac{{\lambda }_3^d}{{\lambda }_j^d}{\mathit{B}}_{row}\left({\lambda }_3^d\right)\right)\cdot {\widehat{\mathit{\omega }}}^{\prime }={\mathit{S}}_{\text{row}}^j\left({\lambda }_j^d,{\lambda }_3^d\right)\cdot {\widehat{\mathit{\omega }}}^{\prime }\\ & \mathrm{f}\mathrm{o}\mathrm{r}j=\mathrm{1,2}\hfill \end{array}(22)$

${\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill 0\hfill \end{array}\right]}_{j,X_{\text{ppal}}}={\left[\begin{array}{c}\hfill {\left[{\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_3^d\right)\right]}^T\hfill \\ \hfill {\left[{\mathit{S}}_{\text{row}}^2\left({\lambda }_2^d,{\lambda }_3^d\right)\right]}^T\hfill \\ \hfill {\left[\mathbf{0}\right]}^T\hfill \end{array}\right]}_j{\left[{\widehat{\mathit{\omega }}}^{\prime }\right]}^T(23)$

$\begin{array}{l}\hfill {\mathit{P}}_j:\nabla {\mathit{h}}_j^i=\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{X}}_j^T\nabla {\mathit{h}}_j^i\right):{\mathit{S}}_j={\left[\begin{array}{ccc}\hfill l_{11}\hfill & \hfill l_{22}\hfill & \hfill 2l_{12}\hfill \end{array}\right]}_{j,X_{\text{ppal}}}^i{\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill 0\hfill \end{array}\right]}_{j,X_{\text{ppal}}}\\ \hfill ={\left(l_{11}\cdot {\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_3^d\right)+l_{22}\cdot {\mathit{S}}_{\text{row}}^2\left({\lambda }_2^d,{\lambda }_3^d\right)\right)}_j^i\cdot {\widehat{\mathit{\omega }}}^{\prime }\\ ={\left({}^{*}\mathit{S}_{\text{row}}\right)}_j^i\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)\cdot {\widehat{\mathit{\omega }}}^{\prime }\end{array}(24)$

Referring again to Eq. 13, internal nodal forces of the element can be written as the product of a known matrix by a vector of unknown vertices:

${\mathit{f}}_{\text{int}}^e={\mathrm{F}}_{\text{int}}^e\cdot {\widehat{\mathit{\omega }}}^{\prime }(25)$

Doing the assembly of the internal nodal forces of the elements in the mesh:

${\mathit{f}}_{\text{int}}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{f}}_{\text{int}}^e=\left(\underset{e=1}{\overset{nel}{\bigwedge }}{\mathrm{F}}_{\text{int}}^e\right)\cdot {\widehat{\mathit{\omega }}}^{\prime }={\mathrm{F}}_{\text{int}}\cdot {\widehat{\mathit{\omega }}}^{\prime }(26)$

From the overdetermined linear system of equations ${\mathrm{F}}_{\text{int}}\cdot {\widehat{\mathit{\omega }}}^{\prime }={\mathit{f}}_{\text{ext}}$ we can solve for ${\widehat{\mathit{\omega }}}^{\prime }$:

$\begin{array}{c}{\left({\widehat{\mathit{\omega }}}^{\prime }\right)}^k={\left({\mathit{A}}_{\text{sys}}^k\right)}^{-1}{\mathit{b}}_{\text{sys}};\\ {\mathit{A}}_{\text{sys}}^k={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathrm{F}}_{\text{int}}+{\mathit{D}}_2^T{\mathbf{\Omega }}_2{\mathit{D}}_2+{\mathit{D}}_1^T{\mathbf{\Omega }}_1^k{\mathit{D}}_1,\\ {\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathit{f}}_{\text{ext}}\end{array}(27)$

The number of vertices has to be enough to capture the complexity of the curve, typically n_vert = 14 suffice, but additional vertices can be added. If more and more vertices are added, the number of equations required to determine them increases and the problem can become ill-conditioned, this is solved by the smoothing term which adds the information of smooth transition between vertices and links the vertices to its neighbours. In the homogeneous case, just a single equation is obtained for every load/displacement step, therefore, in order to obtain information for the function on the considered domain (from the minimum principal stretch to the maximum principal stretch on the simulation), it would be necessary to sweep a whole range of load/displacement steps. On the other hand, for the non-homogeneous case, in principle, it would be possible to employ just a single load/displacement step if the step under consideration is rich enough (in this case, just the last step, u = 4 was used). In case that additional information is needed to determine the function on the considered range, it is also possible to add more steps between u = 0 and u = 4. Since the initial solution was directly monotonically increasing, Ω₁ = 0 and it was not necessary to follow an iterative process, ${\widehat{\mathit{\omega }}}^{\prime }$ was directly obtained from the simple initial system of equations:

$\left({\widehat{\mathit{\omega }}}^{\prime }\right)={\left({\mathit{A}}_{\text{sys}}\right)}^{-1}{\mathit{b}}_{\text{sys}};{\mathit{A}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathrm{F}}_{\text{int}}+{\mathit{D}}_2^T{\mathbf{\Omega }}_2{\mathit{D}}_2,{\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathit{f}}_{\text{ext}}(28)$

With the vertices obtained, the P-Spline could be reconstructed and compared to the original one, see Figure 3.

FIGURE 3

FIGURE 3. Comparison of the Neo-Hookean ${\omega }_{\text{NH}}^{\prime }\left({\lambda }^d\right)$ and the ${\omega }^{\prime }\left({\lambda }^d\right)$ obtained from the non-homogeneous virtual test.

Since the function ${\omega }^{\prime }\left(x\right)$ has been reverse-engineered exactly, the values for the free energy partial derivatives for any state of deformation can be determined, and will be exactly equal to those obtained with the original Neo-Hookean model. If desired, the FEM simulations can be run again with the reversed-engineered derivatives of the energy function and the σ_VM plot could be compared with the one obtained for the original analytical Neo-Hookean, see Figure 4. As it can be seen in Figure 4 the plots are exactly the same as it might be expected from having obtained the exact same ω′ as the one from the original Neo-Hookean.

FIGURE 4

FIGURE 4. Plot comparing the results of σ_VM for the reference Neo-Hookean continuum model and the reverse-engineered one. (A) σ_VM using the Neo-Hookean model. (B) σ_VM using the ${\omega }^{\prime }\left({\lambda }_i^d\right)$ reverse-engineered from the non-homogeneous virtual test.

In this particular example, the matrix of the system ${\mathit{A}}_{\text{sys}}\cdot {\widehat{\mathit{\omega }}}^{\prime }={\mathit{b}}_{\text{sys}}$, presents a null eigenvalue with its corresponding associated eigenvector, i.e., there is a certain subspace of possible solutions ${\widehat{\mathit{\omega }}}_0^{\prime }$ (containing just the direction of the eigenvector associated to the null eigenvalue) such that added to the solution, ${\widehat{\mathit{\omega }}}^{\prime }$, produce null effect on the independent term, b_sys and therefore ${\widehat{\mathit{\omega }}}^{\prime }+\alpha {\widehat{\mathit{\omega }}}_0^{\prime }$, for an arbitrary scalar α, is also a solution of the system. Looking at Eq. 20, to not produce any effect, such function has to satisfy the condition ${\omega }^{\prime }\left({\lambda }_i^d\right)-{\lambda }_3^d/{\lambda }_i^d\cdot {\omega }^{\prime }\left({\lambda }_3^d\right)=0\forall {\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d$. Regarding the functions that satisfy this condition, there is a specific form of function that complies (which at the same time must be smooth due to the regularization) ${\omega }^{\prime }\left({\lambda }_i^d\right){\lambda }_i^d={\omega }^{\prime }\left({\lambda }_3^d\right){\lambda }_3^d\forall {\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d$. The only way two functions of different independent variables can be equal $\forall {\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d$ is that those functions are constant ${\omega }^{\prime }\left({\lambda }_i^d\right){\lambda }_i^d={\omega }^{\prime }\left({\lambda }_3^d\right){\lambda }_3^d={\beta }_0\to {\omega }^{\prime }\left(\lambda \right)={\beta }_0/\lambda$. Since the shape of the eigenvector associated to the null eigenvalue exactly matches the shape of the function, ${\omega }^{\prime }\left(\lambda \right)={\beta }_0/\lambda$, we would like to understand physically the origin of this null contribution to the system. Looking at Eq. 8, we can see that the pressure term is of the form $p\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)/{\lambda }_i^2$, the same form of the term resulting from introducing an additional ${\beta }_0/{\lambda }_i^d$ to the function ${\omega }^{\prime }\left({\lambda }_i^d\right)$. From a more general perspective we can see that this indetermination comes from the split of the energy into a volumetric and a deviatoric contribution:

$\mathit{S}={\mathit{S}}^v+{\mathit{S}}^d;{\mathit{S}}^d=2\frac{d\mathcal{W}}{d\mathit{C}}=2J^{-2/3}\mathbb{P}:\frac{d\mathcal{W}}{d{\mathit{C}}^d}=J^{-2/3}\mathbb{P}:{\mathit{S}}^{\mid d}(29)$

where C = X^TX is the right Cauchy-Green tensor and $\mathbb{P}$ is the fourth order projector tensor, see Supplementary Appendix SB. Since $\mathbb{P}:\left(\alpha {\mathit{C}}^{-1}\right)=0$, from (Eq. 29) it can be seen that there are multiple ${\mathit{S}}^{\mid d}=\partial \mathcal{W}/\partial {\mathit{C}}^{\mathit{d}}$ that produce the same S^d and the undetermined part has the form of a volumetric-like PK2 stress tensor. Any contribution of this form has no effect on the solution since it is going to be finally eliminated by the projector $\mathbb{P}$.

3.2 Micro-Mechanical Approach Based on the Non-Affine Deformation Chain Model

In this section the complete methodology described in Section 2 (Macro-Micro-Macro first, FEM second) is demonstrated. In contrast to the procedure described in Section 3.1, here the P-Spline will approximate the derivative of the chain energy function $P_{ch}\left({\lambda }_{ch}\right)=d{\psi }_{ch}/d{\lambda }_{ch}$ (mechanical behavior of the material constituents at the micro-scale) rather than the continuum contributions to the energy. The experimental data is obtained from a virtual (FEM simulation) using the analytical structure-based material model in Amores et al. (2020). To be more precise, the $P_{ch}\left({\lambda }_{ch}\right)$ that will be employed to perform the virtual test is the one obtained for the Treloar test data for unfilled rubber, Treloar (1944), that function is displayed in (Amores et al., 2020, Figure B.1). The FEM model employed for the virtual tests is again the 2D plate with a hole on plane stress conditions with an imposed displacement in the upper border (u = 8). All the required dimensions and the mesh employed in the simulation is shown in Figure 2. The reader can find further details about the non-affine deformation chain material model in Supplementary Appendix SC and Amores et al. (2020), a discussion about the suitability of the non-affine stretch employed in the model (free-fluctuating network assumption) is presented in Amores et al. (2021).

In Amores et al. (2020) a way was established to compute the strain energy function of the continuum with an homogenization of the micro-structural chain free energy function in polymeric like materials under the assumption of incompressibility (λ₁ ⋅ λ₂ ⋅ λ₃ = 1):

$\mathcal{W}\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)=\frac{1}{S}{\int }_S{\psi }_{ch}\left({\lambda }_{ch}\right)dS(30)$

For FEM simulations even pure incompressible solids are simulated with the quasi-incompressible material framework:

$\mathrm{\Psi }\left(\mathit{A}\right)=\mathcal{U}\left(J\right)+\mathcal{W}\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)$

in which $\mathcal{U}\left(J\right)$ is a penalization function that leads to incompressibility for κ → ∞ like $\mathcal{U}\left(J\right)=\frac{1}{2}k{\left(J-1\right)}^2$ and where $\mathcal{W}\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)$ is the same as the one used in pure incompressibility but with the deviatoric stretches instead of the total ones:

$\begin{array}{c}\mathcal{W}\left({\lambda }_1^d,{\lambda }_2^d,{\lambda }_3^d\right)=\frac{1}{S}{\int }_S{\psi }_{ch}\left({\lambda }_{ch}^d\right)dS;\\ {\lambda }_{ch}^d={\mathit{U}}^d:\left(\mathit{r}\otimes \mathit{r}\right)={\lambda }_1^d{r}_1^2+{\lambda }_2^d{r}_2^2+{\lambda }_3^d{r}_3^2\end{array}(31)$

${\mathcal{W}}_k=\frac{\partial \mathcal{W}}{\partial {\lambda }_k^d}=\frac{1}{S}{\int }_S\frac{d{\psi }_{ch}\left({\lambda }_{ch}^d\right)}{d{\lambda }_{ch}^d}\frac{d{\lambda }_{ch}^d}{d{\lambda }_k^d}dS=\frac{1}{S}{\int }_S{P}_{ch}\left({\lambda }_{ch}^d\right)r_k^{2}dS(32)$

The unknown function that will be approximated using P-Splines is $P_{ch}\left({\lambda }_{ch}^d\right)$, as described above, the P-Splines representation expresses the unknown function in terms of some known basis vector multiplied by a vector of unknown vertices:

$\begin{array}{cc}{P}_{ch}\left(x\right)& =\sum _{i=1}^{n_{\text{vert}}}{B}_i\left(x\right){\widehat{P}}_{chi}\\ & =\left[\begin{array}{ccc}\hfill B_1\left(x\right)\hfill & \hfill \dots \hfill & \hfill B_{n_{\text{vert}}}\left(x\right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\widehat{P}}_{chi}\hfill \\ \hfill ⋮\hfill \\ \hfill {\widehat{P}}_{ch{n}_{\text{vert}}}\hfill \end{array}\right]\\ & ={\mathit{B}}_{row}\left(x\right)\cdot {\widehat{\mathit{P}}}_{ch}\end{array}(33)$

Introducing the expression for $P_{ch}\left(x\right)$ in (Eq. 32) and then all the derivatives of the form ${\mathcal{W}}_k$ in (Eq. 11) a expression for the principal stresses in terms of the vertices that define the Polymeric-chain response function can be obtained:

$\begin{array}{c}{S}_j=J^{-2/3}\frac{1}{{\lambda }_j^d}\left\{\frac{1}{S}{\int }_S{\mathit{B}}_{\text{row}}\left({\lambda }_{ch}^d\right)\left(r_j^2-\frac{{\lambda }_3^d}{{\lambda }_j^d}{r}_3^2\right)dS\right\}\cdot {\widehat{\mathit{P}}}_{ch}\\ ={\mathit{S}}_{\text{row}}^j\left({\lambda }_1^d,{\lambda }_2^d\right)\cdot {\widehat{\mathit{P}}}_{ch}\mathrm{f}\mathrm{o}\mathrm{r}j=\mathrm{1,2}\end{array}(34)$

Note that the previous integral in the microsphere is computed by a numerical quadrature, ${\int }_{S}f\left(\mathit{r}\right)dS={\sum }_{j=1}^{n_{qS}}f\left(\mathit{r}\right)w_j^S$ an that although ${\lambda }_{ch}^d$ depends on ${\lambda }_1^d$, ${\lambda }_2^d$ and ${\lambda }_3^d$, since ${\lambda }_1^d\cdot {\lambda }_2^d\cdot {\lambda }_3^d=1$, there are just 2 independent variables. Again a procedure similar to the one in the toy example, Section 3.1, is followed:

${\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill 0\hfill \end{array}\right]}_{j,X_{\text{ppal}}}={\left[\begin{array}{c}\hfill {\left[{\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_2^d\right)\right]}^T\hfill \\ \hfill {\left[{\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_2^d\right)\right]}^T\hfill \\ \hfill {\left[\mathbf{0}\right]}^T\hfill \end{array}\right]}_j\left[{\widehat{\mathit{P}}}_{ch}\right](35)$

$\begin{array}{ll}\hfill {\mathit{P}}_j:\nabla {\mathit{h}}_j^i& =\mathrm{s}\mathrm{y}\mathrm{m}\left({\mathit{h}}_j^T\nabla {\mathit{h}}_j^i\right):{\mathit{S}}_j\hfill \\ & ={\left[\begin{array}{ccc}\hfill l_{11}\hfill & \hfill l_{22}\hfill & \hfill 2l_{12}\hfill \end{array}\right]}_{j,X_{\text{ppal}}}^i{\left[\begin{array}{c}\hfill S_1\hfill \\ \hfill S_2\hfill \\ \hfill 0\hfill \end{array}\right]}_{j,X_{\text{ppal}}}\\ & ={\left(l_{11}\cdot {\mathit{S}}_{\text{row}}^1\left({\lambda }_1^d,{\lambda }_2^d\right)+l_{22}\cdot {\mathit{S}}_{\text{row}}^2\left({\lambda }_1^d,{\lambda }_2^d\right)\right)}_j^i\cdot {\widehat{\mathit{P}}}_{ch}\hfill \\ & ={\left({}^{*}\mathit{S}_{\text{row}}\right)}_j^i\left({\lambda }_1^d,{\lambda }_2^d\right)\cdot {\widehat{\mathit{P}}}_{ch}\end{array}(36)$

Looking at the previous matrix equation and using again (Eq. 13), it is straightforward to write the internal vector force for an element as

${\mathit{f}}_{\text{int}}^e={\mathrm{F}}_{\text{int}}^e\cdot {\widehat{\mathit{P}}}_{ch}(37)$

where ${\mathrm{F}}_{\text{int}}^e$ is a matrix. Now in the same way that the components that the components of ${\mathit{f}}_{\text{int}}^e$ can be assembled to obtain the global internal forces vector, f_int. The rows of ${\mathrm{F}}_{\text{int}}^e$ can be assembled to obtain a global internal force matrix ${\mathrm{F}}_{\text{int}}$.

${\mathit{f}}_{\text{int}}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{f}}_{\text{int}}^e=\left(\underset{e=1}{\overset{nel}{\bigwedge }}{\mathrm{F}}_{\text{int}}^e\right)\cdot {\widehat{\mathit{P}}}_{ch}={\mathrm{F}}_{\text{int}}\cdot {\widehat{\mathit{P}}}_{ch}(38)$

At that point it is important to note that if compared to the procedure followed in Section 3.1, now the internal force term is linked with the unknown vertices that correspond to the constituents behavior at the micro-scale. With that approach we show that scales can be crossed and information in one scale can be pushed up to other scales, that of course taking into account that the macro to micro connection has already been established.

From the overdetermined linear system of equations ${\mathrm{F}}_{\text{int}}\cdot {\widehat{\mathit{P}}}_{ch}={\mathit{f}}_{\text{ext}}$ a solution can be obtained for ${\widehat{\mathit{P}}}_{ch}$:

$\begin{array}{l}{\widehat{\mathit{P}}}_{ch}^k={\left({\mathit{A}}_{\text{sys}}^k\right)}^{-1}{\mathit{b}}_{\text{sys}};\\ {\mathit{A}}_{\text{sys}}^k={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathrm{F}}_{\text{int}}+{\mathit{D}}_2^T{\mathbf{\Omega }}_2{\mathit{D}}_2+{\mathit{D}}_1^T{\mathbf{\Omega }}_1^k{\mathit{D}}_1,\\ {\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathit{f}}_{\text{ext}}\end{array}(39)$

The discussion about the number of vertices and steps of load/displacement required in Section 3.1 is also applicable here. In this particular case, just the last step of deformation (u = 8) was used. Again, in case that additional information is needed to determine the function on the considered range, it is also possible to add more steps between u = 0 and u = 8. Since the initial solution was directly monotonically increasing, Ω₁ = 0 and it was not necessary to follow an iterative process, ${\widehat{\mathit{P}}}_{ch}$ was directly obtained from the simple initial system of equations:

${\widehat{\mathit{P}}}_{ch}={\left({\mathit{A}}_{\text{sys}}\right)}^{-1}{\mathit{b}}_{\text{sys}};{\mathit{A}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathrm{F}}_{\text{int}}+{\mathit{D}}_2^T{\mathbf{\Omega }}_2{\mathit{D}}_2,{\mathit{b}}_{\text{sys}}={\mathrm{F}}_{\text{int}}^T\mathit{W}{\mathit{f}}_{\text{ext}}(40)$

With the vertices obtained, the P-Spline could be reconstructed and compared it to the original one, see Figure 5.

FIGURE 5

FIGURE 5. Plot comparing $P_{ch}\left({\lambda }_{ch}^d\right)$ extracted from (Amores et al., 2020, Figure B.1) and the reverse-engineered $P_{ch}\left({\lambda }_{ch}^d\right)$ obtained with the methodology described in the article.

As it is shown in Figure 5, the $P_{ch}\left({\lambda }_{ch}^d\right)$ extracted from (Amores et al., 2020, Figure B.1) and the one obtained through the methodology described in the paper are practically identical. Looking closely, there are some small deviations close to ${\lambda }_{ch}^d=1$ on the compression range. As it has already been mentioned, the origin of those deviations is the eigenvector linked to an almost zero eigenvalue on the reduced matrix of the system, this matter has been carefully justified at the end of Section 3.1. With the obtained $P_{ch}\left({\lambda }_{ch}^d\right)$ the partial derivatives of the free energy can be determined for any state of deformation, with them, it is possible to re-run the FEM simulations and compare for example, the σ_VM plot. Figure 6 compares the σ_VM with $P_{ch}\left({\lambda }_{ch}^d\right)$ extracted from (Amores et al., 2020, Figure B.1) and the $P_{ch}\left({\lambda }_{ch}^d\right)$ with the methodology described in this paper. As we can see, although the functions are slightly different, that discrepancy has no effect on the final stresses, this is so because the additional spurious term that appears in the solution corresponds to a pressure term (in the cauchy stress sense) and therefore, it will not appear in the final computed stress due to the projector $\mathbb{P}$.

FIGURE 6

FIGURE 6. Plot comparing the results of σ_VM for the reference micro-mechanical model in Amores et al. (2020) that uses (Amores et al., 2020, Figure B.1) and the reverse-engineered one described in the article. (A) σ_VM using the $P_{ch}\left({\lambda }_{ch}\right)$ obtained for unfilled rubber, Treloar (1944) extracted from (Amores et al., 2020, Figure B.1). (B) σ_VM using the $P_{ch}\left({\lambda }_{ch}\right)$ reverse-engineered from the non-homogeneous virtual test using the methodology described in this article.

4 Conclusion

We have proposed a numerical method for determining both the continuum free energy and the polymer macromolecules behavior from arbitrary non-homogeneous DIC-based tests at the continuum scale. The procedure consists of a combination of our Macro-Micro-Macro approach and a Finite Element model. As a novel contribution of this approach, we show that by crossing scales transferring the microscale unknowns to the finite element formulation it is possible to determine the mechanical chain behavior from non-homogeneous experiments at the continuum scale by simply solving a linear system of equations (i.e., in an even more efficient manner than the subsequent simulations of the polymer behavior). Another key aspect is that Penalized B-splines (P-splines) preserve the general form of the energy function while retaining sufficient tools for enforcing specific desired conditions on the sought functions. The methodology at hand recovers the analytical free energies used as starting point in the virtual tests. Therefore, any finite element simulation performed with the reversed-engineered energy function will provide identical results to the original material model. From the authors perspective, this new approach opens new possibilities for data-driven characterization of the micro-structure from non-homogeneous tests at the macro-scale. Further research has to be conducted to evaluate the generalization of the approach in Amores et al. (2020) to more complex material behaviors from which the applicability of this same methodology has to be assessed.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

VA: Conceptualization, methodology, software, and writing original draft. FM: Conceptualization, methodology, funding acquisition, and supervision. EC: Methodology and supervision. FC: Conceptualization, methodology, funding acquisition, and supervision. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 101007815.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmats.2022.879614/full#supplementary-material

References

Alastrué, V., Martínez, M. A., Doblaré, M., and Menzel, A. (2009). Anisotropic Micro-sphere-based Finite Elasticity Applied to Blood Vessel Modelling. J. Mech. Phys. Sol. 57, 178–203. doi:10.1016/J.JMPS.2008.09.005

CrossRef Full Text | Google Scholar

Amores, V. J., Benítez, J. M., and Montáns, F. J. (2019). Average-chain Behavior of Isotropic Incompressible Polymers Obtained from Macroscopic Experimental Data. A Simple Structure-Based WYPiWYG Model in Julia Language. Adv. Eng. Softw. 130, 41–57. doi:10.1016/j.advengsoft.2019.01.004

CrossRef Full Text | Google Scholar

Amores, V. J., Benítez, J. M., and Montáns, F. J. (2020). Data-driven, Structure-Based Hyperelastic Manifolds: A Macro-Micro-Macro Approach to Reverse-Engineer the Chain Behavior and Perform Efficient Simulations of Polymers. Comput. Structures 231, 106209. doi:10.1016/J.COMPSTRUC.2020.106209

CrossRef Full Text | Google Scholar

Amores, V. J., Nguyen, K., and Montáns, F. J. (2021). On the Network Orientational Affinity assumption in Polymers and the Micro-macro Connection through the Chain Stretch. J. Mech. Phys. Sol. 148, 104279. doi:10.1016/J.JMPS.2020.104279

CrossRef Full Text | Google Scholar

Amores, V. (2022). PSPlines. doi:10.5281/zenodo.6084053

CrossRef Full Text | Google Scholar

Arruda, E. M., and Boyce, M. C. (1993). A Three-Dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials. J. Mech. Phys. Sol. 41, 389–412. doi:10.1016/0022-5096(93)90013-6

CrossRef Full Text | Google Scholar

Bergström, J. S., and Hayman, D. (2016). An Overview of Mechanical Properties and Material Modeling of Polylactide (PLA) for Medical Applications. Ann. Biomed. Eng. 44, 330–340. doi:10.1007/S10439-015-1455-8/TABLES/2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017). Julia: A Fresh Approach to Numerical Computing. SIAM Rev. 59, 65–98. doi:10.1137/141000671

CrossRef Full Text | Google Scholar

Carlsson, K., and Ekre, F. (2021). Ferrite.jl.

Google Scholar

Chagnon, G., Ohayon, J., Martiel, J.-L., Favier, D., Modeling, D. F. H., Chagnon, G., et al. (2017). “Hyperelasticity Modeling for Incompressible Passive Biological Tissues,” in Biomechanics of Living Organs. Academic Press (Elsevier), 3–30. doi:10.1016/B978-0-12-804009-6.00001-8

CrossRef Full Text | Google Scholar

Chagnon, G., Rebouah, M., and Favier, D. (2015). Hyperelastic Energy Densities for Soft Biological Tissues: A Review. J. Elast 120, 129–160. doi:10.1007/S10659-014-9508-Z

CrossRef Full Text | Google Scholar

Crespo, J., and Montáns, F. J. (2018). A Continuum Approach for the Large Strain Finite Element Analysis of Auxetic Materials. Int. J. Mech. Sci. 135, 441–457. doi:10.1016/J.IJMECSCI.2017.11.038

CrossRef Full Text | Google Scholar

Eggersmann, R., Kirchdoerfer, T., Reese, S., Stainier, L., and Ortiz, M. (2019). Model-Free Data-Driven Inelasticity. Comput. Methods Appl. Mech. Eng. 350, 81–99. doi:10.1016/J.CMA.2019.02.016

CrossRef Full Text | Google Scholar

Eilers, P. H. C., and Marx, B. D. (2021). Practical Smoothing: The Joys of P-Splines. Cambridge University Press. doi:10.1017/9781108610247

CrossRef Full Text | Google Scholar

Flaschel, M., Kumar, S., and De Lorenzis, L. (2021). Unsupervised Discovery of Interpretable Hyperelastic Constitutive Laws. Comput. Methods Appl. Mech. Eng. 381, 113852. doi:10.1016/J.CMA.2021.113852

CrossRef Full Text | Google Scholar

Flory, P. J., and Rehner, J. (1943). Statistical Mechanics of Cross‐Linked Polymer Networks II. Swelling. J. Chem. Phys. 11, 521–526. doi:10.1063/1.1723792

CrossRef Full Text | Google Scholar

Gent, A. N. (1989). Rubberlike Elasticity-A Molecular Primer, by James E. Mark and Burak Erman, John Wiley & Sons, New York, 1988, 196 Pp. J. Polym. Sci. C Polym. Lett. 27, 405–406. doi:10.1002/POL.1989.140271012

CrossRef Full Text | Google Scholar

González, D., Chinesta, F., and Cueto, E. (2018). Thermodynamically Consistent Data-Driven Computational Mechanics. Continuum Mech. Thermodyn. 31 (1), 239–253. doi:10.1007/S00161-018-0677-Z

CrossRef Full Text | Google Scholar

Holzapfel, G. A. (2002). Nonlinear Solid Mechanics: a Continuum Approach for Engineering Science. Meccanica 37, 489–490. doi:10.1023/a:1020843529530

CrossRef Full Text | Google Scholar

Ibañez, R., Abisset-Chavanne, E., Aguado, J. V., Gonzalez, D., Cueto, E., and Chinesta, F. (2018). A Manifold Learning Approach to Data-Driven Computational Elasticity and Inelasticity. Arch. Computat Methods Eng. 25, 47–57. doi:10.1007/S11831-016-9197-9/FIGURES/11

CrossRef Full Text | Google Scholar

Ibañez, R., Borzacchiello, D., Aguado, J. V., Abisset-Chavanne, E., Cueto, E., Ladevèze, P., et al. (2017). Data-driven Non-linear Elasticity: Constitutive Manifold Construction and Problem Discretization. Comput. Mech. 60, 813–826. doi:10.1007/S00466-017-1440-1

CrossRef Full Text | Google Scholar

Khiêm, V. N., and Itskov, M. (2016). Analytical Network-Averaging of the Tube Model:. J. Mech. Phys. Sol. 95, 254–269. doi:10.1016/J.JMPS.2016.05.030

CrossRef Full Text | Google Scholar

Kirchdoerfer, T., and Ortiz, M. (2016). Data-driven Computational Mechanics. Comput. Methods Appl. Mech. Eng. 304, 81–101. doi:10.1016/J.CMA.2016.02.001

CrossRef Full Text | Google Scholar

Kirchdoerfer, T., and Ortiz, M. (2018). Data-driven Computing in Dynamics. Int. J. Numer. Meth. Engng 113, 1697–1710. doi:10.1002/NME.5716

CrossRef Full Text | Google Scholar

Korzeniowski, T. F., and Weinberg, K. (2021). A Multi-Level Method for Data-Driven Finite Element Computations. Comput. Methods Appl. Mech. Eng. 379, 113740. doi:10.1016/J.CMA.2021.113740

CrossRef Full Text | Google Scholar

Latorre, M., and Montáns, F. J. (2013). Extension of the Sussman-Bathe Spline-Based Hyperelastic Model to Incompressible Transversely Isotropic Materials. Comput. Structures 122, 13–26. doi:10.1016/J.COMPSTRUC.2013.01.018

CrossRef Full Text | Google Scholar

Latorre, M., and Montáns, F. J. (2014). What-You-Prescribe-Is-What-You-Get Orthotropic Hyperelasticity. Comput. Mech. 53, 1279–1298. doi:10.1007/S00466-013-0971-3

CrossRef Full Text | Google Scholar

Liu, X., Tao, F., Du, H., Yu, W., and Xu, K. (2020). Learning Nonlinear Constitutive Laws Using Neural Network Models Based on Indirectly Measurable Data. J. Appl. Mech. Trans. ASME 87. doi:10.1115/1.4047036/1083320

CrossRef Full Text | Google Scholar

Marckmann, G., and Verron, E. (2006). Comparison of Hyperelastic Models for Rubber-like Materials. Rubber Chem. Techn. 79, 835–858. doi:10.5254/1.3547969

CrossRef Full Text | Google Scholar

Mooney, M. (1940). A Theory of Large Elastic Deformation. J. Appl. Phys. 11, 582–592. doi:10.1063/1.1712836

CrossRef Full Text | Google Scholar

Moreno, S., Amores, V. J., Benítez, J. M., and Montáns, F. J. (2020). Reverse-engineering and Modeling the 3D Passive and Active Responses of Skeletal Muscle Using a Data-Driven, Non-parametric, Spline-Based Procedure. J. Mech. Behav. Biomed. Mater. 110, 103877. doi:10.1016/J.JMBBM.2020.103877

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen-Thanh, V. M., Zhuang, X., and Rabczuk, T. (2020). A Deep Energy Method for Finite Deformation Hyperelasticity. Eur. J. Mech. - A/Solids 80, 103874. doi:10.1016/J.EUROMECHSOL.2019.103874

CrossRef Full Text | Google Scholar

Ogden, R. W. (1972). Large Deformation Isotropic Elasticity - on the Correlation of Theory and experiment for Incompressible Rubberlike Solids. Proc. R. Soc. Lond. A. 326, 565–584. doi:10.1098/RSPA.1972.0026

CrossRef Full Text | Google Scholar

Sáez, P., Alastrué, V., Peña, E., Doblaré, M., and Martínez, M. A. (2011). Anisotropic Microsphere-Based Approach to Damage in Soft Fibered Tissue. Biomech. Model. Mechanobiol 11, 595–608. doi:10.1007/S10237-011-0336-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sussman, T., and Bathe, K.-J. (2009). A Model of Incompressible Isotropic Hyperelastic Material Behavior Using Spline Interpolations of Tension-Compression Test Data. Commun. Numer. Meth. Engng. 25, 53–63. doi:10.1002/CNM.1105

PubMed Abstract | CrossRef Full Text | Google Scholar

Treloar, L. R. G. (1944). Stress-Strain Data for Vulcanized Rubber under Various Types of Deformation. Rubber Chem. Techn. 17, 813–825. doi:10.5254/1.3546701

CrossRef Full Text | Google Scholar

Treloar, L. R. G. (1975). The Physics of Rubber Elasicity. OUP Oxford, 322.

Google Scholar

Volokh, K. (2016). Mechanics of Soft Materials. Singapore: Springer, 337.

Google Scholar

Zheng, Z., Ye, H., Zhang, H., Zheng, Y., and Chen, Z. (2020). Multi-Level K-D Tree-Based Data-Driven Computational Method for the Dynamic Analysis of Multi-Material Structures. Int. J. Mult. Comp. Eng. 18, 421–438. doi:10.1615/INTJMULTCOMPENG.2020035167

CrossRef Full Text | Google Scholar

Glossary

Ω₁ Weight matrix for the first differences

Ω₂ Weight matrix for the second differences

A Green Lagrange strain tensor

b Volumetric forces

C Right Cauchy Green tensor

D₁ First differences matrix

D₂ Second differences matrix

${\mathit{f}}_{\text{ext}}^e$ Element external nodal forces

${\mathit{f}}_{\text{int}}^e$ Element internal nodal forces

f_ext External nodal forces

f_int Internal nodal forces

N_i Principal referential direction of deformation number i

P PK1 (First Piola-Kirchhoff) stress tensor

r Arbitrary direction of the unit sphere in the reference configuration

S PK2 (Second Piola-Kirchhoff) stress tensor

t Surface forces

u Displacement field

u^a Displacement vector of node a

W Weight matrix for the system of equation

X Gradient of deformation tensor, X = ∂^tx/∂⁰x

∇ Nabla operator, ∇ = ∂_ie_i

δu Virtual displacement field

δW Virtual work

δW_ext External virtual work

δW_int Internal virtual work

Γ Area in the reference configuraion

κ Bulk modulus

$\mathbb{P}$ Fourth order projector tensor $\mathbb{P}={\mathbb{I}}^S-1/3{\mathit{C}}^{-1}\otimes \mathit{C}$

⁰ξ Local coordinates for the unit reference element

μ Shear modulus

$\mathrm{\Omega }$ volume in the reference configuraion

$\mathrm{\Omega }$ volume in the reference configuraion

$\mathrm{\Psi }$ Continuum free energy function

ψ_ch Chain free energy function

□ Volume of the unit reference element

$B_i\left(x\right)$ B-Spline basis function corresponding to the vertex i

h^a Shape function of node a in the reference unit element

$I_1^C$ First principal invariant of C, $I_1^C={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2\lambda$

J Volume ratio $J=\det \left(X\right)$

$J_j^e$ Jacobian of the element e at integration point j

MSE Mean square error

n_dofs Number of DOFs per element

n_vert Number of vertices of the BSpline

n_el Number of elements

n_n Number of nodes

n_qp Number of quadrature points per element

n_qS Number of quadrature points on the microspheres

P_ch Chain mechanical behavior function, P_ch = dψ_ch/dλ_ch

r_i Director cosine of the vector r with respect to N_i

S_i Principal value i of the PK2 stress tensor

w_j Quadrature weight for the integration point j of the element

${\mathit{w}}_{\mathit{j}}^{\mathit{S}}$ Quadrature weight for the integration point j on the sphere

${\mathit{\lambda }}_{\mathit{c}\mathit{h}}^{\mathit{C}}$ Affine chain stretch, ${\lambda }_{ch}^C=\sqrt{\mathit{r}\cdot \mathit{C}\cdot r}$

${\mathit{\lambda }}_{\mathit{i}}^{\mathit{d}}$ Deviatoric principal stretch i

λ_ch Non-affine chain stretch, λ_ch = r ⋅U ⋅ r

λ_i Principal stretch i

$\mathcal{U}$ Volumetric contribution to the free energy function

$\mathcal{W}$ Deviatoric contribution to the free energy function