Data-Driven GENERIC Modeling of Poroviscoelastic Materials

Biphasic soft materials are challenging to model by nature. Ongoing efforts are targeting their effective modeling and simulation. This work uses experimental atomic force nanoindentation of thick hydrogels to identify the indentation forces are a function of the indentation depth. Later on, the atomic force microscopy results are used in a GENERIC general equation for non-equilibrium reversible–irreversible coupling (GENERIC) formalism to identify the best model conserving basic thermodynamic laws. The data-driven GENERIC analysis identifies the material behavior with high fidelity for both data fitting and prediction.


Introduction
Thanks to the recent progress in simulation technology and computing power, the mechanical behavior of biological tissues is nowadays one of the most active research topics. However, many biological tissues are biphasic by nature, which renders their effective modeling and simulation a challenging issue, even with the impressive progress achieved recently [1]. Human cartilage, for example, is a biphasic material, where the fluid pressurization is believed to be the main load-carrying phenomenon [2]. Moreover, the microstructure and the fluid-solid interactions in such materials is complex and not fully understood [3,4]. The modeling of biphasic materials is, however, mandatory to design effective replacements of human soft tissues as well as understanding their behavior [5].
Different efforts are undertaken to model soft materials. In [5], the authors attempt to define a viscoelastic model based on only three parameters identified experimentally, using a pure solid mechanics approach. Other works aim to model the contact and lubrication phenomenon when using soft materials in contact mechanics, similarly to human body lubricated contact [6]. In [6], an interface element is developed and identified to transmit contact efforts from contact bodies into the lubrication fluid, soft material film. A study of cutting behavior of soft material films is performed in [7] using energy methods, however, assuming a non-linear hyperelastic neoprene rubber-like material. In [1], the authors model the soft material indentation using a biphasic approach, combining an elastic solid behavior and a fluid pressurization one. Indentation of soft materials is also studied in [8], where the authors performed a theoretical and finite element study using Hertzian contact and generalizing for large hyperelastic deformation including material non-linearities. Suitable nonlinear compressible material model for soft materials is also derived in [9], with its parameters identified for lung parenchyma using experimental results and inverse analysis. and then the GENERIC formalism used in this work is detailed along with the identified model. Finally, the numerical results are illustrated and discussed.

The Experimental Procedure
In this section, we review the experimental setup used to identify the mechanical behavior of a thick hydrogel. The experiment uses atomic force microscopy nanoindentation (AFM nanoindentation) to identify the mechanical response of a hydrogel. AFM is usually modeled as a cantilever beam supporting a semi-spherical indenter of radius ρ = 36 µm. In turn, the cantilever beam is modeled as a spring of stiffness k = 2.88 N/m, as illustrated in Figure 1.  The specimen indentation depth is w(t), while the displacement at the base is z(t). The experimental setup identifies z(t) using optical sensors and the reaction force F at the base of the spring, which corresponds to the force in the spring. Therefore one may easily write: The experiment is repeated for five different indentation rates,ż, equal to 0.5 µm/s, 2 µm/s, 8 µm/s, 40 µm/s and 80 µm/s. The experimental results are illustrated in Figure 2, which provides real experimental data obtained using atomic force microscopy (AFM) non-indentation. The experiments were performed at room temperature or 22 • C ± 1 • C. The normal spring constant of the cantilever was measured using the thermal noise method before attaching the colloidal microsphere indenter. The indenter is a silica microsphere glued with UV-curable glue (Norland optical adhesive 63) to the end of the tipless cantilever by means of a home-built micromanipulator. The exact approach rateṡ z(t) were measured using a Z-piezo sensor to control the approach of the probe as much as possible. The experimental setup measures z and F at every time step ∆t = 0.5 ms. One may note that the indentation depth does not exceed 1.8 µm as shown in Figure 2, and therefore we can approximate the spherical indenter by a flat one considering the radius of curvature of the indenter ρ as very large with respect to the indentation depth. Thus the indentation area can be approximated as A = πR 2 with R the indentation radius defined as: Indentation force F (N) Moreover, we note that the initial hydrogel specimen height H 0 is large enough to avoid any substrate effect during the indentation. For instance, in the tested specimen, H 0 = 6 mm, more than 1000 times the maximum reached indentation depth.

