Learning nonlinear constitutive models in ﬁnite strain electromechanics with Gaussian process predictors

This paper introduces a metamodelling technique that employs gradient-enhanced Gaussian process regression (GPR) to emulate diverse internal energy densities based on the deformation gradient tensor F and electric displacement ﬁeld D 0 . The approach integrates principal invariants as inputs for the surrogate internal energy density, enforcing physical constraints like material frame indifference and symmetry. This technique enables accurate interpolation of energy and its derivatives, including the ﬁrst Piola-Kirchhoff stress tensor and material electric ﬁeld. The method ensures stress and electric ﬁeld-free conditions at the origin, which is challenging with regression-based methods like neural networks. The paper highlights that using invariants of the dual potential of internal energy density, i.e., the free energy density dependent on the material electric ﬁeld E 0 , is inappropriate. The saddle point nature of the latter contrasts with the convexity of the internal energy density, creating challenges for GPR or Gradient Enhanced GPR models using invariants of F and E 0 (free energy-based GPR), compared to those involving F and D 0 (internal energy-based GPR). Numerical examples within a 3D Finite Element framework assess surrogate model accuracy across challenging scenarios, comparing displacement and stress ﬁelds with ground-truth analytical models. Cases include extreme twisting and electrically induced wrinkles, demonstrating practical applicability and robustness of the proposed approach.


Introduction
Electro Active Polymers (EAPs) have emerged as a category of intelligent materials capable of undergoing substantial changes in shape in response to electrical stimuli [9,32,[53][54][55].Among these, dielectric elastomers are particularly noteworthy due to their remarkable actuation capabilities, encompassing attributes such as light weight, rapid response times, flexibility, and low stiffness properties.Notably, these materials can undergo electrically induced substantial strains (with reported area expansions of up to 1962% [37] as observed in research conducted at Harvard's Suo Lab).Their potential is exceptionally promising, with applications spanning bio-inspired robotics [8,11,29,36,47], humanoid robotics, and advanced prosthetics [9,40,68], as well as implications for tissue regeneration [48].
The realm of nonlinear continuum mechanics has reached an advanced stage of development, encompassing the variational formulation, finite element implementation, and principles related to the constitutive modeling of EAPs [6,7,16,17,67].In the context of the latter, the reversible constitutive model for dielectric elastomers is encapsulated within the free energy density, contingent upon the deformation gradient tensor F and the material electric field E 0 .Complementary to this potential, which exhibits a saddle point nature, is the internal energy density, contingent upon the deformation gradient tensor and the electric displacement field D 0 .Building upon this foundation, researchers in [22,51,52] introduced an extension of the concept of polyconvexity, originally from the field of hyperelasticity [2,[12][13][14]64], into this coupled electromechanical scenario.This novel definition of polyconvexity played a pivotal role in establishing the existence of minimizers in this context for the first time [65], serving as a sufficient condition for the extension of the rank-one convexity criterion within electromechanics.
In spite of the considerable inherent potential exhibited by EAPs, a primary limitation arises from their demand for a substantial electric field magnitude to induce significant deformation, rendering them prone to electromechanical instabilities or electrical breakdown [4,57,70].To mitigate the requirement for high voltage operation, some researchers propose the adoption of composite materials as the basis for EAPs.These composites often amalgamate a low-stiffness, low-permittivity elastomer matrix with stiffer, higher-permittivity inclusions distributed randomly as fibers or particles.Experimental investigations have evidenced a substantial enhancement in the coupled electromechanical performance of electroactive composites, thereby reducing the voltage prerequisites for actuation.A noteworthy development in recent years pertains to rank-one laminate composite dielectric elastomers [15,43,44].
The determination of the macroscopic constitutive response of the composite material hinges upon the specific type of microstructure under consideration.In the case of laminated composite materials, the homogenization challenge at the microstructure level is governed by a system of nonlinear equations that implicitly establish the microstructural parameters with respect to the macroscopic strain gradient tensor and the electric displacement field [15,21,58].In the case of more intricate microstructures, such as randomly distributed inclusions embedded within an elastomeric matrix, the determination of the macroscopic constitutive response of the composite material necessitates the utilization of computationally intensive homogenization techniques.However, these methods come with a limitation.EAPs exhibit nonlinear behavior, leading to a non-linear dependency of their macroscopic response on macroscopic deformations and electric or magnetic fields.Essentially, this signifies that a boundary value problem must be solved at the micro level, considering suitable boundary conditions, for every stress and electric/magnetic field combination [63].
In the effort to mitigate the high computational demands associated with methods like computational homogenization, recent developments in the realm of nonlinear continuum mechanics have witnessed the emergence of Machine Learning algorithms.These methods enable the generation of diverse constitutive models through the utilization of data gathered from experimental tests or in-silico (computational)-based computations.This paper aims to explore a specific type of Machine Learning technique, namely Gaussian Process Regression (GPR), to showcase its viability in approximating the constitutive model of analytical constitutive models in electromechanics, particularly within a specific category of composites known as rank-one laminates.By applying rank-n homogenization principles, it becomes feasible to derive the homogenized or effective response of rank-one laminates in an almost analytical manner, without necessitating computational homogenization.Utilizing in-silico data generated based on this model, the objective is to create surrogate models capable of replicating the behavior of the aforementioned types of constitutive models.This initial focus pertains to simpler scenarios, with the intention of demonstrating the feasibility and accuracy of the employed GPR technique.This approach sets the stage for addressing more complex composite cases in near future works, eliminating the need for computational homogenization and facilitating a computationally efficient evaluation of their effective behaviour [63].
Artificial Neural Networks (ANNs) have been employed for learning or discovering constitutive models based on data generated either in silico or in physical laboratories, as demonstrated in studies several [10,20,31,38,39].The work by Klein [31] represents a pioneering effort in the successful application of ANNs for uncovering constitutive laws in nonlinear electromechanics.Additionally, Gaussian Process Regression (GPR) has gained traction and found application in the development of data-based constitutive models for moderate strains in soft tissue applications, as shown by Aggarwal et al. [27].A distinctive feature GPR compared to the ANN approach lies in its inherent probabilistic nature.GPR allows for the specification of prior knowledge, the generation of a distribution encompassing potential predictive functions, and the direct calculation of prediction uncertainties [5,41,56].Moreover, GPR or Kriging offers control over the degree of interpolation between known points through the specification of noise in the correlation function [20].
Kriging [33,46,62] predicts Gaussian random field values using observed data from a finite set of points, finding applications in geostatistics, numerical code approximation, global optimization, and machine learning [28,46,56,61].It employs Gaussian distributions to define a joint distribution based on observations and predictions, utilizing spatially correlated covariance to weigh observation importance.The joint distribution conditions on observed data, yielding a prediction distribution characterized by mean and covariance, facilitating sampling for predictions [49].This emulator type has grown popular due to its nonlinear function capture and statistical output [28,45], yielding confidence intervals and adaptive metamodel refinement strategies.
This paper proposes a gradient-enhanced Gaussian process regression metamodelling technique for emulating internal energy densities characterizing soft/flexible EAP behavior.The method enforces physical constraints upfront by incorporating principal invariants as inputs.Gradient Kriging excels in precise interpolation of energy, first Piola-Kirchhoff stress tensor.Derivative incorporation reduces sampling points while maintaining accuracy.In contrast to neural networks, Kriging's interpolatory nature precisely matches stress tensors at sample points, ensuring stress-free origin compliance.
The structure of this paper unfolds as follows: In Sect.2, we establish the foundational concepts by introducing the essential elements of nonlinear continuum electromechanics, emphasizing constitutive modeling.Moving forward, Sect. 3 provides a comprehensive and self-contained overview of Gaussian Process Regression (GPR) or Kriging.Proceeding to Sect. 4, we undertake the calibration of Krigingbased surrogate models by employing synthetic data derived from well-established ground truth internal energy densities.Lastly, Sect. 5 exemplifies the practical application of these surrogate models within a 3D Finite Element computational framework.A thorough assessment is conducted to gauge the precision of these models across diverse and demanding scenarios, juxtaposing displacement and stress fields against their corresponding ground-truth analytical model counterparts.Notation Throughout this paper, A :

Differential governing equations in finite strain electromechanics
Let B 0 denote a subset of three-dimensional Euclidean space R 3 , representing the initial, undeformed state of an Electro Active Polymer (EAP) material.We postulate the existence of an injective function φ, which uniquely maps each point X in the material configuration B 0 to a corresponding point x in the deformed, spatial configuration B ∈ R 3 .This mapping relationship is mathematically expressed as x = φ(X) (as illustrated in Fig. 1).Associated with the mapping φ, we define the deformation gradient tensor The behavior of the EAP represented by B 0 is governed by the ensuing coupled boundary value problem: where the equations on the left correspond to the purely mechanical physics and those on the right hand side, with the electrostatics equations.In (1), DIV(•) signifies the divergence with respect to the material coordinates X ∈ B 0 , while f 0 represents the force applied per unit volume B 0 .Dirichlet boundary conditions for the field φ are imposed on ∂ φ B 0 , and t 0 represents a force per unit undeformed area, being N the outward normal at X ∈ ∂ t B 0 .Furthermore, on the right hand side of (1) ρ 0 represents an electric charge per unit undeformed volume B 0 .Dirichlet boundary conditions are prescribed on ∂ ϕ B 0 for the field ϕ, and ω 0 represents an electric charge per unit undeformed area ∂ ω B 0 , being N the outward normal at X ∈ ∂ ω B 0 .For both coupled physical problems, the boundaries where Dirichlet and Neumann boundary conditions are prescribed satisfy the following Finally, P and D 0 symbolize the first Piola-Kirchhoff stress tensor and the material electric displacement field, respectively.These tensors are interlinked with the deformation gradient tensor F and the material electric field E 0 by means of an appropriate constitutive law, as described in Sect.2.2.

The internal energy density in electromechanics
The constitutive model of the undeformed solid B 0 is encapsulated in the internal energy density per unit underformed volume, denoted as Taking the derivative of the internal energy density with respect to both F and D 0 gives rise to the first Piola-Kirchhoff stress tensor P and the material electric field E 0 as defined in Eq. ( 1) The internal energy density e(F, D 0 ) is required to adhere to the principle of objectivity, also known as material frame indifference.This entails its invariance with respect to rotations Q ∈ SO(3) applied to the spatial configuration, as follows Moreover, the internal energy density must conform to the material symmetry group G ⊆ O(3), a defining factor in determining the isotropic or anisotropic attributes of the underlying material.This requirement can be succinctly expressed in mathematical terms as follows Furthermore, the internal energy density e(F, D 0 ), along with the first Piola-Kirchhoff stress tensor P and the material electric field E 0 , must all vanish when no deformations are present, i.e.

e(F, D
The conditions in Eqs. ( 5), (6), and (7) embody essential physical criteria.Alongside these, there is a requisite for the internal energy density function to satisfy pertinent mathematical criteria.Specifically, the internal energy density function conventionally adheres to mathematical constraints rooted in the concept of convexity.One of the simplest conditions is that of convexity of e(F, D 0 ), that is which requires positive semi-definiteness of the Hessian operator [H e ], defined as However, convexity away from the origin (i.e., F ≈ I, D 0 ≈ 0) is not practically suitable, as it doesn't encompass realistic material behaviors like buckling [2].An alternative mathematical constraint is the quasiconvexity of e(F, D 0 ) [3].Unfortunately, quasiconvexity is a nonlocal condition that is challenging, and even infeasible, to verify.An implied requirement of quasiconvexity is that of generalized rankone convexity of e(F, D 0 ).A generalized rank-one convex energy density satisfies Remark 1 Notice that the vector V ⊥ in (10) is orthogonal to V .The reason for this choice has its roots in the analysis of the hyperbolicity of the system of PDEs in (10) in the dynamic context.In this case, it is customary to express the fields φ and D 0 as a perturbation with respect to equilibrium states φ eq and D eq 0 , respectively, by means of the addition of travelling wave functions as where V represents the polarisation vector of the travelling wave and c the associated speed of propagation of the perturbation with amplitudes u and V ⊥ .Introduction of the ansatz for D 0 into the Gauss's law in Eq. (1) reveals that and therefore, V ⊥ must be orthogonal to V .
Condition (10) is known as the Legendre-Hadamard condition or ellipticity of e(F, D 0 ).It is associated with the propagation of traveling plane waves in the material, defined by a vector V and speed c [50].Importantly, the existence of real wave speeds ab initio for the specific governing Eq. in (1) is assured when the electromechanical acoustic tensor Q ac is positive definite, with with A sufficient and localized condition aligned with the rankone condition in (10) is the polyconvexity of e.The internal energy density is considered polyconvex [2,22] if a convex and lower semicontinuous function W : where H and J represent the co-factor and determinant of F, and with the vector d a vector in the spatial configuration, being the three of them defined as Polyconvexity of the internal energy density entails the satisfaction of the following inequality

Invariant-based electromechanics
A simple manner to accommodate the principle of objectivity or material frame indifference and the requirement of material symmetry is through the dependence of the internal energy density function e(F, D 0 ) with respect to invariants of the right Cauchy-Green deformation gradient tensor C = F T F and D 0 .Let I = {I 1 , I 2 , . . ., I n }, represent the n objective invariants required to characterise a given material symmetry group G.Then, it is possible to express the strain energy density e(F, D 0 ) equivalently as Application of the chain rule into Eq.( 4) yields the first Piola-Kirchhoff stress tensor P and the material electric displacement field D 0 in terms of the derivatives of U (I) as

Isotropic electromechanics
For the case of isotropy, the invariants required to characterise this material symmetry group, and the first derivatives of the latter with respect to F and D 0 (featuring in the definition of P and E 0 in ( 19)) are Inserting the expressions in ( 20) into ( 19) yields the following expression for the first Piola-Kirchhoff stress tensor P and for E 0

Transversely isotropic elastromechanics
In the context of transverse isotropy, a preferred direction N emerges, perpendicular to the material's plane of isotropy, imparting anisotropic characteristics.Our focus centers on the material symmetry group D ∞h [25], where the structural tensor takes the form N ⊗ N.This group is distinct from C ∞ , also present in transversely isotropic materials, characterized by the structural vector N and encompassing the potential for piezoelectricity.The D ∞h group, beyond the invariants (20), is distinguished by three additional invariants, detailed below (22) In this case, the first Piola Kirchhoff stress tensor P and the electric field E 0 adopt the following expressions

Application to rank-one laminates
Section 2.3 presented the case of phenomenological internal energy densities using principal invariants.In composite materials, computing effective strain energy density requires homogenization.This section focuses on rank-one laminates, composed of two constituents perpendicular to N. Rank-n homogenization theory [15] relates macroscopic F, D0 1 to microscopic F a , F b , D a 0 , D b 0 as where indices a and b differentiate the constituents and c a and c b denote their respective volume fractions, with where α ∈ R 3 and β ∈ R 2 represent the mechanical and electric amplitude vectors, respectively, which need to be determined.Furthermore, being T 1 and T 2 any two perpendicular vectors to N, and Remark 2 Notice that in Eq. in (25), although it might seem a priori not intuitive, the definition of F a in terms of c b and vice-versa (and also for D a 0 and D b 0 ) is necessary in order to comply with Eq. ( 24).The first of these two equations entails 1 We could continue using the notation F and D 0 for the macroscopic deformation gradient tensor and electric displacement, respectively.However, we prefer to adhere to the traditional convention in homogenisation theory [4,24,59], where macroscopic fields are generally referred to as average fields, for which the bar symbol on top of these fields can be used, hence the notation F, D0 .
which is clearly satisfied, and the second which is also satisfied.The same derivations and conclusions are obtained when considering D a 0 and D b 0 as in (25).
The determination of α and β can be done by postulating the effective energy e( F, D0 ) as The stationary conditions of ê with respect to α and β yield which represent two nonlinear vector equations from which {α, β} can be obtained.Finally, computation of {α, β} permits to obtain the effective first Piola-Kirchhoff stress tensor and electric field as

Gaussian process predictors
In the realm of computer experiments, metamodelling or surrogate modelling entails substituting a resource-intensive model or simulator U = M (I) with a computationally efficient emulator M (I).Both the simulator and emulator share the same input space D I ⊆ R n and output space D U ⊆ R.
In our context, M represents the response of an actual internal energy density U , dependent on principal invariants I (as discussed in Sect.2.3).Thus, we replace the common input field x with I and the output y with U .As the internal energy U is scalar, the theoretical developments presented in this paper exclusively pertain to scalar outputs.Our approach employs Kriging models [30,60], also known as Gaussian Process (GP) modelling.Succinctly, the key components of this method will be detailed in Sects.3.1-3.3.

Gaussian process based prediction
GP modelling assumes that the output U = M (I) is characterised by where, g(I) • β signifies the prior mean of the Gaussian Process (GP), representing a linear regression model over a specific functional basis g i , i = 1, . . ., p ∈ L 2 (DI , R).The subsequent component, denoted as Z (I), characterizes a GP with a zero mean, a constant variance σ 2 U , and a stationary autocovariance function, defined as follows where R is a symmetric positive definite autocorrelation function, and θ , the vector of hyperparameters.In this work we employ the Gaussian kernel for the definition of R, defined as The construction of a Kriging model consist of the twostage framework described in the upcoming Sects.3.2 and 3.3.

The conditional distribution of the prediction
The Bayesian prediction methodology assumes that observations gathered in the vector U = U (I (1) ), U (I (2) ), . . ., U (I (m) )

T
The observed values, including any unobserved U (I), constitute a realization of a random vector adhering to a parametric joint distribution.This section seeks to derive a stochastic prediction for this unobserved quantity by harnessing this statistical interdependence.The Gaussian assumption for Z (I) in Eq. ( 31) and the linear regression model's nature enable the inference that the observation vector U is also Gaussian, characterized by (34) being G and R the regression and correlation matrices, defined as and Likewise, a new random vector, encompassing the observed set U alongside any unobserved value U (I), follows a joint Gaussian distribution, presented as where g(I) is the vector of regressors evaluated at I and r(I) is the vector of cross-correlations between the observations and prediction given by Assuming that the autocovariance function given by Eq. ( 32) is known, the conditional distribution of the prediction Û (I) = U (I)|U is governed by the Best Linear Unbiased Predictor (BLUP) theorem [62].As per BLUP, the unobserved quantity U (I) = M (I) in the prior model of Eq. ( 31) follows a Gaussian distribution, represented by the Gaussian random variable Û with a specific mean: and variance where is the generalised least-squares estimate of the underlying regression problem, and For those readers unfamiliar with the previous derivations, we suggest consulting the comprehensive treatment offered in Reference [18], which provides a thorough elucidation of the foundational mathematical principles and methodologies involved.

Joint maximum likelihood estimation of the GP parameters
In the previous sections, we operated under the assumption of a known autocovariance function.However, the specifics of the correlation functions R(I, I , θ ) and the variance value σ 2 U are typically not known in advance.In this study, we predefine the correlation function type (specifically we make use of Gaussian kernels in (33) [18]), and the determination of hyperparameters θ and variance σ 2 U is achieved using the observation dataset via the technique of maximum likelihood estimation (MLE).The outcome of this process yields the empirical best linear unbiased predictors (EBLUP) [62].The estimation of GP parameters involves solving the following minimization problem where L (U|β, σ 2 U , θ ) is the opposite log-likelihood of the observations U with respect to its multivariate normal distribution given by The MLE of β and σ 2 U are obtained from the first order optimality conditions of L (U|β, σ 2 U , θ ), namely from which the following optimal values can be obtained Substituting β * (θ ) and σ 2 U * (θ ) into the log-likelihood function (44) enables it to be re-written as where the reduced likelihood function has been introduced as This entails that the minimisation problem in Eq. ( 43) is equivalent to Unfortunately, an analytical solution for the optimal hyperparameters θ ∈ R n is unavailable.Instead, a numerical minimization approach is typically employed.In our research, we employ the box-min algorithm [69].

Gradient-enhanced Gaussian-process based prediction
In addition to function observations, leveraging output derivatives concerning input variables is possible, aiming to enhance predictor accuracy.This gives rise to what is termed Gradient Enhanced Kriging in the literature [23,35], in contrast to the conventional Kriging detailed in Sects.3.1-3.3.To establish a gradient-enhanced predictor, the observation vector is extended to encompass derivatives of the strain energy density U concerning its input variables I, resulting in: 1) , . . ., U (m) , ∂ I U (1) , . . ., where To interpolate both the variable and its gradient at any unobserved location, the extension of the correlation matrix R is necessary to incorporate the correlation between the variable and its gradient, formulated as where R UU is the correlation matrix presented in (36) for the non-gradient case.R UU includes the partial derivatives of R according to 1) , I (1) , θ ) . . .∂ I (m) R(I (1) , m) , I (1) , θ ) . . .∂ I (m) R(I (m) , given , . . ., ∂R(I (i) , I ( j) , θ) The submatrix R U U contains the second derivatives 1) R(I (1) , I (1) , θ) . . .∂ 2 I (1) I (m) R(I (1) , I (m) , θ) . . . . . . . . .m) , I (1) , θ) . . .∂ 2 where Similarly the vector of cross-correlations between the observations and the prediction is extended as follows r(I) = R(I, I (1) , θ ), . . ., R(I, I (m) , θ ), ∂ I (1) R(I, I (1) , θ), . . ., ∂ I (m) R(I, I (m) , θ) T .(57) Once these adaptations have been made, the revised definitions for the respective quantities can be inserted into the descriptions provided in Sects.3.2 and 3.3.To begin, let us revisit the mean prediction and the variance with Finally, the optimal hyperparameters are achieved by minimizing the log-likelihood function.
Remark 3 The conceptual framework presented in Sects. 2 and 3 has been constructed around the internal energy density e(F, D 0 ).Yet, intertwined with this formulation, lies the ability to define its corresponding dual, referred to as the free energy density, symbolized as Ψ (F, E).This duality is established through the subsequent Legendre transformation (Fig. 2) The free energy density Ψ (F, E 0 ) imposes distinct convexity constraints compared to its dual counterpart e(F, D 0 ).As a consequence, in the proximity of F ≈ I and E 0 ≈ 0, it assumes the character of a saddle point function, exhibiting convexity with respect to F while displaying concavity concerning E 0 .This divergence in convexity/concavity attributes in the context of both mechanics and electro physics can introduce challenges in the application of Kriging or Gradient Enhanced Kriging interpolation models that rely on invariants of F and E 0 (i.e., free energy-based Kriging), as opposed to those formulated in terms of F and D 0 (i.e., internal energy-based Kriging).Notably, our observations underscore that an internal energy-focused approach yields markedly superior outcomes compared to the utilization of the free energy density methodology.

Remark 4
Applying Kriging and Gradient Enhanced Kriging techniques (discussed in Sects.3.1 to 3.4) can be criticized due to the need to include strain energy values (U ) at each observation or training point.This constraint limits their suitability for datasets from physical experiments, unlike those from in-silico or numerical sources, the paper's focus.Quantifying energy measurements poses challenges in such cases.Yet, the Gradient Enhanced Kriging excels in adaptability, accommodating cases with a single observation (F = I) where U often equals zero.Here, derivative information at this point, combined with derivatives from other points, can be seamlessly integrated.This tailored approach is elaborated in Appendix C.

Derivatives of strain energy density for Gradient Enhanced Kriging
As detailed in Sect.3.4, the gradient-enhanced Kriging method incorporates both the internal energy density U and its derivatives concerning the invariants I = I 1 , I 2 , . . ., I n .
In cases involving isotropy or transverse isotropy within material symmetry groups, coupled with a principal invariant approach (see Sect. 2.3), addressing the derivatives of U with respect to I becomes imperative.While obtaining these derivatives for analytical energies like those derived from a Mooney-Rivlin model is straightforward, intricate internal energy densities arising from complex homogenization techniques in composites (e.g., rank-one laminates in Sect.2.4) may lack readily available derivatives.In such scenarios, deriving these derivatives from the first Piola-Kirchhoff stress tensor P and the electric field E 0 can be undertaken using conventional linear algebra principles.To facilitate this, let's revisit the invariant-based expressions for P and E 0 in (19), conveniently restated below where Let us introduce now the following notation where P ∈ R 9 and Ûi ∈ R 9 represent the vectorised expressions for both P and U i .This entails that A can be conveniently written in terms of W i as In (67), W i can be understood as the linear independent vectors of a basis, whilst ∂ I i U represent the coordinates of A along the vectors W i , (i = {1, . . ., n}).As standard in basic courses of linear algebra, given A, the coordinates ∂ I i U can be obtained through projection of the latter over the n vector of the basis, which yields the following linear system of equations Remark 5 By examining the algebraic system of equations presented in (68), it becomes feasible to ascertain the conditions under which a solution for {∂ I 1 U , . . ., ∂ I n U }, representing the components of A with respect to the basis {W 1 , . . ., W n }, can be derived.Notably, the solvability of (68) (i.e., the linearity independence of {W 1 , . . ., W n }) hinges on the determinant of the system, which must not equal zero.Ill-conditioning in the equation system (68) can stem from several factors: identical principal stretches of deformation (found in both isotropic and transverse isotropic models), alignment of principal deformation directions with the preferred direction of transverse isotropy, or alignment of D 0 with one of the principal directions of F. To rectify this numerical issue, we propose a perturbation approach, introducing slight variations to the identical principal stretches and a misalignment of the coinciding principal direction with the preferred direction.This ensures solvability of (68).

Noise regularisation
In cases of substantial training data and the incorporation of derivative information into the training strategy (e.g., in the context of gradient-enhanced Kriging), the correlation matrix R defined in Eq. ( 52) can become ill-conditioned.To mitigate this issue, a customary practice is to introduce regularization by augmenting the correlation matrix with a diagonal matrix as follows: [19,49,66] While our paper primarily highlights the interpolation properties of this technique, we consistently employ sufficiently small values of ε 1 and ε 2 to mitigate potential challenges.It is noteworthy, as elucidated in Remark 3, that Kriging and its gradient counterpart can achieve interpolation when R remains unregularized, specifically for ε 1 = ε 2 = 0.However, when ε 1 = 0 and ε 2 = 0, and in the extreme scenario of both parameters assuming larger values, Kriging transitions from an interpolation technique to a regression technique.Thus enabling to filter noisy data.
To illustrate this technique we explore the performance of Gradient-Enhanced Kriging in a regression context, particularly when confronted with severely ill-conditioned correlation matrices arising from noise-contaminated data.To elucidate this aspect, we employ two designated training samples, denoted as: -Unperturbed sample training sample devoid from noise in the output variables, including the values of the energy Ψ (F) = U (I 1 , I 2 , I 3 , I 5 ) and its derivatives where the ground-truth constitutive model from which these data have been generated in-silico corrresponds with the Mooney-Rivlin/ideal dielectric model described in Appendix A. -Noisy sample this training sample has been obtained by perturbing the deterministic sample according to: with In both datasets, Fig. 3 illustrates the performance of interpolation and regression based approaches.In the case of the unperturbed sample (see Fig. 3a and b), Kriging prefectly reproduces the training data points (represented by circles).Conversely, in the noisy sample, Kriging strives to replicate the perturbed and irregular data to the greatest extent possible.Discrepancies observed at certain points, resulting in the ill-conditioning of the correlation matrix and subsequently the loss of interpolation properties.Notably, it is evident that the condition number of the matrix R experiences a substantial increase when dealing with the noisy sample, as illustrated in Fig. 3c.This observation aligns with expectations and raises concerns regarding the predictive accuracy of Kriging between training points, potentially leading to undesired oscillations.
Alternatively, we have explored a regression-based methodology, as detailed in [19,49].In this context, the regularization parameters {ε 1 , ε 2 } are treated as supplementary hyperparameters.Consequently, both sets of hyperparameters, namely {θ 1 , θ 2 , θ 3 } and {ε 1 , ε 2 }, are optimized through the minimization of the reduced likelihood function ψ( θ) where the augmented set of hyperparameters is defined as Applying this approach to only the the noisy sample yields the outcomes depicted in Fig. 3.The values of {ε 1 , ε 2 } are determined to strike a balance between the interpolation and regression properties of the Kriging response.Naturally, the response does not precisely match the noisy data, thereby avoiding the introduction of undesirable oscillations caused by data perturbations (see Fig. 3f).

Design of Experiments
In this section, we present the procedure used for generating synthetic data, utilizing a diverse set of ground truth constitutive models.The internal energy densities and material parameters for these models can be found in A. To acquire the dataset, we adhere to the procedure outlined in [34], extended to the coupled context of electromechanics.The deformation gradient tensor F is parameterized via a chosen set of deviatoric directions, amplitudes, and Jacobians (J , i.e., the determinant of F).The process of generating sample points for deviatoric directions, amplitudes, and Jacobians is elucidated in Algorithm 1.Similarly, the electric displacement D 0 is also parametrised in terms of unitary directions and amplitudes.Concerning the deviatoric directions for F, denoted as V F we formulate them using a spherical parametrization in R 5 , precisely representing these directions using five pertinent angular measures ( ) within this 5-dimensional space.For the directions employed for the parametrisation of D 0 , denoted as V D 0 , these are created using a spherical parametrization in R 3 , using as angular measures (θ, ψ) ∈ [0, 2π ] × [0, π], namely Once the sample is generated following Algorithm 1, the reconstruction of the deformation gradient tensor F and of D 0 becomes possible at each of the sampling points.This reconstruction process is demonstrated in Algorithm 2, where represents the basis for symmetric and traceless tensors (refer to Appendix B for details on ).
Construct the directions, V F , using an extended Spherical parametrisation in R 5 -detailed in (73); 6: Construct the directions, V D 0 , using Spherical parametrisationdetailed in (73); 7: Evaluate the deformation gradient tensors, F and D 0 -detailed in Algorithm 2; 8: Evaluate the invariants, energy U , P and E 0 ; 9: Get derivatives ∂ I i U according to Eq. ( 68); Algorithm 2 Pseudo-code for construction of the set of deformation gradient tensors and electric displacement fields for j = 1 : n J do 3: for l = 1 : n t D 0 do 6: end for 8: end for 9: end for 10: end for

Calibration and Validation
The synthetic data, generated as per Sect.4.1, calibrates Kriging and Gradient Enhanced Kriging surrogates, following principles in Sect.3. To assess surrogate accuracy at non-observation points, generated evaluation points mirror the procedure in Sect.4.1.These points test the model performance but are not part of calibration.For validation, a substantial validation set of 10, 000 data points is used.This density ensures assessment of the smaller calibration set's accuracy.Validation comprehensively evaluates the surrogate model's performance, verifying its reliability and generalizability.
The calibration and validation process has been carried across a diverse range of constitutive models.These include: (a) Mooney-Rivlin/ideal dielectric model (MR/ID); (b) Arruda-Boyce/ideal dielectric (AB/ID) (see Reference [1]); (c) Gent/ideal dielectric (Gent/ID); (d) Quadratic Mooney-Rivlin/ideal dielectric (QMR/ID); (e) Yeoh/ideal dielectric (Yeoh/ID); (f) Rank-one laminate composite (ROL).Specific expressions for strain energy densities and material parameters are available in A. For each model, 2 training datasets are generated, each containing N = {45, 100} training points.Kriging and Gradient Enhanced Kriging models are calibrated for all 6 ground truth models within each training set.Results include mean squared error (R 2 ( P) and R 2 (E 0 )) for first Piola-Kirchhoff stress tensor P and E 0 , and values of Ê P and Ê E 0 , defined below are presented in Table 1 (for N = 45 training points) and 2 (for N = 100 training points).In Eq. (74), || A|| denotes the Frobenius norm of A, n is the number of experiments, P An i and P Kr i represent the analytical and Kriging-predicted first Piola-Kirchhoff stress tensors, respectively.Similarly, E An 0 i and E Kr 0 i represent the analytical and Kriging-predicted material electric field, respectively.The findings of the analysis, as presented in Tables 1 and 2, offer insights into the performance of Kriging and Gradient Enhanced Kriging techniques.In both tables, the achieved   R 2 ( P) and R 2 (E 0 ) values are notably high, approaching unity, signifying an impressive level of accuracy in predicting the first Piola-irchhoff stress tensor.This is true for all the models except for the rank-one laminate composite material, where the performance of the ordinary Kriging approach is extremely poor.Furthermore, under the consideration of the alternative metric, specifically Ê P , Ê E 0 , Gradient Enhanced Kriging demonstrates a significantly superior accuracy, consistently yielding values approximately an order of magnitude smaller compared to the Kriging counterpart.For comprehensive understanding, Fig. 4 depicts the evolution of the Ê P metric for both conventional and Gradient Enhanced Kriging methodologies.This illustration pertains to two specific constitutive models considered in Tables 1 and  2, namely the Mooney-Rivlin/ideal dielectric model and the rank-one laminate composite model.Notably, as the number of training points increases, the Gradient Enhanced technique adeptly diminishes this metric, substantiating its efficacy.On the contrary, the ordinary Kriging method is incapable of decreasing the metric Ê P as the number of infill points increases.
These observations emphasize the distinct advantage of adopting the Gradient Enhanced technique, as it facilitates precise predictions of the first Piola-Kirchhoff stress tensor P and of the material electric field E 0 even when operating with an considerably small number of training points.This characteristic positions Gradient Enhanced Kriging as an exceedingly expedient and efficacious alternative in comparison to the conventional Kriging methodology.

Remark 6
Appendix A contains the material parameters used in each of the constitutive models (MR/ID, AB/ID, Gent/ID, QMR/ID, Yeoh/ID, ROL) considered for calibration of their respective Kriging and Gradient Krigind predictors.Notice that the values of these material parameters do not correspond with those of typical dielectric elastomers such as the VHB 4910 by 3 M.It is importante to emphasize that our Kriging and Gradient Kriging predictors are flexible to deal with any realistic value of material constants, and in particular, those typical of the popular VHB 4910.These materials exhibit a large disparity between the values of mechanical constants and electrical constants.For instance, the shear modulus μ and electric permittivity ε of VHB 4910 material [26] take values of approximately μ ≈ 10 3 −10 (Pa) and ε ≈ 10 −11 − 10 −12 (F/m).
The enormous gap between the material constants of both physics (i.e.mechanics and electric physics) can ultimately pose challenges for the accurate calibration of the Kriging predictor (i.e.yielding ill-conditioning of the correlation matrix R (36)) or of any other type of machine learning technique.In order to remedy this, instead of considering the data generated by the model with such material constants, and in particular, the first Piola-Kirchhoff stress tensor P, the material electric field E 0 and the material electric displacement D 0 , we can alternatively perform the calibration with their dimensionless counterparts, P, E 0 and D 0 (notice that the deformation gradient tensor F is already dimensionless), respectively, defined as Remark 7 With regards to the rank-one laminate material, in our previous publications (see Reference [42]), we demonstrated that whenever each of the phases a, b of the rank-one laminate comply with the polyconvexity condition (15), and therefore with the ellipticity condition (10), the solvability of α and β in ( 29) is always guaranteed.This entails that at microscopic level there is no apparent difficulty.However, the homogenised response of phases which are elliptic individually does not necessarily inherit this desirable property, and can indeed exhibit loss of ellipticity.We have not discarded this situation for the calibration of the Kriging and Gradient Kriging predictors, and in fact, some points within the data generated violate the ellipticity condition (at the macroscopic level).Despite this, the predictors can handle these situations.

Numerical three-dimensional examples
The analysis in Sect. 4 strongly supports the superiority of gradient enhanced Kriging over its energy-only counterpart, which lacks first derivatives.Inspired by these promising results, the primary objective of this section entails the seamless integration of gradient enhanced Kriging models into an in-house Finite Element computational framework.This assimilation endeavors to establish the accuracy and efficacy of these metamodels through meticulous juxtaposition with the Finite Element solutions provided by their respective ground truth counterparts.Specifically, this evaluation embraces intricate and exacting scenarios including complex bending and wrinkling, thus furnishing a robust appraisal of the metamodels' performance within demanding contexts.

Electrically induced bending example: isotropic ground truth model
The inaugural exemplification within the Finite Element domain revolves around a cantilever beam configuration, as illustrated in Fig. 5.The geometric attributes and boundary conditions underpinning this scenario are succinctly elucidated in Fig. 5. Pertaining to the discretization framework, tri-quadratic Q2 Finite Elements have been judiciously employed to effectuate the interpolation of the displacement field.
In this illustrative instance, we examine a Mooney-Rivlin/ideal dielectric model (as expressed in Eq. ( 76)) as the ground truth internal energy density.We adopt the material parameters specified in Table 3. Upon subjecting the system to a voltage differential V √ ε/μ 1 = 0.5 (see Table 3 for the value of μ 1 ), the ensuing deformation is depicted in Fig. 6a and b for the ground truth Mooney-Rivlin model and its gradient enhanced Kriging counterpart, respectively.Evident congruity emerges between both figures.This congruence is also manifest in the contour plot of σ 13 , where σ denotes the Cauchy stress tensor, i.e., σ = J −1 P F T .These consistent parallels serve to affirm the precision and robustness of the gradient enhanced Kriging models, signifying their potential in effectively capturing the intricate behavior of the underlying physical systems.3 for the value of μ 1 ) Fig. 7 Complex electrically induced actuation.Geometry and electrical boundary conditions.Displacements fixed at X 1 = 0. {a, b, c} = {120, 10, 1} (mm).Electrodes highlighted with colour green (two regions in the lowest surface across the thickness and one region on the top surface) and red colour (one region in the mid surface across the thickness)

Complex electrically induced bending example: rank one laminate ground truth model
The second example considers the same cantilever beam as in the preceding section, subjected to a more complex set of boundary conditions for the electric potential ϕ.This can be seen in Fig. 7. Pertaining to the discretization framework, tri-quadratic Q2 Finite Elements have been judiciously employed to effectuate the interpolation of the displacement field.
In this instance, we consider a more complex constitutive model in comparison to the precedent section, where the selection encompassed an isotropic ground truth model.Our present investigation is directed towards a rank-one laminate composite ground truth model.The homogenized internal energy governing this model is encapsulated within Eq. (81), while the associated material parameters are cataloged in Table 9. Upon the imposition of a voltage gradient denoted as V ε a /μ a 1 = 2.5 across the electrodes (see Table 9 for the value of μ a 1 and ε a ), the intricate phenomenon of electrically induced bending is explored for both the ground truth and Gradient Enhanced Kriging models.The outcomes of this analysis are presented in Fig. 8, effectively showcasing the marked concordance evident in the electrically induced deformations, as well as the alignment in stress distributions, between the two models.
For a more comprehensive evaluation of the electrically induced deformation in the context of both models (the isotropic ground truth model and its corresponding counterpart developed using the Gradient Enhanced Kriging approach), an enhanced comparative perspective is available in Fig. 9.This representation serves to underscore the notable concurrence observed between the deformation predictions of the two models.

Electrically induced wrinkles: rank one laminate ground truth model
The last example considers the geometry and boundary conditions shown in Fig. 10.This example has been previously analysed in other works by the authors [43,44].The squared plate is completely fixed in its borders.The voltage is grounded at the maximum value of coordinate X 3 whilst a surface charge ω 0 = 220/ μ a 1 ε a (Q•mm −2 ) (see Table 9 for the value of μ a 1 and ε a ) is applied at the minimum value of coordinate X 3 .The Finite Element discretisation considers Q2 (tri-quadratic) hexahedral Finite Elements with 80 × 80 × 2 elements in X 1 , X 2 and X 3 directions.
The primary objective of this illustrative instance is to assess the precision of the Gradient Enhanced Kriging model within scenarios characterized by the emergence of wrinkles induced through electrical stimuli.In a specific context, our focus centers on employing the rank-one laminate model, The pertinent material parameters integral to this model can be found in Table 9. Upon the application of an electric charge denoted as ω 0 , the progression of electrically induced wrinkles is portrayed in Fig. 11, these predictions being furnished by the Gradient Enhanced Kriging model across a spectrum of escalating ω 0 values.In addition, this figure represent the evolution of the displacement along Z or X 3 direction of points A and B (see Fig. 11) predicted by both ground truth and Gradient Enhanced Kriging (Emula-  9 for the value of μ a 1 and ε a ).Volumetric force of value 9.8×10 −2 (N/m 3 ) action in along X 3 axis (in the positive direction) tor) models, showing a clear similarity between the paths of both models.
Crucially, Fig. 12 offers a side-by-side comparison of the wrinkles projected by the rank-one laminate composite model, serving as the veritable benchmark, and its concomitant representation through the Gradient Enhanced Kriging methodology.Evidently, the congruence between the patterns of electrically induced wrinkling is remarkable, further corroborated by the similarity observed in the distribution of stress patterns as depicted in the contour plots of both models.9 for the value of μ a 1 ) for = 0.1.a Rank-one laminate composite ground truth model; b transversely isotropic Gradient Enhanced Kriging model

Concluding Remarks
This manuscript introduced an innovative metamodelling technique that leverages gradient-enhanced Gaussian Process Regression or Kriging to emulate a diverse range of internal energy densities.The methodology seamlessly incorporates principal invariants as input variables for the surrogate internal energy density, thereby enforcing crucial physical constraints such as material frame indifference and symmetry.This advancement has facilitated precise interpolation not only of energy values, but also their derivatives, including the first Piola-Kirchhoff stress tensor and material electric field.Furthermore, it ensures stress and electric field-free conditions at the origin, a challenge typically encountered when employing regression-based methodologies like neural networks.The research has indicated the inadequacy of utilizing invariants derived from the dual potential of the internal energy density, particularly the free energy density.The inherent saddle point nature of the latter diverges from the convex nature of the internal energy density, engendering complexities for models based on GPR or Gradient Enhanced GPR that rely on invariants of F and E 0 (free energy-based GPR).This is contrasted with models formulated using F and D 0 (internal energy-based GPR).
Numerical examples within a 3D Finite Element framework have been thoughtfully incorporated, rigorously assessing the accuracy of surrogate models across intricate scenarios.The comprehensive analysis juxtaposing displacement and stress fields with ground-truth analytical models encompasses situations involving extreme bending and electrically induced wrinkles, thus showcasing the utility and accuracy of the proposed approach.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A.1 Mooney-Rivlin/ideal dielectric elastomer model
The internal energy density for this model is: and the material parameters used are (Table 4)

A.3 Gent/ideal dielectric elastomer model
The internal energy density for this model is: and the specific values for the material parameters used are (Table 5)

A.4 Yeoh/ideal dielectric elastomer model
The strain energy density for this model is: and the material parameters used are (Table 6):

A.6 Transversely Isotropic/ideal dielectric elastomer model
The internal energy of the this model is: The material parameters used for this model are (Table 8):

A.7 Rank-One Laminate
We consider Mooney-Rivlin and ideal dielectric internal energy densities for the individual phases a and b within   25).The material parameters used for this composite material are found below (Table 9): where the mechanical and electrical contrasts f m and f e are defined as and with the vector of lamination N spherically parametrised in terms of α and β according to N = cos θ sin ψ sin θ sin ψ cos ψ . (85) and the use of repeated indices implies summation.The tensor product is denoted by ⊗ and the second order identity tensor by I.The tensor cross product operation between two artibrary second order tensor A and B entails [ A B] I J = E I P Q E J RS A P R B QS .Furthermore, E represents the third-order alternating tensor.The full and special orthogonal groups in R 3 are represented as O(3) = {A ∈ R 3×3 , | A T A = I} and SO(3) = {A ∈ R 3×3 , | A T A = I, det A = 1}, respectively and the set of invertible second order tensors with positive determinant is denoted by GL

Fig. 1
Fig. 1 mapping of material quantities to the spatial quantities

Fig. 2 a
Fig. 2 a Convex nature of the strain energy density e(F, D 0 ) and b saddle point nature of the free energy density Ψ (F, E 0 ) in the vicinity of F ≈ I and E 0 ≈ 0

Fig. 3
Fig. 3 Performance of regression-based Gradient-Enhanced Kriging.a, b interpolation using an unperturbed training sample; c, d interpolation using a perturbed training sample; e, f Regression using a perturbed training sample

Fig. 4 Algorithm 1
Fig. 4 Evolution of metric Ê P for both ordinary and Gradient with the number of training of infill points for: a Mooney-Rivlin/ideal dielectric (MR/ID); b rank-one laminate composite model (ROL)

Fig. 8 Fig. 9
Fig. 8 Complex electrically induced actuation.Contour plot distribution of σ 13 /μ a 1 for various of V μ a 1 /ε a (see Table 9 for the value of μ a 1 and ε a ): a Rank-one laminate composite ground truth model; b Transversely isotropic Gradient Enhanced Kriging counterpart

Fig. 11
Fig. 11 Electrically induced wrinkles.Wrinkling patterns for various values of surface charge ω 0 , being the load factor.Top row: results predicted by the transversely isotropic Grandient Enhanced Kriging model, calibrated against rank-one laminate composite ground truth Acknowledgements R. Ortigosa and J. Martínez-Frutos acknowledge the support of grant PID2022-141957OA-C22 funded by MCIN/AEI/10.13039/501100011033and by "RDF A way of making Europe".The authors also acknowledge the support provided by the Autonomous Community of the Region of Murcia, Spain through the programme for the development of scientific and technical research by competitive groups (21996/PI/22), included in the Regional Program for the Promotion of Scientific and Technical Research of Fundacion Seneca -Agencia de Ciencia y Tecnologia de la Region de Murcia.The fourth author also acknowledges the financial support received through the European Training Network Protechtion (Project ID: 764636) and of the UK Defense, Science and Technology Laboratory.Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

5 }
being the effective strain energy Ψ F, D0 e F, D0 = arg minα,β {ê F, D0 , α, β }; ê F, D0 , α, β = c a e a F a F, α , D a 0 D0 , β + c b e b F b F, α , D b 0 D0 , β ,(82)with e a F a F, α , D a 0 D0 , β = U a I a 1 , I a 2 , I a 3 , I a 5 ; e b F b F, α , D b 0 D0 , β = U b I b 1 , I b 2 , I b 3 , and {I b 1 , I b 2 , I b 3 , I b 5 } represent the principal invariants of F a , F b , D a 0 and D b 0 related to the macroscopic deformation gradient tensor F and D0 through Eq. (

Table 1 R
2( P), R 2 (E 0 ), Ê P and Ê E0 for all six models for number of training points N = 45, for both Kriging and Gradient Enhanced Kriging

Table 4
Material parameters used with the quadratic Mooney-Rivlin/ideal dielectric model

Table 5
Material parameters used with the gent/ideal dielectric model

Table 6
Material parameters used with the Yeoh/ideal dielectric model

Table 7
Material parameters used with the Arruda-Boyce/ideal dielectric model

Table 9
Material parameters used with the rank one laminate model