Elías Cueto
email: ecueto@unizar.es

David González

Francisco Chinesta

Learning Corrections for Hyperelastic Models From Data

Keywords: data-driven computational mechanics, hyperelasticity, model correction, GENERIC, machine learning

Unveiling physical laws from data is seen as the ultimate sign of human intelligence. While there is a growing interest in this sense around the machine learning community, some recent works have attempted to simply substitute physical laws by data. We believe that getting rid of centuries of scientific knowledge is simply nonsense. There are models whose validity and usefulness is out of any doubt, so try to substitute them by data seems to be a waste of knowledge. While it is true that fitting well-known physical laws to experimental data is sometimes a painful process, a good theory continues to be practical and provide useful insights to interpret the phenomena taking place. That is why we present here a method to construct, based on data, automatic corrections to existing models. Emphasis is put in the correct thermodynamic character of these corrections, so as to avoid violations of first principles such as the laws of thermodynamics. These corrections are sought under the umbrella of the GENERIC framework (Grmela and Oettinger, 1997), a generalization of Hamiltonian mechanics to non-equilibrium thermodynamics. This framework ensures the satisfaction of the first and second laws of thermodynamics, while providing a very appealing context for the proposed automated correction of existing laws. In this work we focus on solid mechanics, particularly large strain (visco-)hyperelasticity.

INTRODUCTION

In a very recent paper about how construct machines that could eventually learn and think like humans, [START_REF] Lake | Building machines that learn and think like people[END_REF] state that "machines should build casual models of the world that support explanations and understanding, rather than merely solving pattern recognition problems" and that "model building is the hallmark of human-level learning, or explaining observed data through the construction of causal models of the world". Indeed, machine learning of physical laws could be seen as the ultimate form of machine intelligence, and this should be done, of course, from data.

There is a very active field of research around this way of reasoning. For instance, in [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] a method is presented that operates on a bag of terms like sines, cosines, exponentials, etc., so as to find an expression that is sparse (i.e., it incorporates few of theses terms) while still explaining the experimental data. Similar approaches include techniques to find reducedorder operators from data (Peherstorfer andWillcox, 2015, 2016) or the possibility to construct physics-informed machine learning (Raissi et al., 2017a,b;[START_REF] Swischuk | Projection-based model reduction: formulations for physics-based machine learning[END_REF].

In the field of computational materials science, this approach seems to begin by the works of Kirchdoerfer andOrtiz (2016, 2017a). In it, and the subsequent works, they present a method in which the constitutive equation is substituted by experimental data, that could be possibly noisy [START_REF] Kirchdoerfer | Data driven computing with noisy material data sets[END_REF][START_REF] Ayensa-Jiménez | A new reliability-based data-driven approach for noisy experimental data with physical constraints[END_REF]. In them, it is recognized that some equations (notably, equilibrium, compatibility) are of a higher epistemic nature, while constitutive equations-that are often phenomenological and, therefore, of lower epistemic valuecould easily be replaced by data [START_REF] Latorre | What-you-prescribe-is-whatyou-get orthotropic hyperelasticity[END_REF]. The criterion is to establish a distance measure that indicates the closest experimental datum to be employed every time the constitutive law is called at the finite element integration point level.

In some of our previous works, this approach is further generalized by defining the concept of constitutive manifold, a low-dimensional embedding for the stress-strain pairs (see [START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF]. Thus, by alternating between stress-strain pairs that satisfy either equilibrium or the constitutive equation, the solution that satisfies the three families of equations is found, regardless of the non-linearity of the behavior. Several methods have been studied for the construction of this constitutive manifold [START_REF] Ibañez | Data-driven non-linear elasticity: constitutive manifold construction and problem discretization[END_REF].

Another inherent difficulty in trying to machine learning models is that of the adequate level of description. Every physical phenomenon can be described at different levels of detail. In the case of fluid mechanics, for instance, these levels range from molecular dynamics to thermodynamics-in descending order of detail-. In between, different theories have been developed that take care of different descriptors of the phenomenon taking place: from the Liouville description to the Fokker-Planck equation, hydrodynamics, ... to name but a few of the different possibilities [START_REF] Español | Statistical Mechanics of Coarse-Graining[END_REF]. Thus, there should be a compromise between detail in the description and the resulting computational tractability of the approach. This is something very difficult to discern for an artificial intelligence.

The risk of employing an approach based upon pure data regression is to violate-due to the inherent noise in data, for instance-some basic principles such as the laws of thermodynamics: conservation of energy, positive dissipation of entropy. Trying to avoid these possible inconsistencies, in [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] we developed a data-driven method that operates under the framework of the GENERIC formalism [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. I. development of a general formalism[END_REF][START_REF] Öttinger | Beyond Equilibrium Thermodynamics[END_REF]. The General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) constitutes a generalization of the Hamiltonian mechanics. Therefore, under the GENERIC umbrella, the equations satisfy basic thermodynamic principles by construction.

Thus, the problem translates to finding-by means of datathe right expression of the particular GENERIC formalism for the system at hand (or its finite element approximation, if we work in a purely numerical framework). The resulting approximation is thermodynamically sound and very appealing from the numerical point of view. The stability of the GENERIC approach and its thermodynamic consistency-in particular, the conservation of symmetries in the formulation-has been thoroughly investigated in previous works, whose lecture is greatly recommended [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF][START_REF] Romero | Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part i: Monolithic integrators and their application to finite strain thermoelasticity[END_REF].

However, even if the usual parameter fitting procedure from experimental data is often painful and, notably, gives poor fitting of the results in many occasions, we believe that well-known constitutive equations should not be discarded, thus waiting centuries of scientific discovery. Instead, we believe that it is interesting to simply correct those models that sometimes do not fit perfectly the results-sometimes locally, in a delimited region of the phase space-. This is the approach followed in [START_REF] Ibañez | Hybrid constitutive modeling: data-driven learning of corrections to plasticity models[END_REF], where corrections are developed to yield criteria so as to render them compliant (to a specified tolerance level) with the available experimental results. A similar approach has been pursued recently in [START_REF] Lam | Should you derive, or let the data drive? an optimization framework for hybrid first-principles datadriven modeling[END_REF] for a study on the interaction of aircraft wings. In these approaches, the chosen level of description is defined by the (poor) model, so, in principle, no further decision needs to be taken, as will be discussed later on.

The GENERIC formalism is valid for all levels of description, and could also help in deriving corrections from data that still maintain the thermodynamic properties of the resulting model. Hyperelastic models fall within Hamiltonian mechanicsi.e., they represent a purely conservative material-. However, rubbers or foams usually present some degree of viscoelasticity, for instance. In this framework, Hamiltonian mechanics will no longer be the right formalism to develop their constitutive equations. GENERIC should be preferred instead.

In this paper we study how to learn these corrections from data. First, in Section 2 we review the basics of the GENERIC formalism, with an emphasis on hyperelastic and visco-hyperelastic materials. In Section 3 we explain how to employ GENERIC to develop corrections to existing models from data, while in Section 4 we introduce, by means of an academic example in finite dimensions, the basic ingredients of our approach. This will be further detailed for visco-hyperelastic materials in Section 5. The paper ends with a discussion on the just developed techniques and the future lineas of research in Section 6.

A REVIEW OF THE GENERIC FORMALISM

The Basics

The GENERIC formalism was introduced by [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. I. development of a general formalism[END_REF] in a seminal paper in an attempt to give a common structure for non-Newtonian fluid models. The establishment of such a model in the GENERIC framework starts by selecting appropriate state variables. This is not straightforward in a general case in which we have no prior information about the precise behavior of the system at hand. However, for most systems-and specially when we start from known models, as it is the case in this work-simple rules exist for the selection of such variables (Öttinger, 2005). Selecting mutually dependent variables does not constitute a problem, in fact, as most of the literature on GENERIC demonstrates. Let us call these variables z t = z(t) : I → S, z ∈ C 1 (0, T], and emphasize their obvious time dependency in the interval I = (0, T]. S represents the space in which these variables live, which depends obviously on the particular system under scrutiny. The final objective of the GENERIC model is to establish an expression for the time evolution of these variables, ż(t).

The GENERIC equation takes, under these assumptions, the form

żt = L(z t )∇E(z t ) Hamiltonian + M(z t )∇S(z t ) Dissipative , z(0) = z 0 .
(1)

The first sum on the right-hand side term represent the Hamiltonian, or conservative, part of the behavior of the system. In it, the term L(z t ) is the so-called Poisson matrix. The second sum is responsible for the dissipative behavior of the system, with M(z t ) the so-called friction matrix. Here, E(z t ) represents the total energy of the system, while S(z t ) represents its entropy.

For Equation ( 1) to give a valid description of any physical system, it must be supplemented with the so-called degeneracy conditions:

L(z) • ∇S(z) = 0, (2a)
M(z) • ∇E(z) = 0. ( 2b
)
Enforcing these conditions leads to the necessity of L(z) to be skew-symmetric and a M to be symmetric, positive semi-definite. If these conditions are met, then, it holds,

Ė(z) = ∇E(z) • ż = ∇E(z) • L(z)∇E(z) + ∇E(z) • M(z)∇S(z) = 0, (3)
which is, in fact, the equation of conservation of energy for the system. Additionally, these conditions ensure the satisfaction of

Ṡ(z) = ∇S(z) • ż = ∇S(z) • L(z)∇E(z) + ∇S(z) • M(z)∇S(z) ≥ 0,
(4) or, equivalently, the fulfillment of the second principle of thermodynamics.

Noteworthy, Equation (1) constitutes the most general framework to develop a valid constitutive equation in the light of the principles of thermodynamics. A valid constitutive model must satisfy the GENERIC equation, and any possible correction to it should not deviate the result from this framework. For a thorough description of a long list of models under the GENERIC formalism, the interested reader can consult Öttinger (2012). To exemplify the just introduced concepts, consider the simplest case of a conservative mechanical system whose time evolution can be expressed, in the Hamiltonian framework, by resorting to a description of the type żt = {q t , p t }, where q t represents the position and p t the momentum. In that situation, the system is purely Hamiltonian and

L(z) = 0 1 -1 0 ,
with no entropy evolution, i.e., M = 0. In this simple situation, L(z) turns out to be the canonical symplectic matrix and the GENERIC description of the system reduces to that of a Hamiltonian system.

Hyperelasticity Under the Prism of GENERIC

It is important to highlight the fact that, for hyperelastic materials, the expression

żt = L(z t )∇E(z t )
represents the usual hyperelastic problem under the Hamiltonian formalism [START_REF] Romero | A characterization of conserved quantities in non-equilibrium thermodynamics[END_REF]. Indeed, if we choose z(x, t) = [x(X, t), p(X, t)] ⊤ , where x = φ(X)-the deformed configuration of the solid-and p represents the material momentum density, then,

ż = ẋ ṗ = L∇E = L ∂E ∂x ∂E ∂p
.

The total energy of an elastic body can be decomposed as

E = W + K,
i.e., the sum of elastic and kinetic energies. Here, we assume a strain energy density potential w of the form

W = w(C) d ,
where C represents the right Cauchy-Green deformation tensor. While, in general, the strain energy density for an isotropic case would be of the form w = w(X, C, S), in the context of isotropic hyperelasticity-a purely Hamiltonian case-, this dependence is often simply w = w(C). In turn, the kinetic energy will be

K = 1 2ρ 0 |p| 2 d .
In this framework, it is clear that

∂E ∂x = ∂W ∂x = ∇ X • P = ∇ X • [FS],
where P and S represent, respectively, the first and second Piola-Kirchhoff stress tensors and F is the deformation gradient.

Given that

p = ρ 0 V = ρ 0 ∂x ∂t ,
with V the material velocity and ρ 0 the density in the reference configuration so that, finally,

ż = ẋ ṗ = L∇E = L ∇ X • P p ρ 0 .
This implies that

L = 0 3×3 I 3×3 -I 3×3 0 3×3 , (5)
which is fully compliant with the GENERIC framework, see Equation ( 1). This model is readily seen as equivalent to

ẋ = p ρ 0 , ∇ X • P = ṗ,
which correspond to the definition of the material momentum density and the equilibrium equation, respectively. Under this rationale, the possible viscous effects in the material would be described by the second sum in Equation (1).

REMARK. We have stated that, under the GENERIC formalism, an isotropic Hamiltonian or conservative hyperelastic model can be written in the form w = w(C) and therefore will not depend on S. This discussion is strongly related with that of the adequate level of description of the model. In fact, many hyperelastic models exist that depend on different parameters, that can influence its viscous behavior, for instance, see [START_REF] Mihai | How to characterize a nonlinear elastic material? a review on nonlinear constitutive parameters in isotropic finite elasticity[END_REF].

Indeed, by introducing a new potential (entropy) in the formulation, what we are doing is to introduce ignorance on these details, while still taking into account their influence on the results. It is the same process we face if we are not interested in tracking every molecule of a gas in a container but prefer instead a description based on macro-scale magnitudes such as pressure, volume, and temperature. The process of coarsegraining the description in a non-equilibrium setting makes it necessary to introduce a new potential that accounts for the neglected information: entropy [START_REF] Español | Statistical Mechanics of Coarse-Graining[END_REF][START_REF] Pavelka | Multiscale Thermodynamics[END_REF]. Thus, in the correction procedure that we are about to introduce, there will be no need to add new variables to the model, but an adequate entropy potential to the formulation.

The problem of constructing a valid constitutive model under the GENERIC point of view is therefore reduced to that of finding the particular structure of the terms L(z), E(z), M(z), and S(z). The classical approach is to do it analytically, as in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF][START_REF] Romero | Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part i: Monolithic integrators and their application to finite strain thermoelasticity[END_REF], for instance, or Vázquez-Quesada et al. (2009) and [START_REF] Español | Statistical Mechanics of Coarse-Graining[END_REF], to name but a few of the examples in the literature. A different approach is to find the structure of these terms numerically, from data. This will be done possibly with the help of manifold learning techniques such as LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] or isomap [START_REF] Tenenbaum | A global framework for nonlinear dimensionality reduction[END_REF], among others. It is the approach followed by the authors in [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] and, in some sense, it is also the approach followed by [START_REF] Millán | Nonlinear manifold learning for model reduction in finite elastodynamics[END_REF] without even knowing the structure of GENERIC. This approach is also somehow related to the use of compositional rules to construct models [START_REF] Grosse | Exploiting compositionality to explore a large space of model structures[END_REF]. This last reference shares with the approach herein the need of identifying the structure of several matrices that are then used to develop models-in that case, of phenomena that do not even obey the laws of physics, such as voting tendencies, for instance.

CORRECTING MODELS IN A GENERIC FRAMEWORK

In this work we do not pursue to unveil models by means of GENERIC and experimental data. As explained in the introduction, we believe that is simply nonsense to discard models that have demonstrated to be useful for decades. In the case of hyperelasticity, these include, among a wide list of references, the works of [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], [START_REF] Ogden | Non-linear Elastic Deformations[END_REF] or [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]. These models, as analyzed before, already had a GENERIC structure.

Purely hyperelastic materials are strictly conservative. However, soft living matter, for instance, that is often modeled under the hyperelastic theory, present some non-negligible viscous effects [START_REF] Peña | A formulation to model the nonlinear viscoelastic properties of the vascular tissue[END_REF][START_REF] García | Viscoelastic properties of the passive mechanical behavior of the porcine carotid artery: Influence of proximal and distal positions[END_REF]. In that case, in the light of the GENERIC formalism, it is necessary to complement the model with a dissipative part, i.e., to determine the precise form of M(z) and S(z).

What we will do in this work, in fact, is to assume that an inexact model exists, so that a correction is needed,

z corr = z exp -z mod ,
where "corr", "exp" and "mod" stand, respectively, for correction, experimental and model. We will develop a correction in the GENERIC framework so as to guarantee that the corrected model for the experimental results will also have a GENERIC structure. To this end, we cast the correction in the form żcorr = L∇E(z corr ) + M∇S(z corr ).

We do not consider a correction for L nor M, since, in the light of the previous remark, L is assumed to be identical to that of the model (we consider the same state variables). Since the correction of the model could have an important influence on the form of M-recall again the remark in the previous section, we attribute to S the possible presence of fine-grained state variables that are not considered in the Hamiltonian part of the model-, we discard any possible M coming from the inexact model and instead re-compute it from scratch. With these assumptions, the resulting model that fits with the experimental results will have the form żexp = żmod + żcorr = L∇E(z corr ) + M∇S(z corr ) + L∇E(z mod ), so that, finally, żexp = L ∇E(z corr ) + ∇E(z mod ) + M∇S(z corr ), which proves that the corrected model for z exp possesses a GENERIC structure with a correction in the Hamiltonian term.

Consider that a set of n meas experimental measurements Z = {z exp 0 , z exp 1 , . . . , z exp n meas } is available. The predictions of the inexact model are then subtracted from the experimental results. The final objective will be therefore to obtain a discrete approximation

z corr n+1 -z corr n t = L DE(z corr n+1 ) + M(z corr n+1 )DS(z corr n+1 ),
to the GENERIC structure of the discrepancy between data and experiments, by identifying DE(z), and possibly also M(z) and DS(z). DE and DS represent the discrete gradients (in a finite element sense). Therefore, the proposed algorithm will consist in solving the following (possibly constrained by the degeneracy conditions) minimization problem within a time interval J ⊆ I:

µ * = {M, DE, DS} = arg min µ ||z(µ) -z meas ||,
with z meas ⊆ Z, a subset of the total available experimental results. See the discussion in [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] about how to determine the right size of the sample set, the possibility of employing monolithic or staggered strategies, etc.

In the next Section this procedure is exemplified with the help of an academic example in finite dimensions.

AN INTRODUCTORY EXAMPLE

We first consider an example analyzed in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF] and then again in [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF]. The system is a double pendulum, which is connected by thermoelastic springs. It comprises two masses m 1 and m 2 connected by springs of internal energy e 1 and e 2 . They oscillate around a fixed point, see (Figure 1). We employ the classical notation of Hamiltonian mechanics where q i , p i , i = 1, 2 represent position and momenta, respectively. For the springs, their respective entropies are s j , and the longitudes at rest will be denoted by λ 0 j , j = a, b. The set of state variables for this double pendulum will be therefore

S = {z = (q 1 , q 2 , p 1 , p 2 , s 1 , s 2 ) ∈ (R 2 × R 2 × R 2 × R 2 × R × R), q 1 = 0, q 2 = q 1 }.
The GENERIC structure for this problem needs to consider the internal energy of the system. Again, the internal energy is composed by the kinetic energy of the masses and the potential energy in the springs, i.e., E(z) = K 1 (z) + K 2 (z) + e a (λ a , s a ) + e b (λ b , s b ), with λ a = √ q 1 • q 1 , λ b = (q 2q 1 ) • (q 2q 1 ). The temperature in the springs, θ j , is assumed to be originated by the Joule effect,

θ j = ∂e j ∂s j , j = a, b.
The conductivity in the springs will be denoted by κ. Under this rationale, the resulting equations for the double pendulum will be

qi = p i m i , ṗi = - ∂ ∂q i (e a + e b ), ṡj = κ θ k θ j -1 , with i = 1, 2, j = a, b, k = j.
Therefore, the gradients of the GENERIC formalism will look

∇E(z) = f a n a -f b n b , f b n b , p 1 m 1 , p 2 m 2 , θ a , θ b , (6a)
∇S(z) = (0, 0, 0, 0, 1, 1), (6b) with f j , n j , j = a, b, the forces in the springs and their respective unit vector along their direction.

Poisson and friction matrices will result in this case,

L(z) =        
0 0 1 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

        , M(z) =         
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 κ θ b θ a -κ 7) However, we will assume that this description of the system is not available-will be used as a ground truth to determine errors-and that the system is thought to be purely Hamiltonian.

0 0 0 0 -κ κ θ a θ b          . (
In this scenario, the goal of our method will be that of unveiling the dissipative part of the model so as to correct the pure Hamiltonian behavior of the assumed model. In other words, the system will be considered as modeled by żt = L(z t )∇E(z t ), with L as in Equation ( 7) and ∇E as defined in Equation (6a). Results of the ground truth, the assumed (purely Hamiltonian) model and the found corrected model are shown in Figure 2.

The mean squared error of the assumed model with respect to the pseudo-experimental data was initially 0.1732%. Note the little influence of the Joule effect on the results. However, after a correction is found and the dissipative character of the model is taken into account, this error is decreased up to 0.0125%, i.e., one order of magnitude.

CORRECTIONS TO HYPERELASTIC MODELS

In order to show the full capabilities of the proposed method, we consider now an example of a visco-hyperelastic material whose precise constitutive model is to be corrected from experimental data.

Ground Truth. Pseudo-Experimental Data

The pseudo-experimental data is obtained by finite element simulation of a visco-hyperelastic Mooney-Rivlin material in which

W = C 1 (I 1 -3) + C 2 (I 2 -3) + D 1 (J -1) 2 , ( 8
)
with I 1 = J -2 3 I 1 and I 2 = J -4 3 I 2 , and where the invariants of the right Cauchy-Green tensor C are defined as I 1 = λ 2 1 +λ 2 2 +λ 2 3 , and

I 2 = λ 2 1 λ 2 2 + λ 2 2 λ 2 3 + λ 2 3 λ 2 1
, respectively. J represents, as usual, the determinant of the gradient of deformation tensor. In this case, C 1 = 27.56 MPa, C 2 = 6.89 MPa and D 1 = 0.0029 MPa.

To model the viscoelastic behavior of this rubberlike material, it is assumed that the material's shear modulus G and bulk modulus K evolve in time. This evolution is modeled by means of a Prony series in terms of the instantaneous moduli,

G(t) G 0 = 1 - 2 i=1 g P i 1 -exp - t τ i , K(t) K 0 = 1 - 2 i=1 k P i 1 -exp - t τ i , with g P i = [0.2, 0.1] and k P i = [0.5, 0.2].
The relaxation times take the values τ i = [0.1, 0.2] seconds, respectively. With these values, the initial instantaneous Young's modulus takes the value E = 206.7 MPa, with Poisson's ratio ν = 0.45.

Data was generated after a total of 557 different loading processes to the same specimen. It was subjected to a load history of different amplitudes. In every case, a first plane stress state (σ x , σ y , τ xy )-values are not correlated-is applied during a short impulse of 0.021 seconds, then maintained at constant value for one more second, allowing the material to creep. This is followed by a second loading process of 0.021 seconds at a different (σ x , σ y , τ xy ) value, followed by a final plateau of one more second. For each one of the 557 different experiments these two stress states were different. These results are stored in the form of 557 different Z vectors, thus representing a trajectory in time.

Modeling the Results With a Purely Hyperelastic Model

After the generation of the pseudo-experimental data, we tried to reproduce these results with a deliberately wrong model: the material was assumed to be modeled by a Mooney-Rivlin model with no viscous response (and thus purely Hamiltonian or conservative). The comparison of the experimental results and the predictions given by this (poor) model are shown in Figure 3.

It seems obvious that a classical Mooney-Rivlin model can not reproduce the viscous behavior of the reference material. In the next section a correction to this model is developed based on the available data and the procedure introduced in Section 3.

Correction of the Dissipative Part of the Model

Knowing in advance that the pseudo-experimental results come from a viscous modification to a Mooney-Rivlin model, a first attempt is made of finding a correction by incorporating a dissipative part in the GENERIC description of the model. To this end, for each one of the experimental results, a fitting procedure of the dissipative GENERIC terms was accomplished.

In Figure 4 results are shown for one of the 557 essays. Experimental results, Mooney-Rivlin prediction and the subsequent GENERIC correction are shown. As can be noticed, experimental results are captured to a high degree of accuracy. In this case, for the particular test shown in Figure 5, the mean squared error was 0.018%. All the tests showed similar levels of error.

What if Some Terms Need no Correction?

Of course, in general we will not know in advance that a particular model is the best for the Hamiltonian part of the behavior. In a general situation both parts of the model will need to be corrected. To show the robustness of the presented method, we demonstrate here that if we try to correct the Hamiltonian part of the model, the method is able to detect that it is already correct (Mooney-Rivlin) and that it needs no correction. The method proceeds by correcting the dissipative part only, obtaining the same levels of error as the preceding section.

Constructing the Good Model

The final goal of the method is not to reproduce each one of the experimental results, but to be able to construct a true model from data. To this end, we first unveil the underlying manifold structure of the experimental data. The temporal series of z exp (t) is grouped into a high dimensional vector, one for each of the 557 experiments. These are then embedded, by means of Locally Linear Embedding techniques [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] onto a low-dimensional manifold. This permits to unveil the true neighborhood structure between experimental data and, notably, to perform rigorous interpolation on the manifold structureand not on the Euclidean space-among data.

The first step when applying LLE techniques to a set of high-dimensional data is to find the right dimensionality of the embedding space. To do so, the eigenvalues of the projection matrix are usually studied. These are depicted in Figure 5.

The first LLE eigenvalue is always close to zero within machine precision, and is discarded. The next "isolated" eigenvalues represent the true dimensionality of the embedding space (in this case, three). The rest of the eigenvalues are usually much closer to each other and do not represent the right dimensionality of the embedding space. Therefore, it seems that the right dimensionality of the embedding space is three-even two.

Locally Linear Embedding techniques need some user intervention to determine, by trial and error, the adequate number of neighbors for each datum. In this case we assume some 20 neighbors for each one. The key step in finding the good low-dimensional embedding of the data is to find a vector of weights W that minimizes the functional

F (W) = 557 m=1 z - 20 i=1 W mi z i 2 .
Once these weights are found, LLE assumes that they continue to be valid in the low-dimensional embedding, and looks for the new coordinates ξ in this space accordingly, by minimizing a new functional

G(ξ 1 , . . . , ξ 557 ) = 557 m=1 ξ m - 20 i=1 W mi ξ i 2 .
This procedure allows us to find the constitutive manifold, as defined in [START_REF] Ibañez | Data-driven non-linear elasticity: constitutive manifold construction and problem discretization[END_REF]. It is shown in Figure 6. The objective of this validation procedure will be to try to reproduce a control point in the manifold-a complete loading history, in fact-by obtaining its GENERIC model from the neighboring experimental points. This control point is shown in red in Figure 6.

In Figure 7 the result of the interpolated model (in continuous line) and the eight neighboring experimental results (dashed lines) that served to construct the final GENERIC model for the red point in Figure 6 are shown. The mean squared error with respect to the control experimental history resulted be 0.174%.

Full Model Correction

In the preceding sections we assumed that the Hamiltonian part of the model (basically, a Mooney-Rivlin model) was known and that the model needed only some amendment in its dissipative part. In this section we study the performance of the proposed technique if every term in the assumed model is wrong.

To this end, we assume for the solid a Neo-Hookean model with no viscous dissipation. The neo-Hookean model is basically equal to Mooney-Rivlin, see Equation ( 8), with C 2 = 0. To make things even more difficult, we assume a bad calibration of the instruments so that, for this "wrong" model, C 1 = 68.9 MPa (four times the right value for the Mooney-Rivlin model-the actual one-) and D 1 = 0.0016.

Proceeding like in previous sections, we first computed corrections for each one of the 557 different experimental time series. For one of these essays, the prediction given by For this particular case (every experiment provided similar results), the initial error for the prediction given by the "wrong" neo-Hookean model was 13.05%. After correction, the relative mean square error in the time history was 0.092%.

Once the whole 557 experiments have been corrected, the constitutive manifold for this material can be constructed by LLE methods, as detailed in Section 5.5.

With this constitutive manifold thus constructed we can now evaluate the behavior of any new strain-stress state by simply locating it in the manifold, determining its surrounding neighbors, and employing the LLE weights to interpolate its GENERIC terms. This was done for one of the experimental results, that was removed from the manifold for control purposes, and interpolated from its neighbors. The result of this process is shown in Figure 9.

The mean squared error along the time history with respect to the control experiment was 1.057%.

DISCUSSION

From the results just presented, it is clear that the proposed technique presents an appealing alternative for the machine learning of models from data. Instead of constructing data-driven models from scratch, constructing only corrections to existing, well-known models has shown to provide very accurate results that very much improve these models. One key ingredient in these developments is the concept of constitutive manifold, that allows to interpolate experimental results in the right manifold structure. Existing works choose simply the nearest experimental neighbor, but, notably, this neighborhood is found in an Euclidean space [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] or in a Mahalanobis space [START_REF] Ayensa-Jiménez | A new reliability-based data-driven approach for noisy experimental data with physical constraints[END_REF].

The presented method is robust even if some parts of the model need no correction. The final method, as has been presented, has the important property of being sound from the thermodynamic point of view, guaranteeing, thanks to its GENERIC structure, the conservation of energy and positive production of entropy.

From the numerical point of view, the resulting, GENERICbased time integrator schemes have already demonstrated their ability to conserve the right symmetries of the system (see, for instance, [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF][START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF]. In sum, we believe that the just presented technique, that should be extended to other types of systems, presents a promising future.

FIGURE 1 |

 1 FIGURE 1 | Double thermal pendulum.

FIGURE 2 |

 2 FIGURE 2 | Results for the thermal pendulum problem. Results are shown (see the detail in the small window in the bottom figure) for the ground truth (pseudo-experimental data), the uncorrected (purely Hamiltonian) assumed model and the corrected one.

FIGURE 3 |

 3 FIGURE 3 | Loading process for one particular experiment. Pseudo-experimental response and prediction made by the standard (non-viscous) Mooney-Rivlin material.

FIGURE 4 |

 4 FIGURE 4 | Comparison of the Mooney-Rivlin model prediction and its subsequent GENERIC correction with the experimental results for one particular experiment.

FIGURE 5 |

 5 FIGURE 5 | Evolution of the eigenvalues of the projection matrix in the embedding of experimental data. Only the first 17 eigenvalues are shown for clarity.

FIGURE 6 |

 6 FIGURE 6 | Obtained constitutive manifold by embedding the experimental results onto a three-dimensional space. Only a portion of the 557 experimental results are shown for clarity. In red, control point employed to validate the approach. Note that it is surrounded by a user-defined number of neighbors, whose GENERIC model is employed to obtain, by interpolation by means of the LLE weights, the sought model.

FIGURE 7 |

 7 FIGURE 7 | Result (continuous line) of the interpolation on the constitutive manifold of the eight neighboring experimental results (in dashed line). These are the time history of the eight neighboring points in Figure 6.

  the "wrong" (neo-Hookean) model, the experimental results (coming from the Mooney-Rivlin model) and the corresponding corrected model predictions are shown in Figure 8.

FIGURE 8 |

 8 FIGURE 8 | Experimental results (circles), prediction made by the neo-Hookean model (dashed line) and corrected model (red line) for experiment number 85.

FIGURE 9 |

 9 FIGURE 9 | Experimental results (circles) for a new experiment, prediction made by the neo-Hookean model (dashed line) and interpolated corrected model (red line). The interpolation is made by employing the same weights provided by LEE techniques in constructing the constitutive manifold.

Frontiers in Materials | www.frontiersin.org

February 2019 | Volume 6 | Article 14

FUNDING

This work has been supported by the Spanish Ministry of Economy and Competitiveness through Grants number DPI2017-85139-C2-1-R and DPI2015-72365-EXP and by the Regional Government of Aragon and the European Social Fund, research group T24 17R.

DATA AVAILABILITY STATEMENT

The datasets [GENERATED/ANALYZED] for this study will be available upon request.

AUTHOR CONTRIBUTIONS

DG, FC, and EC contributed in the conception and design of the study. DG coded the Matlab program. EC wrote the first draft of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.

Conflict of Interest Statement:

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.


Lagarde Antoine

Josserand Christophe

Protière Suzie
email: protiere@ida.upmc.fr

The capillary interaction between pairs of granular rafts †

When an object is placed at the surface of a liquid, its weight deforms the interface. For two identical spherical objects, such a deformation creates an attractive force, leading to the aggregation of the two-body system. Here, we experimentally investigate the interaction between two granular rafts, formed by the aggregation of dense millimeter-sized beads placed at an oil-water interface. The interfacial deformation created by such a two-dimensional object exceeds by at least an order of magnitude the deformation of a single bead. This leads to unusually high capillary forces which strongly depend on the number of particles. Likewise, because the raft grows in size as more particles are added, the viscous drag experienced increases along with the capillary attraction, leading to a non-trivial dependence of the balance of forces on the number of beads. By studying the relative motion of two granular rafts in relation with the interfacial deformation they generate, we derive a model for the observed speed profiles. With this work, we generalize how the capillary interaction between two non-identical complex structures evolves with their respective geometry.

Introduction

When a few identical objects are deposited on a fluid interface, a daily yet striking phenomenon occurs, during which the passive objects will spontaneously move towards one another. The implications of what could appear as a trivial experiment are countless. In Nature, fire ants use it to survive flood by aggregating into a floating raft [START_REF] Mlot | Proceedings of the National Academy of Sciences[END_REF] , while some animals are able to propel by deforming the water surface in a specific way [START_REF] Hu | [END_REF] . This phenomenon is also used in industrial processes for the fabrication of macroscopic objects with a specific microstructure [3][4][5] .

The self-assembly of objects at a liquid interface relies on the interaction due to the surface deformation around each particle 6 . When the menisci around several objects overlap, capillary attraction (or repulsion) arises. Colloidal particles generate lateral forces due to the undulation of the contact line around their surface [7][8][9][10] , either because of surface roughness and wettability, or due to a non-spherical shape. The interface around the particle is then conveniently seen as a capillary multipole. For anisotropic colloids, the interaction depends strongly on the individual shapes 11 , giving rise to oriented aggregation either through a precise manufacturing of the form of each particle 12 , or thanks to a curved interface acting as a force field 13 .

For larger objects, another contribution comes from the parti-cles' weight which deforms the surrounding interface by pushing the liquid down, an effect that can be taken into account through a generalized Archimedes' principle 14 . Models for such interactions have been developed almost exclusively for two identical rigid bodies 15,16 . In 1949, Nicolson 15 proposed a model using a linear superposition approximation of the interface deformations to describe the force of interaction between two floating bubbles at a liquid interface. This far-field calculation leads to an analytical formulation of the capillary force between identical objects (bubbles, solid sphere, ...) 16,17 . Good experimental agreement has been reported for the relative motion of a small number of large identical spheres 18,19 . However, most of these models and experiments concern individual and identical particles and no detailed investigations have been made so far to characterize the interactions between clusters of particles.

In this paper, we focus on such large objects, where dense millimeter-sized beads act as monopoles. We study the interaction between two granular rafts 20,21 , each formed by the aggregation of dense spherical beads into a single axisymmetric object. Such 2-D structures are flexible and exhibit high deformations of the interface due both to their weight and their geometrical extension. The overall meniscus created by these large objects generates a specific raft-raft interaction that needs to be characterized. We quantify the dependence of the drag and the capillary force experienced by a raft on both its size and the size of the other interacting raft. The present study questions whether two assemblies of particles forming dense objects can interact as single entities, and investigates the influence of the clusters di-mensions on their interactions.

Kinematics of two interacting granular rafts

Experimental setup

In a typical experiment, a thick layer of oil (kinematic viscosity ν = 50.10 -6 m 2 .s -1 , density ρ o = 960 kg.m -3 , oil-water surface tension γ = 38 mN.m -1 ) is carefully poured into a tank of dimensions 0.2 x 0.2 x 0.25 m filled with pure water. A precise number of particles n A are sprinkled from above. Beads are deposited far away from one another, from every direction around the raft. Due to the very strong gravitational interaction between such objects the beads automatically aggregate into a compact axisymmetric raft. When using the beads presented in table 1, we never observe a loosely packed assembly such as the one described in 22,23 . The raft size can be controlled by adding beads progressively (meaning the raft is made of precisely n A particles), as represented in figure 1.

The same procedure is reproduced elsewhere in the tank, with n B particles. The motion of the two rafts thus formed (respectively identified by the letters A and B) is recorded either from above or from the side by a camera at 250 frames per second. The properties of the various beads used in the experiments are summarized in table 1, along with the approximate maximum number of particles n max the corresponding raft can reach before sinking 21 . At the oil-water interface, the capillary length c = γ/((ρ wρ o )g) is greatly increased by buoyancy effects, so that c ≈ 10 mm (with ρ w the density of water, and g the gravitational acceleration). The maximal possible deformation is thus far more important than at an air-water interface. This leads to unusually high long-range capillary forces.

Speed measurement

For a given set of sizes (n A , n B ), we record the dynamics from above to measure the speed of the rafts, and from the side when we want to access the deformation of the interface. A typical experiment can be visualized in figure 2 (see supplementary information for the movies, taken from above and from the side), for two rafts with respectively n A = 30 and n B = 50. The two rafts move towards one another until they come into contact, at which point they will rearrange to form a bigger raft or sink. Here, we focus on the interaction of two rafts before they reach contact.

In a typical experiment such as the one presented in figure 2.a, the fluid interface undergoes strong deformations due to the weight of each raft. This induces a vertical displacement which cannot be neglected in our measurement of the raft speed. For example, in figure 2.a, the amplitude of the vertical motion of raft B (n B = 30) between the last two images is more than half the length of its horizontal motion. In our study, all the speed measurements are derived from top views of the granular raft motion. Yet, a movie taken from above only gives us information on the horizontal projection of the speed. The vertical component is not directly accessible. We deduce the vertical displacement from the radial motion, via measurements of the interface deformation.

The equation for the interface is obtained by a classic equilibrium between the hydrostatic pressure and the pressure jump due to the curvature. For a cylindrical coordinate system centered in the middle of the raft, the height of the interface beyond the raft B is obtained as the solution of the following system of equations:

     γ h 1+h 2 + h r -(ρ w -ρ o )gh √ 1 + h 2 = 0 h(r → ∞) = 0 h(r = R ra f t B ) = h ra f t B (1)
with h the height of the interface (h=0 for a flat interface), r the radial coordinate, h the derivative of h with respect to r, R ra f t B the radius of the raft and h ra f t B the height of the interface at the edge of the raft (see figure 1). We neglect the irregularities at the edge of the raft due to the presence of the particles, making the assumption that the shape of the interface around the raft is isotropic. We can thus write equation ( 1) in an axisymmetric configuration.

Oil

Water h raft B R raft B h r Fig. 1 Schematic side view of a granular raft. Deformation of the interface around a single raft B, described by two quantities: the radius of the raft R ra f t B , and the depth of the interface at the edge of the raft h ra f t B . The radial coordinates are defined from the center of the raft, while the origin of the vertical axis h corresponds to the height of the undisturbed flat interface.

After having measured R ra f t B and h ra f t B from a side view of raft B (without the presence of A), we solve equation (1) numerically, and obtain the expression of h(r) imposed by the presence of the raft B. Then, from the knowledge of the radial position L of the center of A with respect to the center of B, which is measured directly from the movie, we can deduce the expected vertical position of the granular raft A along its motion. Following this procedure, we can calculate the total speed V of a granular raft from its radial displacement. Implicitly here, we use the Nicolson linear approximation 15 ted in figure 2.b, for the two rafts, as a function of the nondimensionalized distance between the two centers. As a consequence, the curves have to be read from right to left, with an increase of V as the two rafts get closer, until it reaches a maximum speed. Then, the two rafts briefly slow down just before contact, because of a hydrodynamic coupling in the drag force: the liquid that separates the two rafts before collision has to be expelled, causing an increase in the drag. Qualitatively, we can already see in figure 2 that the sizes of the two rafts have an effect on their motion. Indeed a raft made of 30 particles will move at a larger speed than the raft made of 50 particles during their interaction. We now need a model to relate these speed profiles to the interacting forces acting on a granular raft.

F cap 1→1 = a(R part , ρ s , ...)K 1 L c (2)
with a a coefficient taking into account the particle and fluid properties ( 16 for more details), L the distance between the centers of the two particles, and K 1 the modified Bessel function of the second kind of order one. We generalize this formulation to the attraction of two granular rafts (designated by the letters A and B), as schematically represented in figure 3:

F cap A→B = f (n A , n B )a(R part , ρ s , ...)K 1 L c ( 3
)
where f is the function we want to determine experimentally, F cap A→B the force exerted by the raft A on B, and n A and n B the number of particles in each raft. To make such a generalization, we assume here that the wavelength and amplitude of the undulation of the edge of the raft are small enough so that they can be neglected at long range and the interface around the raft can be described as isotropic. We neglect here the granular nature of the raft, the edge of which may be roughened by individual particles.

The drag force acting on a raft also needs a careful treatment. For two identical spherical beads, F drag is expressed as a Stokes drag corrected by the mobility function G for two spheres in a single phase, in order to take into account the drainage of the liquid between the particles:

F drag = bV G -1 L R part (4) L l 2R part F drag A F cap B A A B

Oil Water

Fig. 3 Forces acting on a given axisymmetric granular raft. Schematic representation of the interaction between two granular rafts A and B (n A and n B respectively stand for the number of particles in each raft), at an oil-water interface. L is the distance between the centers of the two rafts, while l is the distance between the two nearest particles of each raft. The raft A (respectively B) is subjected to a capillary attraction F cap B→A and a viscous drag F drag A (respectively F cap A→B and F drag B ).

with V the speed of the particle, b = 6π µ o R part k a corrected Stokes drag coefficient on a sphere, where k accounts for the fact that the particle moves along an interface and as a consequence is immersed in two phases 24 . The mobility function G was tabulated by Batchelor 25 , and can be approximated by an interpolation formula 18 . At infinity, there is no coupling and G(+∞) = 1, whereas when the two particles are in contact, G(2) = 0. Again, we generalize equation ( 4) to the motion of a granular raft:

F drag A = g(n A )bV A G -1 l + 2R part R part (5)
where g is the function we want to determine experimentally, and l designates the distance between the edges of the two closest particles of the two rafts, as defined in figure 3. We approximate the hydrodynamic coupling between the two rafts by considering the drainage of the liquid only between the two closest particles of each interacting raft. Such a formulation of the drag is only valid for a sphere, but here we use it for a granular raft, which has the shape of a curved disk. As a consequence, the scaling law for the drag should also differ from the Stokes drag of a sphere, although we expect the general scaling laws for Stokes drag to be valid.

Keeping in mind the assumption lying under the scaling law of equation ( 5), we combine it with equation (3) and finally deduce an expression for the speed of a raft made of n A particles attracted by a second raft constituted by n B beads:

V A = f (n A , n B ) g(n A ) a b G l + 2R part R part K 1 L c ( 6
)
The ratio a/b is measured once and for all for each type of particle (given radius and given density, see table 1) thanks to a two-bead experiment, for which f (1, 1) = g(1) = 1 by definition. Then, f /g can be experimentally determined as a function of n A and n B , giving us information on the ratio of the two forces.

Experimental results

We plot in figure 4 the speed of the raft A as a function of the distance between the centers of the two rafts L, first with n A = 1 and an increasing n B , and then with n B = 60 and an in-creasing n A . The speed of the raft appears to increase both with n A and n B , as expected from our definition of the capillary force. This is due to the deformation of the interface that increases with n B . Similarly, F cap B→A is related to the weight and the size of the raft A, an increasing quantity with n A . Focusing on a single curve of figure 4, we recover the behavior described previously: an increasing speed as the two rafts are attracted towards one another, until a maximum where the speed then starts to decrease due to the drainage of the liquid between the two rafts.

As n A or n B increases, the curves are shifted to the right since L min = R ra f t A + R ra f t B , a value which increases as we add particles.

The black dotted curves of figure 4 represent the fit of equation (6) for each couple of rafts (n A , n B ). There is a good agreement between the experimental data and the theory for n A = 1 (figure 4.a), the fitting line being within the experimental noise. As we increase n A , the fit starts to drift from the measured speed, in particular regarding the position of its maximum. We discuss the limitations of our description in the last section of the paper and propose a number of hypothesis which may explain this result. Keeping in mind that for high n A , our model does not explain the entire dynamics, we observe the evolution of the fitting coefficient f /g as a function of n A and n B for different types of beads (see table 1), and a large variety of combinations of n A and n B (figure 5). All the data collapse onto the same straight line when represented as a function of n B √ n A , on approximately three decades, giving the following experimental result:

f (n A , n B ) g(n A ) = n B √ n A (7)
Equation (7) gives us an empirical relation between the ratio of the two forces and the number of particles in each raft. Even though both forces are functions of n A , n B only intervenes in the capillary force. This immediately proves that F cap A→B ∝ n B . Moreover, applying Newton's third law [26][27][28] , we obtain that this force has to be symmetrical in n A and n B and thus F cap A→B ∝ n A n B , implying also for the drag force g(n A ) ∝ √ n A . The far-field behavior of the capillary force F cap A→B ∝ n A n B follows a Coulomb or gravitational-like law, replacing the electric charge or the mass by the number of grains; since this force only depends on the number of grains, it indicates that the organization of the particles within such compact granular rafts plays no role in the interaction. This is consistent with the superposition approximation, where interface deformation at large distances can be considered as the sum of the individual meniscii of each particle. However we lack a physical understanding of these scalings and need to go back to the definition of the capillary force.

The drag and the capillary force

As explained earlier, the capillary force exerted by a raft, assuming it has an isotropic shape and the deformations of the interface are small, is proportional to the slope of the interface deformed by its presence. In the limit of small deformations we derive an analytical formula for this slope by solving a linearized version of equation ( 1):

     ∇ 2 h = h 2 c h(r → ∞) = 0 h(r = R ra f t ) = h ra f t (8)
which gives us the following expression:

dh dL = h ra f t / c K 0 (R ra f t / c ) K 1 L c . (9)
with K 0 the modified Bessel function of the second kind of order zero.

Equation ( 9) can be reconciled with the scaling F cap A→B ∝ n A if we show that the modified aspect ratio of a raft

h ra f t / c K 0 (R ra f t / c
) is linear with its number of particles n.

Experimental and numerical description of a raft shape

Experimentally, we measure h ra f t from side views of the granular raft (blue arrows, figure 6.a). We then measure the area of the raft A ra f t from top views (see video 2 of the supplementary information) and calculate the equivalent radius by R ra f t = A ra f t /π. The experimental modified aspect ratio is plotted in black diamonds in figure 6.b for different types of beads and raft sizes, as a function of the effective weight of a raft F vert , that we define as the difference between the weight of the particles and their buoyancy:

F vert (n) = nF vert (1) = n 4 3 πR 3 part g(ρ s - ρ o + ρ w 2 ) ( 10
)
with n the number of particles in the raft. In equation (10), both the weight and the buoyancy of the particles are considered linear with n, with the supplementary assumption that a bead is equally immersed in oil and water. However, we are aware that the generalized buoyancy of the raft, which takes into account surface tension effect, is not strictly additive, as discussed in previous collective sinking experiments 14 .

All the data collapse onto a single curve, and the modified aspect ratio appears linear with F vert (a power law regression on the experimental data gives an exponent of 1.01), which is by definition proportional to n. As a consequence, the modified aspect ratio is linear with the number of particles.

To this experimental evidence, we add a numerical confirmation for the morphology of a granular raft, using the code described in 21 . The raft is simulated as an axisymmetric continuous membrane of thickness 2R part , with a given volume fraction of a material of density ρ s , a complementary volume fraction of a material of density (ρ o + ρ w )/2, and a tension T directed tangentially to the raft surface. As long as the tension remains positive, a static solution is found, and the shape of both the granular raft and the interface beyond can be solved numerically.

In figure 6.a, the simulated interface is plotted on top of experimental visualisations of the corresponding granular raft, with a visually good agreement for the top three photos. For the last two granular rafts shown in figure 6.a, no numerical result is available: the simulations predict the sinking of the raft for such a deformation. This is probably due to the simplicity of the model used for the simulations, the limitations of which are discussed in our previous work 21 . As a consequence, the numerics can only give us information for part of the sizes accessible experimentally. However, they allow us to explore a far broader range of densities and radii of particles, and appear in that sense complementary to the experiments. The numerical measurement of the modified aspect ratio is represented in colored circles in figure 6. Here again, the data collapse onto the same curve when represented as a function of F vert , over more than four decades, with an exponent of 1.08 for a power fitting law. In conclusion, we recover the following result both experimentally and numerically:

h ra f t / c K 0 (R ra f t / c ) ∝ n (11)
Interestingly, in figure 6.b, we find that beyond a critical vertical force, the linear relation between the modified aspect ratio and F vert no longer holds (see the right portion of the blue curve). Such a curve corresponds to a situation where the raft reaches a maximal depth which corresponds to a balance between the hydrostatic pressure and the raft weight, and then extends indefinitely as more particles are added 21 . As a consequence, h ra f t is bounded, leading to the saturation of the capillary force. Since

K 0 (R ra f t / c ) ∼ π 2R ra f t e -R ra f t and F vert ∝ R 2
ra f t , it explains the vertical deviation in figure 6.b for large rafts.

Equation (11) can also be recovered theoretically in the limit of small deformations (see supplementary information for more details). All these results validate that:

F cap B→ A = n A n B F cap 1→ 1 . ( 12
)

The drag force on each raft

Finally, we need to explain how the drag force depends on the raft geometry:

F drag A = √ n A F drag 1 (13)
One could understand the drag dependency with n A by considering the equivalent sphere of radius R ra f t . For such a sphere, the Stokes drag is proportional to R ra f t , which for a large range of granular raft sizes evolves as √ n. This argument does not hold for the largest granular rafts (see for instance the last photo of figure 6), for which the high curvature of the raft makes R ra f t deviate from this law. However at first order, R ra f t ∝ √ n, leading to the expected drag on the granular raft. The proportionality between the drag force and the raft radius is still quite surprising, since the shape of a granular raft is not a perfect sphere, but more a curved disk.

Discussion

A number of hypotheses underlying the approach described in this paper are no longer valid for a granular raft. They give us some clues to understand the limits of our model. In particular we observe an increasing discrepancy between the fit and the experimental data measured for the velocities when the rafts become large (the two highest velocities of figure 4.b for n A = 20 and 40). The two-sphere canonic example we use to derive equations (3) and ( 5) strongly relies on the assumption of small deformations around spherical beads (small in relation to the capillary length). With granular rafts, we can reach a vertical deformation of the interface up to almost half the capillary length. Moreover, the morphology of a raft, a two-dimensional circular monolayer of beads, cannot be approximated as a spherical shape, and one could therefore expect the drag to be drastically different. Finally, the viscous hypothesis used for the calculation of the drag can be questioned for the bigger rafts. In this case, the maximum speed reached during their motion is almost a hundred times higher than for two single beads, leading to a Reynolds number of order one. All these differences may change the amplitude of the forces through f and g and even the form of these functions.

We also assumed that the capillary force is proportional to K 1 (L/ c ), a result that comes from the linearized solution of equation (1). However, for some experiments, when two dense granular rafts interact, as in figure 7, equation (1) becomes highly non-linear, and the approximate linear solution is no longer sufficient to account for the curvature of the interface. Figure 7 illustrates this limitation of our model: on top of the experimental visualisation of the interface when two granular rafts are almost in contact, we plot the interfacial shape we would get if the deformations could be superimposed, deduced from the experimental measurement of the shape of the interface around a single granular raft of the same size. The superposition principle seems to account rather well for the curvature on the left and right side of each raft, but clearly overestimates the depth of the interface between the two rafts, leading to an error on the calculation of the capillary force. This could be part of the explanation for the disagreement between the fit and the measured speed in figure 4 for the motion of the two biggest rafts.

Another important feature that we did not take into account is the lateral extension of a granular raft. As already emphasized, the derivation of equation ( 2) for the capillary force between two spheres at an interface relies on various hypotheses, including the size of the particles: R part c . However, the largest rafts, as il- Comparison between a side view of two granular rafts during their motion towards one another (background photo) and the expected interfacial shape obtained if the assumption of small deformations was valid (blue dotted curve). Each of the rafts is made of 50 ceramic spherical beads (density ρ s = 3,800 kg.m -3 , radius a single particle R part = 0.45 mm), in an axisymmetric configuration. The photo is taken 0.1 s before contact. Scale bar: 5 mm. lustrated in figure 7 for instance, can exhibit a diameter of the order of the capillary length. In this situation equation (2) should not be valid anymore. Here again, this could explain the discrepancy between our experimental data and the corresponding fitting curve for large rafts.

Despite all these limits, our model still accounts for the major part of our data, proving that in a first approximation, a granular raft does behave as a heavy membrane which is not deformed by the interaction with other rafts.

Conclusion

We have investigated the interaction dynamics of a two-body system formed after the clustering of beads into granular rafts at an oil-water interface. In particular we described the capillary forces generated by each granular raft, as well as the drag due to their individual motion. Because the deformation around a granular raft strongly depends on its size (depth and width), the capillary interaction can surpass by several orders of magnitude the forces created by individual beads, with a very strong dependence on the number of particles. A generalization of the two-particle model is found to account very satisfactorily for the motion of the granular rafts, and only starts to deviate from the experimental data for objects of dimensions close to the capillary length. Despite the clear non-linearity of the equations describing the shape of the interface for such large rafts, the capillary force appears linear both with the size of the attracting and the attracted raft. The drag, on the other hand, increases linearly with the radius of the considered raft (or in other terms, as the square root of the number of particles).

In this paper we have completely neglected the intrinsic discrete nature of a granular raft by modeling it as a twodimensional membrane. However, in some cases the granular aspect of the raft may play an important role. We have observed that rafts made of smaller particles (diameter < 400 µm) show a weaker cohesion between the grains within the raft, which may lead to erosion and internal motions for a fast enough raft displacement. Further work is thus required to characterize this phenomenon and investigate the role of shape changes on the overall interaction of granular objects at a liquid interface. Discrete aspects of the raft, as well as elasticity, could be also included in future developments of the model.

Fig. 2

 2 Fig.2Aggregation of two granular rafts. a) Time lapse of the motion of two granular rafts at an oil-water interface, filmed from the side (see supplementary information for the movies, taken from above and from the side). Each raft is made of ceramic spherical beads (density ρ s = 3,800 kg.m -3 , radius of a single particle R part = 0.45 mm), in an axisymmetric configuration, with 50 particles for the raft on the left (raft A), and 30 particles for the one on the right (raft B). The final image corresponds to the exact moment of contact between the two rafts. Time between two images: 0.13 s. Scale bar: 10 mm. b) Speed of the same two rafts as a function of the distance between the centers of the two rafts (made dimensionless using the radius of a bead R part ). Red circles: speed of raft B; black squares: speed of raft A. The errors bars are of the order of the thickness of the curves.

Fig. 4 Fig. 5

 45 Fig.4Agreement between experimental results and theoretical model. Speed of the raft A for different number of particles in A and B, as a function of the distance between the centers of the two rafts (made dimensionless using the radius of a bead R part ), fitted by the theoretical prediction of equation (6) (black dotted curves). For all these experiments, the particles are ceramic spherical beads (density ρ s = 3,800 kg.m -3 , radius of a single particle R part = 0.45 mm). a) Speed of a single particle A for different sizes of the other raft, n B = 1, 5, 10, 30, 60. b) Speed of a raft of increasing size (n A = 1, 5, 10, 20, 40), attracted by a raft of fixed size (n B = 60). The errors bars are of the order of the thickness of the curves.

1 Fig. 6

 16 Fig. 6 Morphology and weight of a granular raft. Relation between the deformation of the interface and the weight of the granular raft, for a broad variety of densities, radii and numbers of particles. a) Side view of a granular raft for an increasing number of ceramic particles (density ρ s = 4,800 kg.m -3 , radius of a single particle R part = 0.35 mm). In the top three photos, the red line represents the shape of the interface calculated by the numerical simulation for a membrane of a radius similar to the background photo (dotted line: oil-water interface, solid line: granular raft). The blue arrows indicate the measurement of R ra f t and h ra f t . Scale bars: 3 mm in each photo. b) Modified aspect ratio h ra f t / c K 0 (R ra f t / c ) as a function of the effective weight of the raft F vert , for experiments (black diamonds), and numerical simulations (colored circles). The solid line has a slope of 1. The vertical errors bars for the experiments are of the order of the size of the individual points, and are thus not displayed for clarity reasons.