GENERIC Formalism
Under the framework of non-equilibrium thermodynamics, the GENERIC formalism establishes a completely general equation of the dynamics of a system under reversible and irreversible conditions, as dictated by the evolution of energy and entropy, respectively [22,23,27]. Models constructed using the GENERIC formalism preserve the symmetries of the system and therefore guarantee the conservation of energy and the increase of entropy.
The GENERIC structure of the evolution equations for an arbitrary problem iṡ where: 1.
[z t ] is the vector of state variables of the problem at time t. The choice of z t is irrelevant in the sense that different sets of z t lead to different GENERIC formalisms, all of them thermodynamically consistent. However, z t should contain variables-such as position, momentum, stress, energy, or entropy, for instance-able enough to evaluate the energy conservation and the dissipation therms. L is the so-called Poisson matrix and will be responsible for the reversible (Hamiltonian) part of the evolution of the system. E represents the energy of the system, as a function of its particular state at time t, z t . M represents the friction matrix, responsible for the irreversible part of the evolution of the system. S represents the entropy of the system for the particular choice of variables z.
Equation (3) is supplemented with the complementary degeneracy conditions, i.e., In what follows, if there is no risk of confusion, and for the sake of readability, we will omit the subscript z in ∇ z .
By choosing L skew-symmetric and M symmetric, positive semi-definite, one ensures the conservation of energy: and the fulfillment of the second principle of thermodynamics:

Data-Driven Characterization of the GENERIC Description of a Hydrogel
The objective is to model the biphasic hydrogels using GENERIC formalisms, thus to obtain the GENERIC Equation (3) for the experimental data. In [18,28] an approach has been developed by the authors so as to obtain a numerical description-by means of manifold learning techniques-of the different constituents of the GENERIC equation.
Since data z i are obtained at discrete time increments, Equation (3) is first discretized in time, where, for simplicity, we denote z n+1 = z t+∆t and where L and M are the discrete version of the Poisson and friction operators, respectively. Also, DE and DS represent the discrete gradients. In general, matrix L is constant over the process, while matrix M is frequently a function of z.
The set of variables chosen for each z i are the specimen indentation depth w(t), the specimen indentation velocity v(t) =ẇ(t) and the average normal stress at the indentation point σ(t) at time t i = i × ∆t . Also, in this paper, we assume L known, being: As can be seen in [28] , different definitions of z, L and M can be done, and also the particular structure of L can be considered itself as an unknown. The ones selected here have shown the best convergence for this particular problem. Furthermore, the discrete gradients are discretized in a piecewise linear, finite element, manner: where A and B can be time dependent operators or constant with respect to the time. For the selected problem, operators were shown to be almost constant. Therefore, we decided to perform the regression procedure within a single step. The proposed algorithm will thus consist in solving the following minimization problem: subjected to: with z(µ) given by Equation (7) and z exp all the experimental results of each experiment.
The just-introduced procedure must not be seen as a model fitting procedure. Indeed, the number of values to determine is much higher than the usual number of parameters in a suitable model. What we seek with this procedure is a method of machine learning the GENERIC expression of the problem. The result is the numerical value of these GENERIC building blocks and not the particular values of any model parameter. Those are considered constant for each experiment, although a piece-wise linear variation along time is equally possible, in general.
Once the minimization problem has been solved, and given an initial z 0 value of the variable at the beginning of the experiment, it is possible to reconstruct the solution for each experiment, in order to check the accuracy of the method:

Equivalence With Traditional Ways of Phenomenological Model Fitting
For practitioners used to employ experimental results for constitutive law fitting purposes, the procedure outlined above could seem intricate and unclear. Particularly, regression of the terms in Equation (3) could somehow obscure the true philosophy behind the suggested method.
Note, however, that some terms in Equation (3) should seem familiar to us. Noteworthy, the term ∇E(z) represents the usual form of deriving a constitutive law in a hyperelastic framework. In turn, the term ∇S(z) introduces a second potential, of dissipative nature, taking care of the viscous terms in the constitutive law.
From a practical point of view, the development of finite element time integration schemes deriving directly from GENERIC expressions is something already developed by different authors, particularly I. Romero and coworkers. Equation (3) could be employed advantageously to derive time integration schemes directly that conserve energy and dissipate entropy, showing great robustness from the numerical point of view, even for low orders. The interested reader is referred to [29][30][31] and references therein for details.

