Comptes Rendus Mécanique

A hybrid modeling combining the proper generalized decomposition approach to data-driven model learners, with application to nonlinear biphasic materials


Introduction
Data-driven techniques are becoming increasingly popular in many fields, engineering is not an exception.Although science always relied on data, it is only until the beginning of the 21st century that the eruption of data-driven modeling took place [1].This phenomenon was motivated by both the increase of data availability and sharing, as well as the increase in computing power.Many algorithms are developed based on linear or nonlinear regressions, decision trees, random forests, neural networks, and so on.All these models tend to uncover a relation between an input and an output while neglecting the physical phenomenon involved.As a result, the modeling interpretation becomes too complicated if not possible in many situations.
On the other hand, the development of simulation techniques for solving partial differential equations (PDEs) made a tremendous progress in the filed of developing real-time simulations on the fly, with minimal calculation time, thanks to the domain decomposition techniques [2][3][4].Moreover, many model reduction techniques perform an a posteriori model reduction by compressing the available data using the most relevant information [5][6][7][8].The available techniques allow, after a pre-possessing offline calculation, to use the computed parametric solutions in realtime in what is called the online phase for real-time optimization and control [9,10].The model reduction techniques were successfully used in many optimization problems since the gradients with respect to the process parameters are readily available in the treated problems [4,11].
Using model reduction techniques or reduced order basis methods leads in general to a good approximation, which tends to the finite element approximation when increasing the size of the projection basis [12][13][14].However, model reduction involves truncation errors in general, on top of the modeling error.Therefore, a room for improving model reduction techniques is possible.In this work, we attempt to use advanced machine learning algorithms to correct modeling and simulation errors.Moreover, and despite the impressive progress in computational power, many problems remain intractable using the classical model reduction method [15] or remains too intrusive to be included in an existing commercial software [16].Despite numerous developments in both intrusive and nonintrusive manners [17], the calculation of solutions using model reduction techniques faces major challenges when changing a parameter induces drastic changes in the solution map [18].
With the recent development of artificial intelligence and machine learning algorithms, some combinations of model reduction techniques with machine learning are addressed to account for error reduction and reliability improvements, either by error correction [19] or by improvement of modal selection using the quantities of interest [20].
In this work, we aim to combine the physical models generated in a reduced basis with machine learning algorithms to create a physically informed machine learning algorithms.Eventually, with the multiparametric optimization involved in any machine learning technique involving many iterations and derivatives, computations, any classical simulation algorithm will fail to create a physically informed neural network due to the multiple simulations required at each iteration.Therefore, we leverage the availability of a parametric simulation using the proper generalized decomposition (PGD) model reduction technique to learn the physical problem parameters along with the machine learning algorithms hyper parameters.The selected application for this work is a nonlinear viscoelastic biphasic material.In fact, 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, often biological tissues are biphasic by nature, which renders their effective modeling and simulation a challenging issue, even with the impressive progress achieved recently.Human cartilage, for example, is a biphasic material, where the fluid pressurization is believed to be the main load carrying the phenomenon [17,21].The microstructure and the fluid-solid interactions in such materials is still not fully understood [22,23].Previous works attempted to model soft permeable matter as a superposition of a fluid force contribution to an elastic solid one [17].The previous modeling does not take into consideration the possibility of having a nonelastic behavior of the solid counterpart of the biphasic material.In this work, we apply the previous model to the atomic force microscopy nanoindentation of a biphasic material with a large thickness when compared to the indentation depth.
First of all, we will revisit the biphasic soft material modeling in Section 2. Later on, we use the PGD to simulate the indenter behavior and explain how the obtained reduced order model results can lead to an improvement of the model in Section 3. Once the PGD results  are available, we leverage a multidimensional parametric solution and the experimental results shown in Section 4 to identify the permeability of a biphasic hydrogel using an elastic behavior of the solid component of the hydrogel in question, in Section 5.1.The results being none fully satisfactory, a correction of the model is performed twice, once using an analytical modeling, and next time using a fully coupled, physically informed, machine learning algorithm, in Sections 5.2 and 5.3, respectively.It is finally shown that the modeling coupling the simulation and artificial intelligence neural network performs better than both other models.