Fig. 7

 7 Fig. 7 Example of the failure of the linear superposition assumption.Comparison between a side view of two granular rafts during their motion towards one another (background photo) and the expected interfacial shape obtained if the assumption of small deformations was valid (blue dotted curve). Each of the rafts is made of 50 ceramic spherical beads (density ρ s = 3,800 kg.m -3 , radius a single particle R part = 0.45 mm), in an axisymmetric configuration. The photo is taken 0.1 s before contact. Scale bar: 5 mm.

Table 1

 1 Characteristics of the beads used in experiments

	Type of particle	ρ s (kg.m -3 )	R part (mm)	n max
	Plastic	1,420	2	16
	Plastic	1,420	2.5	7
	Ceramic	3,800	0.35	160
	Ceramic	3,800	0.45	60
	Ceramic	4,800	0.35	60
	Ceramic	4,800	0.45	25

1-8 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Dynamics of two interacting granular rafts3.1 Theoretical modelWe study the influence of the number of particles on the respective motion of two granular rafts. To that end, we use the theoretical framework derived for the dynamics of two identical spheres interacting at an interface, in the limit of small deformations16 , and see how the amplitude of the forces is affected by the size of the raft.The force between two objects at an interface depends only on the shape of the interface around them, and is obtained by integrating the capillary forces around their contact line. This shape results from how the particles deform the interface, which in the case presented here is the consequence of the combined effect of their weight, their wettability and the induced buoyancy. In the case of small spheres (R part c ) where the shape of the interface is isotropic around each bead and its deformations are small, the calculation leads to an energy of interaction equal to the prod-uct between the effective weight of one particle and the vertical deformation created by the other. For two identical spheres, the force exerted by one bead on the other is then obtained by deriving the interaction energy, leading to:

| 1-8 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Conflicts of interest

There are no conflicts to declare.

Notes and references


E A Bufi

P Cinnella

Efficient uncertainty quantification of turbulent flows through supersonic ORC nozzle blades

Keywords: method of characteristics, uncertainty quantification, dense gas, organic rankine cycle"

This work aims at assessing different Uncertainty Quantification (UQ) methodologies for the stochastic analysis and robust design of Organic Rankine Cycle (ORC) turbines under multiple uncertainties. Precisely, we investigate the capability of several state-of-the art UQ methods to efficiently and accurately compute the average and standard deviation of the aerodynamic performance of supersonic ORC turbine expanders, whose geometry is preliminarily designed by means of a generalized Method Of Characteristics (MOC). Stochastic solutions provided by the adaptive Simplex Stochastic Collocation method, a Kriging-based response surface method, and a second-order accurate Method of Moments are compared to a reference solution obtained by running a full-factorial Probabilistic Collocation Method (PCM). The computational cost required to estimate the average adiabatic efficiency, Mach number and pressure coefficient, as well as their standard deviations, to within a given tolerance level is compared, and conclusions are drawn about the more suitable method for the robust design of ORC turbines.

Introduction

The design under multiple uncertainties is a technique called Robust Optimization (RO), which allows achieving designs with stable performances under random variations of the design parameters. Here we consider an Organic Rankine Cycle (ORC) supersonic impulse expander working under stochastic operating and geometrical conditions. In order to carry out a robust design, which takes into account the non-deterministic nature of inputs, the sensitivity of the 2D turbulent flow through a geometry representative of a typical supersonic ORC turbine nozzle cascade is preliminarily investigated via different uncertainty quantification (UQ) algorithms coupled with an in-house finite-volume dense gas flow solver. The supersonic nozzle is designed by means of an in-house software based on the method of characteristics extended to dense gases [START_REF] Elio | Sensitivity of supersonic ORC turbine injector designs to fluctuating operating conditions[END_REF], which are peculiar of ORC applications. Dense gases are molecularly complex gases with a strongly non-ideal thermodynamic behavior in a thermodynamic region above the liquid/vapor saturation curve, where the speed of sound exhibit a behavior opposite to that of standard gases when the fluid undergoes isentropic transformations. Unfortunately, experimental comparisons are not possible since for dense gas flows no experimental data of any kind are available in the literature up to now [START_REF] Congedo | Numerical investigation of dense-gas effects in turbomachinery[END_REF]. This is a well-known difficulty in the ORC community. Attempts to build experimental setups to provide data are currently underway at several universities. ORC systems are subject to multiple uncertainties, like fluctuating operating conditions, geometric tolerances and ill-known fluid properties, which should be taken into account at an early stage of the design process (see, e.g. [START_REF] Elio | Sensitivity of supersonic ORC turbine injector designs to fluctuating operating conditions[END_REF]). Problems involving a high number of uncertain parameters suffer from the so-call "curse of dimensionality" problem [START_REF] De Baar | Efficient uncertainty quantification with gradient-enhanced kriging: Applications in fsi[END_REF][START_REF] Congedo | Efficient robust optimization techniques for uncertain dense gas flows[END_REF], since the number of code runs required to approximate the statistical moments of the probability density functions associated to the output quantities of interest increases exponentially with the number of parameters. Besides, due to the complexity of the geometry and computational cost associated with the ORC geometries and working fluids, no advanced UQ method has been applied so far to these turbines and little work has been done on uncertainty quantification in turbomachinery in general [START_REF] Zou | Stochastic analysis of a radial-inflow turbine in the presence of parametric uncertainties[END_REF]. A few applications of non-intrusive, sparse grid Generalized Polynomial Chaos methods to simple ORC turbine simulations exist [START_REF] Congedo | Efficient robust optimization techniques for uncertain dense gas flows[END_REF][START_REF] Hercus | Robust shape optimization of uncertain dense gas flows through a plane turbine cascade[END_REF]. For realistic applications such the ORC turbines considered in this paper, the overall computational cost of an UQ calculation becomes prohibitive. Computational cost is especially crucial when the aim is using robust design techniques in an industrial context. To overcome these limitations, the selection of efficient sampling techniques of the parameter space is of vital importance. Then, the aim of the paper is to assess and compare very different UQ techniques, in terms of accuracy and efficiency, for a realistic ORC nozzle configuration and to suggest useful guidelines to engineers and designers for selecting an efficient UQ method for this kind of problems.

Dense gas solver

In the numerical simulations presented in this work, the viscous governing equations are discretized using a cell-centered finite volume scheme for structured multi-block meshes of third-order accuracy, which allows the computation of fluids modelled by the Peng-Robinson-Stryjek-Vera cubic equation of state in order to take into account the real gas effects [START_REF] Cinnella | Numerical solver for dense gas flows[END_REF]. The closure of the RANS equations system is achieved by implementing the Spalart-Allmaras turbulence model. The scheme is obtained by correcting the dispersive error of the classical second-order accurate Jameson's scheme [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF]. To preserve the high accuracy on non-Cartesian grids, the numerical fluxes are constructed by using weighted discretization formulas, which take into account the mesh deformations [START_REF] Rezgui | Third-order accurate finite volume schemes for Euler computations on curvilinear meshes[END_REF]. This ensures third-order accuracy on moderately distorted meshes and second-order accuracy at least on highly deformed mesh. The equations are then integrated in time using a four-stage Runge-Kutta scheme. Local time stepping, implicit residual smoothing and multi-grid acceleration are used in order to drive the solution to the steady state. The accuracy of the numerical solver has been already demonstrated in previous works [START_REF] Cinnella | Numerical solver for dense gas flows[END_REF][START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF][START_REF] Rezgui | Third-order accurate finite volume schemes for Euler computations on curvilinear meshes[END_REF][START_REF] Cinnella | Aerodynamic Performance of Transonic Bethe-Zel'dovich-Thompson Flows past an Airfoil[END_REF].

Uncertainty quantification models

In this Section a short description of the UQ models used in this work is provided. Two non-intrusive models (Simplex Stochastic Collocation (SSC) and Bayesian-Kriging (BK) methods) and one deterministic model (Second Order Second Moment (SOSM) method) have been coupled with the dense gas solver in order to evaluate the response of the ORC injector to the stochastic variations of the inputs. As described in the next Section, the results are compared with those provided by a full-factorial analysis performed with a Probabilistic Collocation Method (PCM).

Probabilistic Collocation Method (PCM)