GENERIC Model
The experimental measurements are taken every ∆t = 0.015 seconds, displacements are measured in micrometers µm, velocities in micrometers per second µm/s and stresses in nanoNewtons (nN) per square micrometer nN/(µm) 2 . A regression for A, M, and B is performed using different sets of measurements at different indentation rates. The experiment is performed on two different sets of materials. Figure 3a,b show the indentation depth w as a function of the number of measurements taken from the beginning of the experiment, which is proportional to the time t since measurements are taken every 15 ms. We can clearly see an excellent comparison of the derived GENERIC model with the experimental indentation depth, with an average relative error of 0.4%.     Using the generic formalism, we can also show the indenter velocityẇ obtained from the model and the experimental data through time derivation of the obtained data points w using explicit Euler's derivation. Figure 5a,b show, for both materials, the experimental values of the indenter velocityẇ for differentż, along with their corresponding identified GENERIC models, obtained with the GENERIC regressions. In Figure 5, we can clearly see experimental artifacts of the velocityẇ at high indentation ratesż and near the end of the experiment, for both materials. This is mainly due to the low machine precision at high speeds and potentially induced vibration of the cantilever arm. The relative errors for different GENERIC formalism modeling, for the three quantities of interests, are illustrated in Figures 6 and 7 for the MICA silicate and polyethylene fibers hydrogels respectively.  Indentation velocity error (%)     Stress error (%) All experiments have low relative errors in general. We can clearly see negligible relative error for the indentation depth w for both materials in Figures 6a and 7a. We can also identify very low relative errors on the indenter's velocityẇ in Figures 6b and 7b, except by the end of the experiments atż = 40 µm/s andż = 80 µm/s, which is explained by the experimental artifacts seen by the end of the experiments. The stresses exhibit the highest values of the relative error while being in acceptable ranges eventually. The high relative errors at the beginning of the experiments illustrated in Figures 6c and 7c are explained by the numerical and experimental amplified errors by the derivation of the results, since computing the velocities requires a numerical derivation of the experimental measurements.
Tables 1 and 2 illustrate the average relative errors for the MICA silicate and polyethylene indentations experiments respectively, calculated using: where N is the number of measurements in the experiment using a given indentation rate.

Predictive Capabilities for New Experimental Results
In this section we use the aforementioned GENERIC regression to predict the experimental results for an indentation velocityż = 2 µm/s, while not being present in the training database. For that aim, we will use the GENERIC formalism matrices A, M and B obtained atż = 0.5 µm/s andż = 8 µm/s. Using these results, we predict the matrices atż = 2 µm/s using a linear interpolation approach. For instance: Using the linear interpolation approach, we predict the GENERIC formalism matrices. Eventually, using a more accurate interpolation scheme is an option for interested users (Kriging for example). The prediction results are illustrated in Figures 8-10. The results show good predictions bounded between the two fitted quantities of interest atż = 0.5 µm/s andż = 8 µm/s illustrated in large dashed lines in Figures 8-10.

Discussion
This work consists of a first attempt to explore the possibilities of using the GENERIC formalism with data-driven identification, to model the indentation of thick hydrogels. The results show a good match with the experimental data. The use of three different matrices to identify in the model gives an advantage and flexibility of identifying a model with 27 parameters. The 27 parameters (3 matrices each containing 9 values) are eventually enough to reproduce the best behavior of the material. Indeed, not all the identified parameters are independent. For instance, the symmetry of the M matrix and the skew-symmetry of the L matrix reduces the number of independent variables to 18. Other restrictions may reduce even further the number of independent variables. Classical modeling of the indentation is possible, either using solid mechanics approach [32], or a combination of solid and fluid phases [1].
Different matrices A, M, and B are identified for different indentation velocities. This fact shows a highly non-linear behavior of the thick hydrogel in question. To simulate the behavior of the thick hydrogel at an indentation rate different than the studied ones, an SSL-PGD interpolation, for example, can be performed, for each one of the 27 identified parameters, among many nonlinear interpolation techniques [33]. The result of the identified matrices is illustrated in Appendices A and B for the two selected hydrogels.
The choice of using the GENERIC formalism appears to be a suitable approach for identifying models with challenging mechanical behavior but without prior exact knowledge of the constitutive equations. It also shows good predictive abilities for unfitted experimental results, as illustrated in Section 3.2.

Conclusions
In this work, we investigate the possibility of using the GENERIC formalism for modeling non-trivial material behavior. For instance, the work focuses on modeling thick hydrogels using only the conservation of thermodynamic laws. The results show a good correlation with the experimental ones for different experiments. The predictive ability of the model is illustrated in Section 3.2, where linear interpolation of the model matrices is shown to suffice for predicting new experimental models. Extrapolation of the results requires further thorough investigation and developments to yield good results. Funding: This project has been partially funded by the ESI Group through the ESI Chair at ENSAM Arts et Metiers Institute of Technology, and through the project "Simulated Reality" at the University of Zaragoza. The support of the Spanish Ministry of Economy and Competitiveness through grant number CICYT-DPI2017-85139-C2-1-R and by the Regional Government of Aragon and the European Social Fund, are also gratefully acknowledged.

Acknowledgments:
The authors would like to acknowledge the useful comments and help provided by R. Simik and C. Mathis regarding the experimental results.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Polyethylene Hydrogels Matrices
In this appendix we illustrate the final identified matrices by the considered algorithm for the Polyethylene hydrogels: