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, δ W = δ u ⋅ f int − f ext = 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 δ W = δ u ⋅ f int − f ext = 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: ∫ Ω δ u ⊗ ∇ : P u d Ω = ∫ Ω δ u ⋅ b d Ω + ∫ Γ δ u ⋅ t d Γ (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 δ u ⊗ ∇ : P u = P : ∇ δ u is the internal virtual work density.

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.

A typical assumption made on the phenomenological approach for isotropic incompressible solids is the Valanis Landel decomposition, W λ 1 d , λ 2 d , λ 3 d = ω λ 1 d + ω λ 2 d + ω λ 3 d . With that assumption, the expressions for the principal PK2 stresses can be obtained in terms of an unknown function ω ′ x : S j = J − 2 / 3 1 λ j d ω ′ λ j d − λ 3 d λ j d ω ′ λ 3 d f o r j = 1,2 (20)

Once ω ′ x 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: ω ′ x = ∑ i = 1 n vert B i x ω ̂ i = B 1 x … B n vert x ω ̂ 1 ⋮ ω ̂ n vert = B r o w x ⋅ ω ̂ ′ (21)

With the previous expression, the principal stresses could be written as: S j = J − 2 / 3 1 λ j d B r o w λ j d − λ 3 d λ j d B r o w λ 3 d ⋅ ω ̂ ′ = S row j λ j d , λ 3 d ⋅ ω ̂ ′ f o r j = 1,2 (22) S 1 S 2 0 j , X ppal = S row 1 λ 1 d , λ 3 d T S row 2 λ 2 d , λ 3 d T 0 T j ω ̂ ′ T (23) P j : ∇ h j i = s y m X j T ∇ h j i : S j = l 11 l 22 2 l 12 j , X ppal i S 1 S 2 0 j , X ppal = l 11 ⋅ S row 1 λ 1 d , λ 3 d + l 22 ⋅ S row 2 λ 2 d , λ 3 d j i ⋅ ω ̂ ′ = S row * j i λ 1 d , λ 2 d , λ 3 d ⋅ ω ̂ ′ (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: f int e = F int e ⋅ ω ̂ ′ (25)

Doing the assembly of the internal nodal forces of the elements in the mesh: f int = ⋀ e = 1 n e l f int e = ⋀ e = 1 n e l F int e ⋅ ω ̂ ′ = F int ⋅ ω ̂ ′ (26)

From the overdetermined linear system of equations F int ⋅ ω ̂ ′ = f ext we can solve for ω ̂ ′ : ω ̂ ′ k = A sys k − 1 b sys ; A sys k = F int T W F int + D 2 T Ω 2 D 2 + D 1 T Ω 1 k D 1 , b sys = F int T W f ext (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, ω ̂ ′ was directly obtained from the simple initial system of equations: ω ̂ ′ = A sys − 1 b sys ; A sys = F int T W F int + D 2 T Ω 2 D 2 , b sys = F int T W f ext (28)

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

Since the function ω ′ x 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.

In this particular example, the matrix of the system A sys ⋅ ω ̂ ′ = b sys , presents a null eigenvalue with its corresponding associated eigenvector, i.e., there is a certain subspace of possible solutions ω ̂ 0 ′ (containing just the direction of the eigenvector associated to the null eigenvalue) such that added to the solution, ω ̂ ′ , produce null effect on the independent term, b _sys and therefore ω ̂ ′ + α ω ̂ 0 ′ , 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 ω ′ λ i d − λ 3 d / λ i d ⋅ ω ′ λ 3 d = 0 ∀ λ 1 d , λ 2 d , λ 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) ω ′ λ i d λ i d = ω ′ λ 3 d λ 3 d ∀ λ 1 d , λ 2 d , λ 3 d . The only way two functions of different independent variables can be equal ∀ λ 1 d , λ 2 d , λ 3 d is that those functions are constant ω ′ λ i d λ i d = ω ′ λ 3 d λ 3 d = β 0 → ω ′ λ = β 0 / λ . Since the shape of the eigenvector associated to the null eigenvalue exactly matches the shape of the function, ω ′ λ = β 0 / λ , 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 λ 1 , λ 2 , λ 3 / λ i 2 , the same form of the term resulting from introducing an additional β 0 / λ i d to the function ω ′ λ i d . 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: S = S v + S d ; S d = 2 d W d C = 2 J − 2 / 3 P : d W d C d = J − 2 / 3 P : S ∣ d (29) where C =  X ^T X is the right Cauchy-Green tensor and P is the fourth order projector tensor, see Supplementary Appendix SB. Since P : α C − 1 = 0 , from (Eq. 29) it can be seen that there are multiple S ∣ d = ∂ W / ∂ C 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 P .

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 c h λ c h = d ψ c h / d λ c h (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 c h λ c h 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): W λ 1 , λ 2 , λ 3 = 1 S ∫ S ψ c h λ c h d S (30)

For FEM simulations even pure incompressible solids are simulated with the quasi-incompressible material framework: Ψ A = U J + W λ 1 d , λ 2 d , λ 3 d in which U J is a penalization function that leads to incompressibility for κ → ∞ like U J = 1 2 k J − 1 2 and where W λ 1 d , λ 2 d , λ 3 d is the same as the one used in pure incompressibility but with the deviatoric stretches instead of the total ones: W λ 1 d , λ 2 d , λ 3 d = 1 S ∫ S ψ c h λ c h d d S ; λ c h d = U d : r ⊗ r = λ 1 d r 1 2 + λ 2 d r 2 2 + λ 3 d r 3 2 (31) W k = ∂ W ∂ λ k d = 1 S ∫ S d ψ c h λ c h d d λ c h d d λ c h d d λ k d d S = 1 S ∫ S P c h λ c h d r k 2 d S (32)

The unknown function that will be approximated using P-Splines is P c h λ c h d , 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: P c h x = ∑ i = 1 n vert B i x P ̂ c h i = B 1 x … B n vert x P ̂ c h i ⋮ P ̂ c h n vert = B r o w x ⋅ P ̂ c h (33)

Introducing the expression for P c h x in (Eq. 32) and then all the derivatives of the form 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: S j = J − 2 / 3 1 λ j d 1 S ∫ S B row λ c h d r j 2 − λ 3 d λ j d r 3 2 d S ⋅ P ̂ c h = S row j λ 1 d , λ 2 d ⋅ P ̂ c h f o r j = 1,2 (34)

Note that the previous integral in the microsphere is computed by a numerical quadrature, ∫ S f r d S = ∑ j = 1 n q S f r w j S an that although λ c h d depends on λ 1 d , λ 2 d and λ 3 d , since λ 1 d ⋅ λ 2 d ⋅ λ 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: S 1 S 2 0 j , X ppal = S row 1 λ 1 d , λ 2 d T S row 1 λ 1 d , λ 2 d T 0 T j P ̂ c h (35) P j : ∇ h j i = s y m h j T ∇ h j i : S j = l 11 l 22 2 l 12 j , X ppal i S 1 S 2 0 j , X ppal = l 11 ⋅ S row 1 λ 1 d , λ 2 d + l 22 ⋅ S row 2 λ 1 d , λ 2 d j i ⋅ P ̂ c h = S row * j i λ 1 d , λ 2 d ⋅ P ̂ c h (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 f int e = F int e ⋅ P ̂ c h (37) where F int e is a matrix. Now in the same way that the components that the components of f int e can be assembled to obtain the global internal forces vector, f _int. The rows of F int e can be assembled to obtain a global internal force matrix F int . f int = ⋀ e = 1 n e l f int e = ⋀ e = 1 n e l F int e ⋅ P ̂ c h = F int ⋅ P ̂ c h (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 F int ⋅ P ̂ c h = f ext a solution can be obtained for P ̂ c h : P ̂ c h k = A sys k − 1 b sys ; A sys k = F int T W F int + D 2 T Ω 2 D 2 + D 1 T Ω 1 k D 1 , b sys = F int T W f ext (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, P ̂ c h was directly obtained from the simple initial system of equations: P ̂ c h = A sys − 1 b sys ; A sys = F int T W F int + D 2 T Ω 2 D 2 , b sys = F int T W f ext (40)

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

As it is shown in Figure 5, the P c h λ c h d 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 λ c h 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 c h λ c h d 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 c h λ c h d extracted from (Amores et al., 2020, Figure B.1) and the P c h λ c h d 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 P .