Modeling of AFM nanoindentation of soft matter
Classically, the atomic force microscope (AFM) nanoindentation of a poroelastic brush is performed using a rigid spherical probe [17,24], rigidly attached to a cantilever, as illustrated in Figure 1.The indentation process is controlled by the motion of the base of the cantilever, z(t ), as shown in Figure 1.An optical sensor measures z(t ), the displacement of the cantilever, right above the colloidal probe (the indenter).In this work, the cantilever is modeled as a linear spring with a constant resulting stiffness k s = 2.88 N/m.An increasing force F (t ) is applied on the spring to indent the specimen.This force is measured during the indentation.
The indentation domain is cylindrical of initial height H 0 , the indenter is a semi-spherical probe of radius ρ = 12 µm.The equilibrium equation is written as F fluid being the fluid contribution to the reaction force on the indenter, overcoming the spring compression force, while F solid is the solid phase contribution to the reaction, and F spring is the measured force F (t ) imposed on the specimen.Knowing k s , the spring force is written as where w(t ) represents the indentation depth into the specimen.The reference of z(t ) = 0 is considered as the position of the spring base at the initial contact of the indentation probe with the specimen.Note that the force of the spring is time-dependent during the penetration phase.In the following, for the sake of notation simplicity, the force in the spring is simply denoted F (t ).Next, we model the individual contributions of the fluid and the solid in the biphasic material.

The fluid contribution to the reaction force
The pressure in the fluid can be found using Darcy's law, modeling the fluid behavior in a porous media.One can easily demonstrate by a simple projection that the integral of the pressure over a spherical probe is equal to the integral of the same applied pressure over a disk of the same radius R(t ).If ρ is the radius of the indenter, the radius of the contact disk is written as [17] (3) R(t ) changes as a function of time since the probe penetration depth w(t ) is changing.Using Darcy's law we may write v being the fluid velocity, µ the fluid viscosity, and K the effective permeability, which also takes into account the material porosity.In this work, we are using a highly porous media, with a fluid percentage above 90%, and therefore the conservation of mass is written as [17] ∇ Replacing ( 4) into ( 5) and assuming an homogeneous permeability K and fluid viscosity µ in the domain, we find Darcy's pressure equation to solve along with the boundary conditions specified in Section 3. One may note that the permeability of the domain K, the radius R, the height of the domain H , and the compression velocity U are time-dependent.The indentation velocity can be determined by deriving the indentation depth w(t ) at a given generic time step n where ∆t is the time step.Moreover, one may define the actual domain height H (t ) under the indenter at a generic time step by During indentation, the height H n under the indenter is changing, as well as the indenter radius R n of contact with the specimen and the indentation velocity U n .While the permeability K is an unknown parameter to identify, we will need to solve the problem for different values of H n , U n , R n the radius at a time step n and permeability K. Therefore, having a parametric solution of the problem P (x, y, z, H ,U , R, K) becomes an appealing approach.

Solid contribution to the reaction force
In this section, we define the solid contribution to the reaction force by considering the polymer network as an equivalent homogeneous solid.We define σ z as the normal stress in the equivalent solid specimen along the vertical direction z, z is the strain in the same direction.
In this work, we first indent the biphasic material at a very low indentation speed ( ż = 0.5 µm/s ≈ 0 for instance), which is considered slow enough to neglect any contributions of the fluid or solid viscoelastic effects.Therefore, the measured force in the spring is considered as the elastic contribution of the solid component.Thus, we define the elastic contribution of the solid to the reaction force by In addition to the force defined in ( 9), there may be another contribution coming from the viscoelastic component, for example, as well as other modeling inaccuracy that will be discussed in Section 5.