In the PCM [START_REF] Loeven | Probabilistic collocation: an efficient non-intrusive approach for arbitrarily distributed parametric uncertainties[END_REF] the multi-dimensional random space parameters is discretized in a full-factorial way by means of multiple one-dimensional tensor products. If we consider the generic partial differential equation system (where x are the spatial variables and is the random solution depending on time, space and the random parameter ω), the solution is decoupled in a deterministic part and in a stochastic part , where is the input random variable. The random solution is expanded by means of a Lagrange polynomial chaos according the equation ( 1):

(1)

The number of terms Np for the expansion increases exponentially as function of the number of random variables nξ according to the relation , where P is the degree polynomial. This feature is typical of the full-factorial sampling and introduces the "curse of dimensionality" problem, which has a dramatic impact on the computational costs. The statistics (mean and variance) are evaluated by means of a Galerkin projection of the solution on the polynomial basis, with the collocation points calculated as the nodes of the Gaussian quadrature integration performed by implementing the Golub-Welsch algorithm [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF].

Simplex Stochastic Collocation with Extremum Diminishing (SSC-ED) method

To alleviate the high computational cost associated to a full-factorial sampling, a more efficient technique is proposed in [START_REF] Witteveen | Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces[END_REF]. The hypercubic space parameters is discretized in a unstructured way under the form of simplexes by means of an adaptive Delaunay triangulation. A random sampling is performed in order to avoid clustering of points and high order polynomial interpolation allows to reconstruct the random solution for each simplex. This method is suitable for dealing with discontinuities, such as shocks, thanks to the local smoother based on a Local Extremum Conserving criterion. A super-linear convergence is reached during the simplex refinement process. The statistics are evaluated using the definition of expectation by means of a Monte-Carlo integration on the polynomial interpolation. The refinement process for a generic simplex is carried out according to a local error, which depends exponentially from the local polynomial degree and the number of random variables. The adaptive random sampling of the space parameters is expected to improve the efficiency, leading to a lower overall number of deterministic calculations than the full-factorial sampling.

Bayesian Kriging (B-K)

The third approach we compare is a Kriging-like surrogate model of the random output, derived by using a Bayesian framework [START_REF] De Baar | Improvements to gradient-enhanced Kriging using a Bayesian interpretation[END_REF]. We are interested in the value of a quantity of interest (QoI) as function of certain random variables at n locations contained in the hypercube space parameters. At m of these locations the QoI is observed and known (by running a deterministic calculation) while at the remaining mn points it needs to be estimated. Following a Bayesian inference approach, the probability distribution of the QoI can be calculated as the posterior distribution resulting from the combination of a multivariate normal prior distribution and likelihood function. In the latter, the covariance matrix takes into account the spatial correlation among the random variables by means of a geometrical correlation range. Then, for each location , the solution is approximated by a Gaussian distribution with a certain mean and variance. The error of the surrogate model respect to the observation is provided at the random locations i by the standard deviation, as

, where is the covariance matrix.

Mean Value Second Order Second Moment method (MVSOSM)

The Mean Value Second Order Second Moment method (MVSOSM or simply SOSM) is a so-called deterministic UQ model which computes the statistics of the random QoI by means of a Taylor expansion around the mean value. In this work the mean and variance are obtained with a second and a first order truncated expansion respectively, according to equations ( 2)-( 3):

(2)

(3)

The cross derivatives are neglected in the Taylor expansion, i.e. the interaction among the random parameters are not considered here. If a centered scheme is used for all the derivatives, the number of deterministic calculations Np varies linearly with the number of random variables , as . The method requires a relatively low number of samples and is expected to be cheap in terms of computational costs. An additional deterministic calculation is required in order to evaluate the QoI at the mean value of the random inputs.

Results

The objective of this work is to assess which is the more suitable UQ model for performing a sensitivity analysis of the ORC injector designed with the MOC under stochastic input parameters. After discussing the deterministic numerical results at nominal conditions, we analyse the stochastic solutions provided by the different UQ methods under investigation.

Nominal operating conditions results

First of all, the performances at nominal operating conditions have been evaluated by means of the dense gas solver. A deterministic run requires about 10 hours of CPU time on a single processor machine. The Reynolds-Averaged-Navier-Stokes (RANS) equations where solved numerically on a structured Cshaped mesh with cells. The inlet total thermodynamic conditions, periodicity in peripheral direction and a supersonic outlet are imposed as boundary conditions. The working fluid is penta-fluoropropane R245fa, largely used in ORC applications, and the operating conditions are chosen close the saturation curves in the supercritical region, characterized by significant dense gas effects. The injector geometry has been preliminarily designed by using the MOC, with the inlet total reduced pressure and temperature , the target massflow G and the Mach number Me at exit section (see table 1) as input conditions.

Table 1. Design parameters of the ORC expander and nominal thermodynamic conditions Figure 1 shows the results at nominal operating conditions in terms of Mach number distribution. The flow expands up to a Mach number of 2 in the exit cross section area, as expected from the MOC design. Afterwards it continues to accelerate guided by the suction side wall and then is decelerated by a weak oblique shock departing from the trailing edge. A viscous wake is visible at the exit and an interaction between the oblique shock and the boundary layer can be noticed close to the suction side wall. These phenomena are expected to be cause of high variability in the flow and, then, a UQ analysis is carried out.

Full-factorial sensitivity analysis with PCM

Five uncertain parameters are considered in this analysis: inlet total reduced pressure and temperature (pr, Tr), with a variability of 5% around the nominal point, and inlet flow incidence angle β, blade thickness ε, stagger angle θ, with a variation of 1% around the design values. First of all, a full-factorial analysis via PCM has been carried out. By considering a second order polynomial chaos expansion and the five uncertain parameters, the number of deterministic calculations is , requiring 3000 hours of total calculation time. Second-order polynomial chaos was found to provide a good compromise between accuracy and computational cost [START_REF] Hercus | Robust shape optimization of uncertain dense gas flows through a plane turbine cascade[END_REF]. In figure 2a the distribution of the Mach number variance shows the presence of regions with high uncertainty close to the wall of the divergent part and to the trailing edge, which can be addressed to the viscous effects. In a second step the random space parameters cardinality has been lowered in order to take into account only the most influential random parameters by means of an analysis of variance (ANOVA) [START_REF] Tang | Global sensitivity analysis for stochastic collocation expansion[END_REF]. This allows reducing the number of the deterministic calculations and, then, saving computation time. Figure 2b compares the contributions of the five uncertain parameters to the global variance of the Mach number and shows a great influence of the blade thickness, the inlet total pressure and temperature, while the incidence and stagger angles have a negligible effect. As consequence, the uncertain parameter space is reduced to these three uncertain variables and the four UQ models are applied to the reduced space in order to compare their performances in terms of accuracy and calculation time requirements.

UQ models comparison

The reduced parameter space has been sampled by using the UQ models described in the previous Section. The results are shown in figures 2c-2d under the form of pressure coefficient mean (along the blade wall) and variance (on the nozzle blade axis). The mean distributions provided by the four models show similar results, with variations below 2% respect to the PCM, while the analysis of the variance diagram reveals that SOSM results are not very accurate, due to the first order approximation of variance and neglecting cross interaction. The B-K model, trained with a set of 20 observations, under-estimates the variance in the leading edge region respect to SSC-ED and PCM which, instead, show very close results. However, the three models all agree about the higher variance on the trailing edge, where the viscous effects are important. The analysis of the UQ models performances (see table 2) in terms of calculation time provides that the SSC-ED model is the more expensive, with a computational cost higher and a number of deterministic calculations required Np slightly lower than PCM, but with higher accuracy. Indeed, a third order polynomial reconstruction is obtained on 83% of the total simplexes after refinement, against the 2th global order of PCM. In order to compare these two UQ models based on a polynomial expansion with the B-K surrogate model, a B-K error based criterion has been used. The nodes of the PCM full-factorial and SSC-ED multi-dimensional grids of the parameter space have been used as training set for the B-K model and the cell center values as prediction set. In this way, a criterion to compare the efficiency of the discretization of the space parameters is provided. The maximum error value of the B-K predictions with the PCM and SSC-ED grids is compared with that provided by the B-K for 20 randomly distributed observations (see figure 3). The SSC-ED provides the lower error, with a value slightly below the B-K one. This analysis shows that the random adaptive unstructured SSC-ED grid refinement criterion is very accurate for this application, but its computational cost remains relatively high. Finally, the SOSM method is the cheapest one, with only 7 deterministic calculations required, however its low accuracy for calculating second moments makes it suitable only for fast preliminary analyses.

Conclusions

In this work a comparison of different UQ methods, coupled with a numerical solver suitable for dense gas calculation for ORC applications, has been carried out in order to choose the more efficient and accurate one. The PCM, SSC-ED and B-K provide similar quantitative results. PCM and SSC-ED use approximately the same number of samples and then lead to a similar computational cost, however the SSC-ED is more accurate thanks to the higher order polynomial approximation and better distribution of points in the space parameters grid. Also PCM and B-K lead to almost the same results, but B-K is 70% faster for a slightly higher error. Then, for the present application, the B-K results to be a good compromise between low computational cost requirement and accuracy. The SOSM, with cross-derivatives neglected and first order expansion of the variance, shows low accuracy but gives useful information for engineering design purposes with a very low computational cost.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Mach number distribution for a viscous deterministic calculation at nominal operating conditions

Fig. 3 .

 3 Fig. 3. Comparison of the maximum Kriging error among the PCM, SSC-ED and B-K grids

Table 2 .

 2 Performance analysis of the UQ models