PGD solution of the multidimensional Darcy equation
The problem is solved in a cylindrical domain of height H and radius R with the following boundary conditions which consists of a zero velocity on the bottom support, the indentation velocity U at the top of the domain and zero pressure at r = R.The shown modeling results in a steady-state simulation problem.Since we aim to identify the permeability of the domain, one will have to solve the problem several times during the optimization algorithm for each value of the permeability.Moreover, at each time step, the same problem should be solved again, at each change in the value of the parameters U , R, or H . On top of this complication, computing the gradients numerically during an optimization problem will require several more solutions of the direct differential equation, when using a classical approach.Therefore, disposing of a parametric solution of the form P (x, y, z, H ,U , R, K) becomes a leverage, which we will use during the optimization problem.
In this simulation, we adopt a quasi-static approach while solving the Darcy equation, assuming the equilibrium is reached at each time step.Therefore the Darcy equation is written by where P is the pressure, µ the fluid's viscosity and K the domain permeability.
Another challenge is faced in this simulation, using a domain of initial height H 0 = 3 mm, while the indentation depth has the order of magnitude of few micrometers.Therefore, all the pressure variation exists only at the surface of the domain, within a thickness of few micrometers.Such problem is complicated using classical simulation techniques since it requires an extremely refined mesh in the thickness direction near the indented surface.The result is a large number of required degrees of freedom leading to unrealistic calculation time using classical calculation techniques.
However, using the PGD technique, we solve the problem using an in-plane-out-of-plane decomposition.Such decomposition results in solving the problem in (x, y) domain separately from the thickness z domain.Thus refining the mesh in the thickness direction as much as needed is possible without compromising the computation time [10,25].We may also use the PGD ability to solve problems using geometrical and configuration parameters as extra coordinates of the problem [26].This approach leads to a large problem dimensionality but is tractable using the PGD since it allows circumventing the curse of dimensionality by using the separation of variables [2].The solution is a parametric seven-dimensional solution of the pressure P defined in the following form.
We use in this simulation 100,000 nodes in the out-of-plane thickness direction z, logarithmically distributed, with an extremely refined mesh on the top surface of the domain, to avoid oscillation near the top surface, whereas we use 363 nodes in the in-plane domain.The permeability domain consists of 10,000 nodes logarithmically distributed to capture the behavior at low permeabilities.Moreover, we use 10,000 nodes in the thickness dimension H , 5000 in the radius dimension R and 1000 nodes in the imposed velocity dimension U .The PGD simulation resulted in a three-dimensional (3D) pressure field for any indentation velocity U , domain dimensions H and R and permeability K. Therefore, the results can be used on the fly, online, at any time step t ; once R, H , and U are computed from measured variables, to identify the domain permeability K, for example.The simulation is performed within 2 h on a normal portable PC, core i7 with 8 Gb RAM.The result obtained is the equivalent of 5 × 10 14 3D simulations, each with 36,300,000 degrees of freedom.The solution for one combination of these parameters is illustrated in Figure 2. Readers unfamiliar with the construction of parametric solutions and domain decomposition using the PGD are directed to [27,28] and the numerous references therein.What follows consists of a small introduction of the solution of the problem using the PGD.We first assume the function P as having a separated form known until an order n, and we attempt to define P n+1 using with R j , j = 1 . . .6 being the new functions to be identified defining the unknown terms in P n+1 , for instance the functions of (x, y), z, U , K, R and H . Since the geometrical dimensions are extra parameters of the PGD solution, the (x, y) domain is mapped into an apparent domain (ξ, ζ), a circle of center (0, 0) and radius 1.The thickness domain z is mapped into φ ∈ [0; 1].The mappings are defined by The Jacobian of the transformation becomes therefore Finally, the integral form of the problem becomes after replacing the (x, y, z) coordinates with the apparent domain ones Therefore, the exact notation of P is now Eventually the problem becomes nonlinear and its solution is obtained using a fixed point iterative algorithm.To compute R k , we first define the test function k noted P * k using Then replacing (18) into the weighted integral for of the problem leads to a two-dimensional problem in case k = 1, and a 1D problem for all other values of k.Once the fixed point iterative algorithm is converged for P n+1 , the newly found functions R j will be used to compute P n+2 and so on until the convergence of the residual of the governing differential equation of the problem defined in (6).
One may note that using the PGD, the gradient of P with respect to the process parameters are found in a straightforward manner, simply by deriving the function of the permeability L i (K).For instance

Experimental study
The indentation of the hydrogel poroelastic material was performed at five different probe velocities ż.The domain consisted of a relatively thick hydrogel, of thickness H 0 = 6 mm.The results are obtained using atomic force microscopy [24], and illustrated in Figure 3.As illustrated in Figure 3, the indentation depth has the order of magnitude of few micrometers.Therefore no sensible variation in the material properties are expected to appear in the studied domain.

Identification of the permeability and comparison of two models
In this section, we solve the minimization problem defined by where p solid is whatever property of the solid component to be identified from the minimization problem, F exp is the experimentally measured force in the spring, and F sim is the simulated force, which is defined as the sum of the contribution of the solid component and the fluid component of the reaction force.Two scenarios are considered for this work.

Modeling using a fluid in an elastic solid mesh
In this section, we solve (20) considering only the contribution of the elastic part of the solid component as defined in (9).Therefore, Equation ( 20) is a function of only one parameter, the permeability of the domain K since the solid properties are known for instance The minimization is performed using a classical Newton algorithm using the cost function to minimize C = (F exp − F sim ) 2 , and the unknown parameter to identify X = K are updated at each iteration using The Newton algorithm leads to the permeabilities shown in Figure 4.
As one may note from the results illustrated in Figure 4, that a small transition phase in the first 0.6 µm exists but later on the results are steady as a function of w.These results are expected since the indentation depth is small with respect to the thickness and therefore w should not affect the permeability values.On the other hand, one may notice a dependency of the permeability on the indentation velocity ż.Such result is unexpected and is a possible sign of the presence of a viscoelastic solid and/or other factors that may affect the material reaction.These factors are clearly linear with respect to the indenter velocity ż as illustrated in Figure 5.In fact, Figure 5 illustrates, as a function of ż, the average steady-state permeability value for each indentation velocity, from w = 0.6 µm until the end of the indentation experiment.To explain the speed  dependence of the apparent permeability shown in Figure 4, another term, proportional to the rate of deformation ˙ or ż, can be added to (21).

A nonlinear correcting term
In this section, we compare the experimental results to the simulation of a biphasic material with a nonelastic correcting term in the equilibrium equation.Using the results shown in Figure 5, one might assume that lumping the effects of material viscoelasticity into the permeability term generates inaccuracies that reflect as an apparent permeability that increases with speed.Assuming a constant permeability, we added another term that represents a rate to deformation dependent response of the solid polymer network into the equilibrium equation.Therefore, and since ż is proportional to ˙ , we suggest the following correction where γ is the correction factor representing the nonelastic response of the solid, σ z is the normal stress on the solid, and z is the normal strain in the z direction.The dependency of γ • ˙ z should be linear as a function of ż using this model, since w is comparable to z [17], and z is linear as a function of w.The cost function (24) to minimize is now written as The force in the solid component is now therefore defined as The cost function now has two unknowns to be identified.Therefore we use the Gauss-Newton algorithm to identify the two missing properties.Thus, we first define the following notation Figure 6.The identified viscosity of the solid part of the biphasic material as a function of ż.
A Top being the indented area on the top surface of the indenter.We may redefine the cost function to minimize by C such as where f i is the error function defined by (26) for the given indentation speed ż = i .Now the optimization problem is rewritten as Problem ( 28) is solved using Gauss-Newton algorithm.The expected result is a constant permeability K. Indeed, the identification yields a constant value of the permeability K ≈ 10 −16 but a negative viscosity γ as shown in Figure 6.This shows that the constant permeability is, in fact generating higher contributions to the reaction force than the real experimental one, and negative terms shall be added to go back to the correct value.

An informed neural network
In Section 5.2, the modeling error was corrected with a nonlinear term using a viscoelastic approach in the modeling of the solid mesh behavior.However, the results were not physical and the final model error using a constant permeability was still large.The modeling therefore, although did enhance the end outcome, remained suboptimal.With the presence of this human ignorance with respect to the behavior of the indented specimen, an informed artificial intelligence neural network fitting the modeling bias, becomes an appealing solution.In this section, we combine the fluid-reaction force model of Section 5 to a shallow neural network detailed in what follows, defining, therefore, a hybrid-twins approach, coupling the simulation results to a data-driven modeling technique, enhancing, therefore, the numerical simulation with experimental data.Eventually one can choose any fitting technique with enough parameters to represent the "ignorance", as neural networks are one of the fitting techniques.The author's choice is motivated only by the convenience and ease of use of the considered technique.
In this section, we use a shallow neural network made of 6 neurons in the first layer with Sigmoid activation function and one neuron in the output layer with a linear activation function.The model representation is given by Figure 7.In Figure 7, the first block to the left represents the model depicted in Section 5, which output is the model estimation of the fluid force F fluid using as hyper parameter the permeability K only.We define this modeling error in the section by where F fluid is computed using the PGD results as F fluid = A Top P • dA.Equation (29) represents the error the model defined in Section 5.1, where the simulated reaction force was defined as F sim = F spring@ ż≈0 +F fluid .Therefore E is the error of the elastic model.We will therefore attempt to use a neural network to model E .We define Y 1 as the output of the neural network, which should consist an estimation of E .We also define Y as an estimation of the experimental reaction force such as With this approach, eventually the fitting Y 1 of E with the neural network requires the knowledge of F fluid , the output of the PGD algorithm, and therefore a fitting of K. Therefore, a simultaneous fitting of both K and the hyper parameters Θ i of the defined neural network is the most suitable approach for this process.
To compute the problem's hyper parameters, the following cost function is defined as In (31), λ(1/2) n i =1 Θ 2 i consists of a regularization term to avoid over fitting, the parameters Θ i are the neural network hyper parameters, and λ is the regularization term.Later on, the minimization algorithm identifies the hyper parameters using To achieve convergence, we first converge the fluid force using a regularized Newton algorithm, as the gradient with respect to the fluid permeability K is much larger than the gradient  with respect to the neural network hyper parameter.For instance, we minimize the cost function J * defined by using an iterative regularized Newton update algorithm defined by The constant c 1 is a regularization constant to avoid numerical explosion of the algorithm, used equal to 1 in this work.The value of K converges within several iteration to 1.98×10 −16 Pa•s, a quite comparable value to the previous modeling effort.
Once the value of K is estimated, it will be used to the initiation of the minimization problem defined in (32).The problem is minimized using Levenberg-Marquardt gradient decent algorithm with a nonuniform learning rate.For instance, the learning rate of the hyper parameter K is lower than the learning rate of the neural network hyper parameters Θ i or biases b i , one bias per layer in the neural network.For training the neural network, the data is split into 70% for training set, 15% on the cross-validation set, and 15% for the test set.The hyper parameter λ is fitted to the cross-validation set.A batch gradient decent minimization is used to accelerate the convergence process.
The final estimated permeability K = 3.86 × 10 −16 , and the performance of the neural network fitting algorithm is illustrated in Figure 8.The mean square error is MSE = 0.48% on the best performance on the cross-validation set.As expected, the permeability is larger than the one identified in Section 5.2, which was a bit overestimating the fluid-reaction forces.

Discussion
Many factors may affect the modeling.Moreover, simulating all the local behavior around the indenter may become unacceptably expensive from a computational time point of view, taking into consideration the local effects and the thickness of the domain with respect to the indenter.Therefore, the adopted suitable approach is an empirical modeling of the error by including a  correcting factor in the equilibrium equation proportional to the strain-rate ˙ .In a first attempt, the parametric simulation of the Darcy equation is performed using the PGD and the Gauss-Newton optimization algorithm is used to identify the parameter γ, as identified in Section 5.2.
The obtained result γ is linear as a function of ż, which was expected from the reported error in Figure 5.The obtained permeability is constant, therefore, as a function of ż after the first transient region K ≈ 10 −16 , which is in fact less than the expected values shown in Figure 5.
However, the final results estimates a single value of the permeability.On a second approach, a physically informed neural network is used to correct the modeling error using a set of 25 hyper parameters (18 weights and 7 biases) to fit simultaneously the permeability using the model developed in Section 5 and the neural network parameters.The algorithm converges within few minutes on a normal portable PC, leading to an excellent performance of the model.The estimation modeling of the error E is illustrated in Figure 9, with a much better performance than the viscoelastic solid modeling explained in Section 5.2.
One may note that the availability of the PGD parametric solution in the separated form defined in (12) made the fully coupled minimization of J possible, and therefore possible the identification of the permeability K along with the hyper parameters of the neural network.
Otherwise, each iteration of the Levenberg-Marquardt algorithm would require at many direct solutions of the PDE problem defined in (6), which least to prohibitive and unrealistic calculation times.

Conclusion
In this work, we derived an expression of the error with respect to the previous work modeling biphasic materials as a combination of an elastic solid and a pressurized fluid.The suggested error modeling first identifies a strain-rate term and a parameter γ linear with respect to the indentation speed ż.Later on, a neural network is used for a superior fitting of the modeling error.The identification is performed by comparing a parametric simulation to an experimental atomic force microscopy nanoindentation of a biphasic soft hydrogel with a thickness several order of magnitude larger than the indentation depth.The simulation is performed by using an in-plane-out-of-plane decomposition allowing an impressive refinement of the mesh in the thickness direction and therefore identifying with high fidelity the pressure distribution inside the studied part.
The used hybrid-twins approach is only possible when all the solutions of the direct problem are readily available, which is the case in a model reduction framework, for instance, the PGD framework.The absence of the parametric multidimensional solutions leads to prohibitive calculation times during the optimization process.

Figure 1 .
Figure 1.The reference for z = 0 is taken at the point of contact of the indenter with the specimen, w(t ) is the penetration depth of the probe into the specimen.

Figure 2 .
Figure 2. The pressure (in Pa) of the 3D solution found using the PGD framework for a given combination of parameters.

Figure 3 .
Figure 3.The experimental AFM nanoindentation reaction force in the spring F (t ) as a function of the indentation depth w.

Figure 4 .
Figure 4.The identified as a function of the indentation depth w for different indentation velocities ż.

Figure 5 .
Figure 5.The identified permeability as a function of the indentation depth w for different indentation velocities ż.

Figure 7 .
Figure 7.An informed neural network to fit the ignorance of the model or the modeling error, Y being the experimentally measured force, F exp , while Y 1 is the output of the neural network consists of a fitting of the modeling error E .

Figure 8 .
Figure 8. Neural network fitting performance on the three sets: training, validation, and test sets.

Figure 9 .
Figure 9. E fitted by the neural network.