Reduced and All-At-Once Approaches for Model Calibration and Discovery in Computational Solid Mechanics Open Access

Within the framework of continuum solid mechanics, initial boundary value problems involving partial differential equations (PDEs) and appropriate initial/boundary conditions are solved to determine the fields of interest, including mechanical, thermal, magnetic, and/or electric fields. To be solvable, the PDEs need to be completed by closure models known as material models (or constitutive equations) that describe the relationships between kinematic and kinetic quantities (and possibly additional variables), i.e., strain and stress, or temperature gradient and heat flux, etc., for the materials under consideration. Constitutive equations are often of algebraic, ordinary differential (ODE), or differential-algebraic (DAE) form—or they can even be defined by PDEs, whereby the space-discretized form of the PDEs leads to similar mathematical structures as the aforementioned forms. Models of hyperelasticity lead to algebraic equations, whereas models of viscoelasticity or viscoplasticity often imply ODEs—although they can also be written as integro-differential equations. Rate-independent plasticity, which is based on a yield condition, often delivers a DAE-system, but, in the variational context, can also be formulated as a differential inclusion problem. Nonlocal models such as gradient-based elasticity, plasticity, or damage are described by PDEs. Once again, in the space discretized form the main structure of the equations resulting from PDEs is similar to that stemming from ODEs or DAEs; hence, for simplicity, in this paper, we limit ourselves to local constitutive laws.

Since constitutive models depend on various parameters, one of the fundamental issues in the theory of materials is the identification (or calibration) of these parameters on the basis of given experimental data. A number of identification approaches have been proposed, among which those based on least squares (LS) and the finite element method (FEM) have been especially prominent due to their high flexibility. Further methods have been proposed in the more recent past, including so-called full-field approaches such as the virtual fields method (VFM), and stochastic and machine learning approaches, including those based on physics-informed neural networks (PINNs). A detailed compilation of the related references is provided in the following sections.

Even more recently, approaches that bring identification one step further have been proposed and are drawing significant attention. They advocate a new paradigm, denoted as (material) model discovery, in which the structure of the material model itself is estimated at the same time as its unknown parameters, again based on the given experimental data. A possible strategy consists of forming a library of candidate basis functions, which are used in a model to approximate unknown functions describing the material behavior. In this way, the model discovery problem is reformulated as a parameter estimation problem, but with a typically much higher dimension of the parameter space. It is also possible to apply VFM, stochastic methods, and machine learning to model discovery.

These developments call for a new unified perspective that is able to cover both traditional and novel parameter estimation and model discovery approaches. The purpose of this paper is, on the one hand, to compile an overview of the available methods along with their formalization and, on the other hand, to provide a unified perspective. In contrast to existing attempts to unify parameter estimation methods in mechanics [1,2], our perspective is based on the all-at-once approach [3], which so far has been mainly applied in the inverse problems community. Adopting concepts discussed in this community, we distinguish between all-at-once and reduced approaches. The all-at-once approach employs a weighted combination of a model-based and a data-based objective function, and we demonstrate that both the reduced approach and the class of VFMs can be recovered as limit cases. With this general framework, we are able to structure a large portion of the literature on parameter estimation in computational mechanics, and we can identify combinations and settings that have not yet been addressed. Moreover, by including stochastic methods, the quality of the identified parameters can be addressed as well. Among the various available methods for parameter inference, we compare Bayesian and frequentist approaches and derive new results for two-step identification methods, which are needed to calibrate complex models.

The remainder of this paper is structured as follows. In Sec. 2, we start with an overview of the basic equations of the initial boundary value problems, as well as the structure of the material models. We then discuss a few aspects of the experimental possibilities, followed by a brief review of parameter identification in elasticity, viscoelasticity, and viscoplasticity.

Section 3 treats the numerical approaches to identify the material parameters. We formulate the LS method based on the FEM, the equilibrium gap, and the VFM, as well as more recent methods such as PINNs and the combination of model selection and parameter identification (i.e., model discovery).

As stated earlier, one of the major aims of this paper is to formulate a unified treatment for most of the available parameter identification procedures, which is the subject of Sec. 4. Section 5 focuses on the quality of the identified material parameters—which is addressed in a statistical setting, including the topics of identifiability and uncertainty quantification. In Sec. 6, examples are provided to compare the performance of the various schemes for selected applications.

This section addresses the basic problem of parameter identification. First, the fundamental equations of solid mechanics are summarized. Then, we provide an overview of the basics of experimental observations and constitutive modeling using elastic, isotropic, or anisotropic, and inelastic material models. The notation in use is defined in the following manner: geometric vectors are symbolized by $a$ and second-order tensors $A$ by bold-faced letters. Furthermore, we introduce matrices and column vectors by bold-faced italic letters $A,a$⁠.

In this paper, we thus do not treat parameter estimation in dynamical systems, where properties of wave propagation in solid media are used for identification purposes, see Ref. [4], for example.

where $Div$ defines the divergence operator with respect to $X$⁠. Moreover, $ρR=(detF)ρ$ is the density in the reference configuration, and $P=(detF)TF−T$ symbolizes the first Piola–Kirchhoff stress tensor, where $F=Grad χR(X,t)$ represents the deformation gradient. Here, $Grad$ denotes the gradient operator with respect to $X$⁠.

The balance of angular momentum implies the symmetry of the Cauchy stress tensor, $T=TT$⁠.² Thus, six independent stress components have to be determined, whereas only three scalar PDEs are available. Hence, the system of Eqs. (2), or equivalently Eq. (3), is closed by formulating

1. the kinematics, i.e., a relation between the displacement vector $u(X,t)=χR(X,t)−X$ (or the motion $χR(X,t)$⁠) and the deformation gradient $F$⁠, and

2. a constitutive model describing the stress state in dependence of the deformation and possibly of additional, so-called internal variables describing the process history.

For isotropic materials, $Ŝ(C)$ is an isotropic tensor function [6–8]. In the case of anisotropy, special invariants define the stress state, see Ref. [7]. For a so-called hyperelastic constitutive model, a strain energy density function $Ψ(C)$ defines the stress state as $S=2ρRdΨ(C)/dC$⁠. To limit the scope of the presentation, only the concise notation in Eq. (4) will be used here. Of course, other representations are possible, e.g., it is also common to express an elastic material law as $T=T̂(B)$ with the left Cauchy–Green tensor $B=FFT$⁠, however this choice does not essentially influence the following considerations of material parameter identification.

if the last equation defines the yield condition, see Ref. [11]. In this case, Eq. (5) represents a DAE-system, see Refs. [11] and [12]. Moreover, in the field of yield-function-based, rate-independent modeling, the right-hand side $rq(C,q)$ contains case distinctions (loading-unloading or Karush–Kuhn–Tucker conditions) to be able to reproduce a nonsmooth material behavior. It should also be noted here that we are using variables defined in the reference configuration. When using internal variables with respect to the current configuration or models based on the multiplicative decomposition of the deformation gradient, if the evolution equations of the internal variables depend on relative derivatives with respect to some intermediate configuration, special attention must be paid to time integration in order to obtain objective quantities. To avoid a cumbersome discussion, we assume here that the relative derivatives can be expressed by the material time derivative of a tensorial quantity relative to the reference configuration with a transformation usually referred to as a pull-back operation. This implies a consistent notation in view of the interpretation of considering DAE-systems later on, for details see Refs. [13] and [14].

Different approaches in constitutive modeling derive systems of the form (5) from two scalar potentials, namely, the strain energy density and the dissipation rate potential, so that (under certain conditions on the properties of such potentials) thermodynamic consistency is automatically fulfilled, see Refs. [15–18], for example. Alternative approaches evaluate only the strain energy density function within the Clausius–Duhem inequality and motivate the evolution Eq. (5)₂ in a different manner [7].

Of course, there are also other constitutive models that do not fit into the aforementioned structure—such as models described by integro-differential equations, fractional derivatives [19,20], or the endochronic plasticity formulation in Ref. [21]. Alternatively to the original formulation in Ref. [21], the stresses can also be given by a rate equation using a particular kernel function of the endochronic formulation, so that the integral formulation can be transformed into a rate formulation [22]. Models described by variational formulations, see, for example, Ref. [23], do not fit directly into the explained structure as well. However, very similar equations are obtained after their discretization, so that the schemes discussed later on are directly transferable.

$t¯$ is the given traction vector and $u¯$ the prescribed displacement, whereas $u¯0(x)$ and $q¯0(x)$ are the known initial conditions. Mixed boundary conditions, where $t=t̂(u)$⁠, are allowed as well. Consistent initial conditions imply that Eq. (7)₁ is also fulfilled for the initial conditions. Formulation (7) covers a large number of models.

In the following, we discuss special cases of material models and issues when determining material parameters from experimental data. We start with experimental possibilities and outline measurement techniques in Sec. 2.2, while parameter identification for different classes of constitutive models is discussed in Sec. 2.3.

The constitutive models in Eqs. (7)$2,3$ contain unknown parameters $κ$⁠. They have to be determined in such a manner that they reflect experimental observations. Since the models contain stresses and strains, the question of how to deduce these quantities from experimental measurements is a crucial one. The principal difficulty is to find experiments where the stress state under given (i.e., controlled) or measured external loads is known. Thus, before addressing material parameter identification, we briefly discuss the possible options for mechanical testing and the related detectable experimental data.

First, we consider experiments where the stress tensor $T$ is uniform in a certain area of the specimen, and we start with the uniaxial tension/compression test as the simplest example. In this case, only Eq. (7)$2,3$ have to be solved, since Eq. (7)₁ is automatically fulfilled (provided that $ρb$ is negligible and the material properties are homogeneously distributed). In the data analysis, it is necessary to consider the specific structure of the tensorial quantities due to the assumed boundary conditions, such as the absence of stress components in the lateral direction for the present example of a uniaxial test. It should also be noted that, at any time, testing machines such as the ones depicted in Fig. 1 can either control the displacement and measure the resulting force—or vice versa. Accordingly, in order to obtain information on strains and stresses, it is necessary to know the relationship between the displacements and the strains as well as the forces and the stresses. For a uniaxial test, both relationships are known, but care must be taken to account for rigid body displacements (which occur in a specimen due to the finite machine stiffness) as well as to perform a sufficiently accurate measurement of the cross-sectional area of the specimen. A further experimental choice could be to directly control the strain as a function of time rather than the displacement. This choice would require additional technical equipment and call for a synchronization of the strain measurement device with the testing machine. Of course, a temporal quasi-static process can be divided into two (or more) stages with different controls, e.g., with displacement control during loading and force control during unloading.

Torsion tests on thin-walled tubes are an attractive alternative to tensile tests, since a uniform stress state over the wall thickness can be assumed if the wall thickness is small in comparison to the tube radius. During testing, either the torsion angle or the torsional moment is controlled. Once again, in order to compute strains and stresses, we need the relations between the torsion angle applied by the testing machine and the shear strain, as well as between the torsional moment and the shear stress. For the latter, the uniform stress assumption can be exploited.

Three-point or four-point bending tests are another simple test alternative. However, these experiments generate nonuniform stress fields, so that Eq. (7)₁ needs to be solved. One option is to directly solve Eq. (7)₁ as a local equilibrium equation (viewing the specimen as a three-dimensional solid), and another one is to recast it as a global equilibrium equation (adopting the approximations of beam theory). Respectively, the material parameters are calibrated based on full-field measurements or on discrete deflection information, see Refs. [24] and [25], for example.

For large deformation cases, there exist deformations—called universal deformations—which fulfill the local balance of linear momentum under the assumptions of isotropy and incompressibility, see, e.g., Refs. [8] and [9]. In all other cases, the entire problem (7) has to be solved.

Apart from the experiments we just discussed, there are only very few other alternatives, such as a specific type of shear test, see Refs. [26] and [27] for an overview, or experiments on specific membranes, see Ref. [28], from which the stress state can be extracted, for example. Even a biaxial tensile test, see Fig. 1, does not provide a uniform stress state, see Ref. [29] for a discussion on uniformity and a measure of deviation from it. In general cases, only the forces or torques of the testing device are directly known, and assumptions on the stress state must be provided.

On the other hand, if optical access to the sample surface and the ability to identify surface patterns are given, displacements (and strains) can be determined locally on the surface using, e.g., digital image correlation (DIC, see Refs. [30] and [31]), as shown in Fig. 2. With this technique, out-of-plane strains cannot be obtained without introducing further assumptions. In some situations, strain determination is only possible in an integral sense and in one direction, employing strain gauges, clip-on strain transducers, or optical methods. Fiber–Bragg grating sensors, see Refs. [32–34], are also used to determine strains in one direction inside a component, although care must be taken to ensure that these measurement systems do not influence the local strain state. It is also possible to apply microcomputed tomography (μ-CT) in combination with digital volume correlation (DVC) for volumetric determination of the strains, but the recording times are currently so long that long-term changes in the material, such as creep or relaxation effects, can only be ruled out by estimation [35,36]. Additional information can be obtained by the “single-valued” displacement or angle of the specimen holder.

In some situations, there is no optical access to the specimen, e.g., in triaxial tension/compression testing machines [37], indentation tests, Refs. [38–40], or some metal-forming processes where the specimens are not visible in the forming press. In these cases, only a single displacement component and a single force component, both as functions of time, can be evaluated by interpolating discrete data. In some cases, of course, a combination of full-field information and single-valued data can be extracted and considered in the material parameter identification concept [41].

The parameter identification procedure depends on the particular constitutive model class at hand. Thus, the choice of the constitutive model class, which is based on the observed mechanical response of a specimen, e.g., the presence of plastic deformations or viscous effects, directly influences the approach to identifying the sought parameters. In the following, we provide an overview of material parameter identification for a few important classes of constitutive models—elasticity, viscoelasticity, elastoplasticity, and viscoplasticity.

For linear elasticity at small strains, we distinguish between isotropic and anisotropic cases.

so that the identification of tensorial relations reduces to that of scalar-valued functions. The axial stress is given by $σ:=T11=e1·Te1=EE11=Eε$ (with $ε:=E11=e1·Ee1$⁠). The lateral strains $εq:=E22=E33=e2·Ee2$ are given by $εq=−νε$⁠, i.e., the Poisson's ratio ν can be determined if the lateral strain $εq$ is measured. Thus, we obtain a decoupling of the identification of Young's modulus E and Poisson's ratio ν, whereby E can be determined from the stress–strain information in the axial direction and ν can be obtained by measuring the transverse strain. In both cases, furthermore, the material parameters can be determined from linear expressions.

Then, the identification of the parameters has to be carried out, for example, by a nonlinear LS approach, see Sec. 3.3.1. Moreover, the lateral strain data are of importance since the axial information alone is insufficient to determine K and G uniquely. A detailed discussion of the reduction of the general class of material models (7)$2,3$ from the three-dimensional (3D) case to the uniaxial tensile case, which involves incorporating the boundary conditions into the strain and stress state, is provided in Ref. [42].

Whether material parameters, such as K and G in this case, can be uniquely identified should be linked to a criterion. A first attempt to investigate whether material parameters can uniquely be identified is provided in Ref. [43]. There are situations where infinite combinations of parameters can yield a very accurate reproduction of the experiments, but this raises the question of the physical meaning of the parameters. This is discussed in the general case of optimization by Ref. [44] and later on by Ref. [45], referred to as local identifiability. This concept is transferred to the case of problems in solid mechanics by Ref. [43] and further discussion is provided by Ref. [46]. A similar approach using the terminology of local stability is discussed by Refs. [47–49]. In a more general context, a first attempt to consider quality measures of the parameter identification in solid mechanics can be found in Ref. [50], where the correlation between the parameters is investigated. Here, the local identifiability concept, which can be extended to general models, is discussed in Sec. 5.1.1.

An extension of parameter identification for linear elastic materials to the case of anisotropy, and particularly of transverse isotropy and orthotropy (requiring five and nine material parameters, respectively, for a 3D model), is discussed in Ref. [51], unfortunately only at the theoretical level and with no supporting experiments. In fact, it is not possible to devise an experiment that covers one assumed boundary condition that is connected to a volumetric deformation. If applying a hydrostatic pressure using a fluid, it is very difficult to measure the strains in the three different directions. Inserting the specimen into a chamber whose walls are force gauges, the sensitivity to friction may become a dominating factor. In a 3D-testing machine, a homogeneous stress/strain state becomes questionable. Thus, not all parameters can be uniquely inferred by the procedure in Ref. [51]. This issue is addressed in Ref. [52] for transverse isotropy. It is shown by analytical considerations that there is one volumetric deformation process required—apart from tensile and shear deformations—to uniquely obtain all parameters. To implement this deformation process, a compression tool is developed [52]. However, the evaluation of the data shows that the measured lateral stresses are very sensitive to the geometric accuracy of the specimen and of the tool itself, and also to the friction conditions within the testing equipment. Consequently, the measurement results are not reliable and show large fluctuations. For the case of transverse isotropy with a specific preferred direction (perpendicular to the isotropy plane) in axisymmetric problems, the identification is discussed in Ref. [53].

The case of orthotropy is addressed in Ref. [54], where μ-CT data are evaluated to obtain the material parameters appearing in the analytical results obtained by assuming a homogenized material response. The difficulty is again to perform a reliable compression test to identify one of the parameters. For plane problems, the number of parameters can be reduced, see Ref. [55], for example. Again, the main difficulty is the need for special compression tests and their evaluation, since these tests provide only a very small amount of measurement data, which is also uncertain. Even DIC, which delivers only the surface deformation, does not circumvent the problem of making all material parameters “uniquely” identifiable. The problem is partly solved by using representative volume elements with individual material properties of the constituents to numerically estimate the homogenized material behavior.

The identification of the material parameters in hyperelasticity (elasticity at large strains, whereby the existence of a strain energy density function is postulated) has to be discussed in more detail. First, hyperelasticity relations are chosen either as part of the models of viscoelasticity, rate-independent plasticity, and viscoplasticity, or they are chosen to model purely elastic material behavior. The latter is typically the case for rubber and some biological materials, which are commonly modeled as weakly compressible. Thus, we discuss the case of isotropy and incompressibility first.

Incompressibility is modeled either by an undetermined pressure, which is calculated by the geometrical constraint equation of no volume change [6,7], or by introducing a bulk modulus that is infinitely large (in practice, very large) in comparison to the shear modulus. In the latter case, the bulk modulus can be interpreted as a penalty factor that is defined by a large number—and not determined by parameter identification. This choice of the penalty factor, i.e., the bulk modulus, is made by the user and can lead to lateral deformations that do not precisely represent incompressibility. Further, the resulting lateral stretches might even become nonphysical at very large strains [56]. For a possible modification requiring appropriate forms of the strain energy density function for the compressible part, see Ref. [57]. With the choice of incompressibility, the experimental data of the lateral stretch from a uniaxial tension/compression test are not necessary.

In the most common hyperelastic material models—apart from the micromechanically motivated approaches, e.g., Ref. [58]—the strain energy density is expressed by a polynomial function of the first and second invariants of the right Cauchy–Green tensor [57,59], or by a polynomial in the eigenvalues of the right stretch tensor (eigenstretches) [60]. The latter are called Ogden-type models. The main problem in identifying the parameters from uniaxial tensile tests, torsion tests, or combined tension-torsion tests is that some parameter values might lead to highly oscillatory tension-compression results or to highly sensitive results, which is discussed in Ref. [61] for a number of submodels of Rivlin–Saunders-type. These phenomena occur especially outside the range of the measurement data used for calibration.

One possible approach to solve this issue is to assume a priori that all polynomial coefficients are positive [62]. This leads to a linear LS identification problem with inequality constraints for uniaxial tension and torsion (universal deformations). The inequality constraints (enforcing positivity of the parameters) can be connected to stability requirements [63], and they raise the general question of the existence of a solution in relation to the constraints on the material parameters. This is discussed in Ref. [57] for a new class of polynomial models in the framework of polyconvex strain energy functions. Here, the assumption of non-negative polynomial coefficients (material parameters) is even necessary.

Further models permit arbitrary signs for the parameters, see Refs. [64] or [65]. So far, it has not been proven that negative material parameters exclude a physical behavior in all deformation stages. However, special attention is suggested when applying such models with negative parameters.

As mentioned earlier, in Ogden-type models, the strain energy density function depends nonlinearly on the eigenstretches. At first consideration, this has advantages for simple tests where the loading directions coincide with the eigendirections of the stretch tensor. Here, however, the material parameters occur nonlinearly in the stress–strain relation, necessitating the solution of a nonlinear optimization problem even for simple cases [66–68]. It is hardly possible to guarantee the sensitivity of the result in material parameter identification (due to the choice of the initial guess of the parameters in an iterative approach) or, consequently, the uniqueness of the estimated parameters.

For a recent overview of identification in hyperelasticity, see Ref. [69], where appropriate weighting functions are chosen to identify the parameters. A comparison of hyperelasticity models can be found in Refs. [70–72]. An overview of analytical expressions for determining the material parameters for specific isotropic hyperelasticity models is given in Ref. [73].

The treatment of material parameter identification for large strains in connection with anisotropic materials is similar as with small strains. The applications here are primarily designed for biological tissues, as these have preferred directions due to the presence of collagen fibers. In Ref. [74], the anisotropic material behavior of a soft tissue is experimentally investigated and calibrated using a finite strain, anisotropic hyperelasticity relation. Even the quality of the parameters is studied, which is very difficult to be guaranteed in layered tissues as it is the case of arteries, see Ref. [75]. Reference [76] considers experiments on various tissues and calibrates the material parameters based on data given in the literature. Unfortunately, the quality of the parameters is not studied. In Ref. [77], a polyconvex model for transverse isotropy is calibrated based on virtual experiments generated with a different constitutive model. In this respect, we refer to Ref. [78] as well. Alternative applications of anisotropic hyperelasticity are woven fabrics, see Ref. [79]. There, however, the simpler case of a plane stress problem is considered which essentially reduces the number of unknown parameters.

In the case of viscoelasticity, we have to treat either integral equations, ODEs, or DAEs to consider fading memory properties in the material. For this type of model, after infinitely long holding times of the applied load, the stress state coincides with the equilibrium stress state (which has no hysteresis). One common model structure goes back to overstress-type models, i.e., the stress state is decomposed into an equilibrium stress state, formulated by an elasticity relation, and an overstress part. This model type is sometimes also called hyper-viscoelasticity. The overstress part can be formulated as a sum of Maxwell models, i.e., ODEs of first order. Each Maxwell element contains two material parameters. Since all the parameters are highly correlated, a suitable procedure must be developed to determine the parameters in a reproducible way. Another approach takes a fractional derivative model as a surrogate model to determine the individual parameters for given relaxation spectra [80]. Unfortunately, the consistent reduction from the three-dimensional modeling to the uni-axial tension is not addressed. An alternative scheme is proposed in Ref. [81], where inequality constraints are used to successively determine the parameters. For a recent investigation and literature survey, see Ref. [82].

Since the model has a modular structure, the parameters related to the equilibrium stress and to the overstress can be determined successively, i.e., one arrives at the identification of the parameters of an elasticity relation (discussed in the previous subsections) and of the parameters of an ODE system. In Ref. [9], the successive identification is addressed under the assumption of a universal deformation (tension–torsion). However, no quality measures for the identification are discussed. Even the specification of the identification procedure is missing. Because of problems in identification, singular value decompositions are employed in Ref. [83].

Two questions arise here: First, how can the boundary conditions for homogeneous deformations be included in three-dimensional material models (7)$2,3$? Second, which methods are available to determine the parameters in ODEs and DAEs in the context of LS approaches? The first issue is treated in Ref. [42], where a procedure for arbitrary models is proposed. The aspect of determining the parameters in transient problems is summarized in Ref. [84], where three types of schemes are proposed, namely, simultaneous simulation of sensitivities [85–87], internal numerical differentiation [88], and external numerical differentiation—see Ref. [89] for an application in nonlinear FEM. For the case of viscoelasticity using full-field measurements, which is discussed in Sec. 3.1, the reader is referred to Refs. [90] and [91].

There is a large amount of articles on parameter identification in the context of plasticity. However, the quality of the estimated parameters is not frequently addressed. Among the first papers to investigate the correlation between parameters are Refs. [92] and [93]. Models based on yield functions contain elastic parameters, whose identification is discussed in Sec. 2.3.1. The determination of the yield stress is much more difficult for most metals. Either there are Lüders bands resulting in a spatial motion of dislocations, so that the resulting stress–strain response is not deterministically predictable, or the yield point is not very pronounced. Thus, the first yield stress is a very uncertain value. In the plastic region, a number of material parameters can describe the nonlinear hardening behavior (isotropic and kinematic hardening, for example), and these parameters can be strongly correlated. A linear correlation between the driving and the saturation term of the Armstrong & Frederick kinematic hardening model, see Ref. [94], can be concluded by analytical considerations of this model [42,54]. For pressure-dependent yield functions in soil mechanics, the identification using a triaxial compression testing device is discussed in Ref. [95]. The quality of the estimated parameters is also examined. The contribution [96] discusses the identification of material parameters for a complex viscoplasticity model using a gradient-based method. A detailed explanation of the calculation of the sensitivities is provided, which fits into the internal numerical differentiation procedure [88].

The more complex the models are, the more difficult it is to determine the parameters from simple experiments. Therefore, material parameter identification can serve as a means of model improvement and possibly model reduction (model identification). However, this field is often not documented, as only the final result of a model is published. When devising new material models, attention must be paid to their identifiability from experiments already during the development phase. This is discussed in detail in Ref. [54]. The authors of Ref. [97] discuss the practical identifiability of the parameters for a more complex yield function, whereas Ref. [98] focuses on a nonlinear kinematic hardening model that is calibrated from uniaxial tensile tests on the basis of one-dimensional considerations by applying an evolutionary algorithm. A further discussion on parameter identification using stress data from tensile experiments is provided by Ref. [99]. This is extended by investigations on noisy data by Ref. [100]. A review regarding the calibration of various plasticity models is also provided by Ref. [101].

As mentioned earlier, only very few experiments lead to uniform stress and strain states. If nonuniform stress/strain distributions occur within a specimen, problem (7)—consisting of (a) the balance equations (PDEs), (b) the constitutive equations, and (c) the kinematic relations, in addition to the initial and boundary conditions—has to be evaluated to obtain the parameters. Since the resulting systems of equations are similar for uniform deformations and spatially discretized nonuniform deformation problems, we will now discuss the latter, more challenging case. Also, the emphasis here is on the identification procedure and not on the numerical optimization schemes.

This section is structured as follows. First, details on the space and time discretization using finite elements are given, together with a detailed description of the three representative problems of linear elasticity for small strains, hyperelasticity, and inelasticity (these problems are later discussed from the perspective of identification). Then, a broad classification of the different identification methods is provided in Sec. 3.2, while an in-depth comparative study is deferred to Sec. 4. We start our review of identification methods with the nonlinear LS method using the FEM (Sec. 3.3) and continue with the equilibrium gap method and the VFM (Sec. 3.4). Surrogate models and PINNs are the topic of Sec. 3.5. Finally, an overview of model discovery and Bayesian approaches is provided in Secs. 3.6 and 3.7, respectively.

where the matrices $Zq e,j ∈Rnq×nQ$ represent data management matrices containing only zeros and ones (Boolean matrices). Also, $nGPe$ is the number of integration points (commonly Gauss points) within element e, and $nQ=(∑e=1nelnGPe)×nq$ defines the total number of internal variables of a mesh with $nel$ elements.

which represent the nodal internal forces, whereas $p¯(t) ∈Rnu$ defines the given equivalent nodal force vector. Here, $Z e ∈Rnue×nu$ and $Z¯ e ∈Rnue×np$ symbolize incidence matrices assembling all element contributions into a large system of equations (representing the assembly procedure). Moreover, $Ee,j ∈R6$ denotes the column vector representation of the element strain tensor (Voigt notation), which depends linearly on the element nodal displacement vector $ue ∈Rnue, Ee,j=Be,jue$⁠, with $ue=Z eu+Z¯ eû$⁠. The number of element nodal displacement degrees-of-freedom is denoted with $nue$⁠, and $we,j$ denote the weighting factors of the spatial integration in an element. Furthermore, $Be,j ∈R6×nue$ is the strain-displacement matrix of element e evaluated at the jth Gauss point, $j=1,…,nGPe$⁠. $Je,j ∈R3×3$ symbolizes the Jacobian matrix of the coordinate transformation between reference element coordinates and global coordinates. The symmetric stress tensor (7)₂ is recast into vector $T=T̂(Ee,j,qe,j) ∈R6$⁠, which is evaluated at the jth Gauss point. It depends, via the strain vector, on the displacements $ua$ and the internal variables q.

is commonly solved, whereas Eq. (28)₂ is used to compute the nodal reaction forces during postprocessing. Although it is common to state that a Newton–Raphson method is chosen to solve problem (29), it was shown in Refs. [11] and [104] that the Multilevel Newton algorithm proposed in Ref. [105] is applied (if an iterative stress algorithm is chosen at the local level, i.e., at each Gauss point). This method leads to the classical structure of local iterations, often called stress algorithm (although it is the internal variable computation), and to global iterations, where the increments of the displacement vector are obtained (on the basis of the consistent linearization stemming from the usage of the implicit function theorem).

Clearly, for rate-independent material models such as elasticity or rate-independent plasticity, time is only used to parameterize the external load and has no physical meaning (for this reason, it is sometimes referred to as “pseudo-time”). However, in experimental investigations time is also supplied, as the responses (such as displacements and forces) measured at different times are incorporated in the parameter identification. Therefore, the time dependence is also included here. Frequently, the stepwise increase of the load is chosen to be close to the solution of the Newton–Raphson method applied to Eq. (30)₁, see Ref. [91] for details. In this sense, time is considered at least in the time-dependent boundary conditions.

Here, the explicit time (or load) dependence is omitted since proportionality is given. The stiffness matrices are defined by

and the elasticity matrix $Ce,j ∈R6×6$⁠.

In the following, we formally restructure the equations discretized above with the aim of material parameter identification. For this purpose, we redefine $y ∈Rns$ as the state vector and introduce the vector of material parameters $κ ∈Rnκ$⁠. The structure of the state vector depends on the problem under study.

In the forward problem, $κ$ is given, i.e., $κ=κ¯$ and $F(y,κ¯)=0$ has to hold. However, in parameter identification, $κ$ has to be determined. Usually, no time dependence is assumed within the framework of statics. To indicate that different load values are used for parameter identification (and going back to the time-continuous setting for notational simplicity), Eq. (38) should actually be $F(t,y(t),κ)=0$⁠. While $Aup(κ)$ in Eq. (40) is time-independent, since the stiffness matrices do not depend on time, $f(t,κ)$ is time-dependent since it contains the time-dependent prescribed displacements $u¯(t)$ and equivalent nodal force vector $p¯(t)$⁠.

In the time-discrete setting, Eq. (46) has to be solved for each load or time-step, yielding different parameters $κ$⁠. With this approach, further considerations have to be made to finally determine a unique set of parameters (e.g., by averaging, which is the simplest possibility).

There are hyperelastic material models that are linear in the material parameters, such as the class of Rivlin–Saunders or Hartmann–Neff models, see Refs. [57] and [59]. This leads to the same representation as in Eq. (42). Unfortunately, this is not the case for Ogden-type models [8,60], where the fully general case (47) has to be considered.

is given, see Eq. (15). If an implicit time discretization is applied, $y˙$ disappears from the equation. However, y contains the unknowns at all discrete times, and F consists of all nonlinear systems to be solved at each of those time steps. Now, as introduced in Eq. (17), y contains the unknown nodal displacements u, the nodal reaction forces p, and the internal variables q evaluated at the Gauss points. In this case, it is also possible to compute the reaction forces, an exercise which is left to the reader.

where Eq. (24) defines F in the inelastic case. Thus, once $ŷ(t,κ)$ is known, the parameters alone have to be determined, which is denoted as the reduced approach, see Sec. 4.3.

Calibration is used to determine model parameters with which the model can best approximate available measurement data. A criterion for the approximation accuracy is provided in terms of a closed-form mathematical expression, i.e., the objective or loss function. The latter is a norm of the difference between the available measurement data and the model predictions, obtained from the solution of the governing equations via a numerical (discretization) method or a surrogate model. One important characteristic of the different numerical methods presented here is the type of spatial (and/or temporal) discretization. After formulating an objective function, the optimal model parameters that minimize this function are identified using optimization algorithms, often referred to as optimizers. Alternatively, in a statistical setting, a single best guess of the unknown material parameters is formulated as a point estimator, i.e., a function of the observed data [106]. Here, the likelihood function plays a central role, and Bayesian approaches additionally consider a prior density. With statistical models, the uncertainty of an estimate can be assessed through confidence/credible intervals, which are typically computed with sampling approaches.

In Table 1, the methods presented in the remainder of this paper are classified by the aforementioned features, namely,

1. Discretization: Galerkin or collocation methods, as well as local or global ansatz functions.

2. Parametrization: The objective or likelihood function depends on parameters $κ$ of the material model, and possibly also on the PDE solution or on a parametrization of the latter.

3. Optimization/sampling: While deterministic calibration approaches mostly make use of different gradient-based optimization schemes, both optimization and sampling are common choices in a statistical setting.

If the material parameters cannot be identified on the basis of simple experiments, where the stress and the strain are known, the entire boundary-value problem (7) of the experiment must be solved. This can be done by numerical approximation methods such as finite differences, boundary elements, finite volumes, or the FEM. In the following, we assume that the latter is the method of choice. The nonlinear LS (NLS) method is applied to minimize the difference between the results of the finite element simulation of the experimental boundary conditions and both pointwise and full-field measurement data. First, we discuss proposals for parameter identification of phenomenological constitutive models (macroscale), followed by schemes addressing models with scale separation, specifically materials with a microstructure.

The combination of the FEM and a LS approach to determine parameters $κ$ in the finite element model goes back to Ref. [107]. Later, Ref. [108] linked LS and FEM with the goal to detect the location and the Young's modulus of an inclusion in a matrix material. Discrete data in combination with finite elements were also used in Refs. [40] and [109–113]. This approach can be followed for experiments with nonuniform stress states when access is limited to, e.g., resultant forces, such as in the case of indentation tests or metal forming processes, see Refs. [114] and [115].

A conceptual and practical step beyond the use of pointwise data has emerged with the availability of optical methods. Here, a special focus lies on DIC methods, see, e.g., Refs. [30] and [116], where the surface displacements of the specimen are measured during loading. Approaches taking advantage of full-field displacement data, pioneered by Mahnken [47,117,118], were further extended by Refs. [29], [50], and [119–126], see also Ref. [127] for biaxial tensile tests and Refs. [128] and [129] for more complex deformation cases. In Ref. [1], an approach of this type was denoted finite element model updating (FEMU), a terminology that has been used continuously since then. Unfortunately, this terminology can be misleading since the same name is well-established in structural dynamics, see Ref. [130] for a review. Following Refs. [130] and [131], the FEMU approach is applied when the mathematical structure of the problem, here a finite element equation (e.g., the entries in the stiffness or mass matrix), is changed during the system identification process, see Ref. [132], for example. In contrast, when calibrating constitutive models, the finite element equation is specified once—and only the material parameters are updated. As a result, the calibration of constitutive models has to be clearly classified as parameter identification rather than FEMU. In Ref. [91], the combination of NLS with DIC data, where the boundary-value problem is discretized using the FEM, is denoted as NLS-FEM-DIC—a denomination that includes the objective function, the boundary-value problem solution technique, and the experimental measurement technique.

Alternative approaches to account for optical information are considered in Ref. [133], where gratings on the specimen surface are incorporated, in Ref. [134] using Moiré-patterns, or in Ref. [135] and [136] using contour data and discrete points on the surface. Reference [1] provides an overview of other schemes, such as various “gap methods” mainly applied in model updating using vibrational data, and a comparison to determine the elastic parameters, see Ref. [1]. Further contributions propose the “constitutive relation error” and “modified constitutive relation error” approaches, Refs. [137] and [138]. A study of the integrated DIC method and its references is provided by Ref. [139]. Overall, the testing and identification paradigm—which takes advantage of full-field measurement data rather than simple strain transducers or strain gauges—in conjunction with the FEM is referred to as “Material Testing 2.0” in Ref. [31].

The previously mentioned approaches bear a strong relation to the mathematical literature on the least-squares method applied to the solution of ODE- or DAE-systems [84]. This is explained in detail in Ref. [89] and will be briefly summarized in the following.

Since the number of individual entries, their magnitude, and their physical units vary within the residual r, it is common and reasonable to introduce a diagonal weighting matrix $W ∈RnD×nD$ and to consider the weighted residual $r˜(κ)=W r(κ)=W{s(κ)−d}$⁠. For an approach to account for different amounts of data and different magnitudes of the physical quantities within the data, see Ref. [61]. The steps to preprocess the experimental full-field data as well as the data obtained from the numerical simulation are explained in detail in Appendix  A.

represents the sensitivity matrix (also denoted as functional matrix or Jacobian matrix). The solution of problem (56) can be computed by various methods, see Refs. [140–144]. In the case of gradient-based algorithms, the derivatives (57) are required to compute $κ*$⁠. Parameter identification with sensitivity assessment is also covered by Ref. [145]. The computation of the sensitivities is compiled in Appendix  D.

A systematization of the NLS approach with DIC data is presented in Ref. [146], where various large deformation problems are treated. However, the problem of local identifiability is not addressed, see Sec. 5.1.1 for details.

In some cases, materials possess a heterogeneous microstructure (an example being fiber-reinforced composites), and the goal of the identification is to determine their macroscopic properties, i.e., the parameters describing a homogenized material. For unidirectional fabrics or for orthogonal fabrics within the linear elastic regime, for example, appropriate macroscopic models are transversely isotropic elasticity (with five parameters) and orthotropic elasticity (with nine parameters), respectively. The difficulties that arise here lie in the experiments required to determine these parameters uniquely (see Sec. 2.3.1). One alternative is to use a representative volume element (RVE) and, assuming that the material parameters of the constituents are known, to determine the homogenized material parameters with analytical mixing rules, as discussed in Refs. [147] and [148], or with computational homogenization. However, care must be taken to ensure that the appropriate deformation modes are excited in order for the related parameters to be determined; see Ref. [54] for a specific application of an RVE to determine orthotropic material parameters. Further questions are treated in Refs. [83] and [149], discussing the issue of identifying the inelastic properties of the matrix and the interphase parameters for composite materials.

Another task might be to identify the material parameters of the microscopic constituents of a heterogeneous material. This can be done either by considering only the RVE and applying to it appropriate homogeneous loading conditions, or by relying on macroscopic experiments. The latter option is discussed by Ref. [150] in a FE² context on the basis of the LS approach. Unfortunately, local identifiability of the parameters is not discussed. Since FE² computations are very expensive, an iterative identification procedure using standard FE² codes is very time-consuming. This problem is pointed out in Ref. [151], however, no connection to an identifiability assessment is made. The use of integrated DIC for parameter identification within a multiscale approach is discussed in the recent contribution [152].

An example of RVE-based identification approach for a damaging material is discussed in Ref. [153]. In order to incorporate the microvoid behavior under loading into the model, a unit cell is considered and the material parameters controlling damage evolution are identified. For that application, it is important to emphasize that considering damage at RVE-scale must be handled with care during homogenization. In addition to ensuring that appropriate deformation modes are excited for the parameter identification, it is important to recognize that the assumption of separation of scales is no longer applicable. Consequently, higher-order homogenization procedures should be used.

As discussed in Sec. 3.3, the NLS method with DIC data seeks to minimize the square of the weighted residuals $r˜(κ)$ for the unknown parameters $κ$⁠, thus minimizing the mismatch between FEM predictions $s(κ)$ and experimental observations d.

which implies that, in general, $F(d,κ)=0$ is not fulfilled (“there is a gap to the equilibrium conditions”). Moreover, d must be known. It has to be remarked that $u=du$ (where $du$ are the experimental displacement data) is very restrictive, as it only allows the use of the region of the DIC measurements (which is a subset of the domain) and limits analysis to plane problems if no three-dimensional data are available. While the scheme is applicable for problem classes I and II, it requires additional considerations for problem class III since the internal variables are not measurable. Thus, the scheme is only applicable for very specific problems in the proposed form (59). Theoretically, the scheme is extendable to force data as well. However, this may be difficult when using DIC, as the forces corresponding to the imaged regions may be unknown in this case.

Often, the number of virtual fields in $V$ is chosen such that it coincides with the number of unknown material parameters. These functions must be linearly independent.

which is a linear system with $nκ$ unknowns. The same property also holds for problem class II (hyperelasticity) if the material parameters are linearly embedded, e.g., with Rivlin–Saunders and Hartmann–Neff models [57,59], see Ref. [159]. Thus, material parameter identification with the VFM requires only the solution of one linear system of equations, whereas the NLS–FEM approach requires running a forward FEM simulation at each iteration of the optimization algorithm. Note, however, that the full displacement field is here assumed to be available from experimental measurements. A discussion on exceptions is provided in Ref. [160].

For problem class III (models with internal variables), the scheme explained above is not applicable. Details regarding the application of the VFM to some specific models with internal variables are provided by Refs. [157] and [161]. For general considerations, we refer to Sec. 4.4.1. Please note that, given the conceptual similarities between equilibrium gap and VFM, we will only refer to the VFM in the following.

The evaluation of $s(κ)$ using, e.g., a FEM simulation may be costly, especially if many iterations are required and the model to be calibrated is computationally demanding. Therefore, it can be efficient to replace the simulation with a so-called surrogate or metamodel. The vector of the corresponding results is denoted as $ssurr(κ)$⁠.

Surrogate models require a flexible parametric regression ansatz, such as, e.g., multivariate polynomials. Recently, machine learning approaches such as artificial neural networks (ANNs) and Gaussian processes, often referred to as Kriging, have gained increasing attention. A brief introduction to ANNs and the notation used in the paper is given in Appendix  C. One of the first attempts to apply ANNs as surrogates for the identification of material parameters goes back to Refs. [38,39], and [162]. Recent work in this direction can be found in Refs. [163] and [164]. When using surrogate models in the context of parameter identification, the overall procedure is the following:

1. To train a surrogate model $ssurr(κ;θ)$ as a regression function through the data sampled from a given number of evaluations of the physical model $s(κ)$⁠. In the context of Gaussian processes, $θ$ are the kernel parameters that can be fitted to the data using Bayesian model selection or by minimizing the negative log-likelihood. After fixing the kernel parameters, the GP surrogate is conditioned on the available data to obtain the regression function. The training of ANNs, on the other hand, requires the identification of trainable parameters $θ$⁠, i.e., the weights and biases, via optimization.

2. To identify the physical parameters by solving the minimization problem (63), whereby $s(κ)$ is replaced with its surrogate $ssurr(κ;θ=θ*)$⁠. Herein, $θ*$ denotes the nonphysical parameter set of the surrogate identified in step 1. Note that surrogates can also be used in the context of statistical inference in combination with sampling approaches, see Secs. 3.7 and 5.

Step 1 and 2 can be interconnected, which is then referred to as adaptive sampling, see Ref. [165] for an extensive review in the context of Kriging.

the ANN acts as a function approximator (or ansatz function) of the PDE solution. The above ansatz $U$ is parameterized in the ANN weights and biases $θ$⁠. Figure 3(a) shows a schematic representation of ansatz (64).

Spatial (and, in the case of dynamic problems, temporal) derivatives of the ansatz (64) can be calculated using automatic differentiation, see Refs. [168,170], and [171]. For dynamic problems, this approach is referred to as continuous-time PINN [172]. Alternatively, temporal derivatives can be treated by means of discrete time-integration schemes, see Refs [172] and [174]. As mentioned earlier, dynamic problems are out of the scope of the representation. Instead, the time dependence in Eq. (64) expresses that different load (time) steps within one experiment or various experiments are used to identify the set of material parameters $κ$⁠.

Figure 3(b) shows a schematic representation of a parametric PINN ansatz in comparison to the one in Eq. (64). The training of parametric PINNs, i.e., the identification of the ANN parameters $θ$⁠, is discussed in detail in Sec. 4.3.2. Parametric PINNs have been deployed for thermal analysis [175], magnetostatics [176], or for the optimization of an airfoil geometry [177]. To the best of our knowledge, they have not yet been used for the calibration of material models in solid mechanics.

An alternative to the parametric PINN approach discussed above are so-called inverse PINNs, which are related to the all-at-once formulation introduced in Sec. 4. In inverse PINNs, material, and network parameters $κ$ and $θ$ are identified simultaneously, and the optimization problem (63) does not hold anymore. Details are provided in Sec. 4.5.2.

The following contributions deal with material model calibration from full-field displacement and force data using PINNs. In Ref. [178], the authors propose a multinetwork model for the identification of material parameters from displacement and stress data for linear elasticity and von Mises plasticity. For the analysis of internal structures and defects, the authors in Ref. [179] present a general framework for identifying unknown geometry and material parameters. In Ref. [180], a framework is developed for the calibration of hyperelastic material models from full-field displacement and global force–displacement data. Unlike in the conventional PINN approach [172], the physical constraints are imposed by using the weak form of the PDE. In Ref. [181], variational forms of the residual for the identification of homogeneous and heterogeneous material properties are formulated. The method is demonstrated using the example of linear elasticity. The study in Ref. [182] considers the identification of heterogeneous, incompressible, hyperelastic materials from full-field displacement data. It uses two independent ANNs, one for the approximation of the displacement field, and another one for approximating the spatially dependent material parameters. A comparison of different data sampling strategies as well as soft and hard boundary constraints is reported in Ref. [183]. In Ref. [184], PINNs are further enhanced toward the calibration of linear elastic materials from full-field displacement and global force data in a realistic regime. The realistic regime here refers to the order of magnitude of material parameters and displacements as well as noise levels.

The previously discussed methods for material characterization focus on the calibration of material parameters in an a priori known material model. The a priori choice of the material model, i.e., the combination of mathematical operations that describes the material behavior, typically relies on the intuition and modeling experience of the user. Inappropriate assumptions in the material model are reflected in a poor fitting accuracy of the model even after parameter calibration. To avoid such modeling errors, strategies have been proposed recently to automatically discover a suitable symbolic representation of the material model, while simultaneously calibrating the corresponding material parameters. Thus, the problem of material parameter calibration is generalized to the more intricate problem of material model discovery.

Algorithms for discovering symbolic models from data can be broadly classified into two types: symbolic regression methods based on genetic programming [185], which repeatedly mutate a model until its agreement with the data is satisfactory, and sparse regression methods [186] which select a model from a potentially large set (also called library or catalogue) of candidate models based on data. Such methods have first been used in the physical sciences to deduce symbolic expressions of the governing equations of dynamical systems [187,188]. In the field of material modeling, symbolic regression based on genetic programming has been used since the early work of Ref. [189], see also the more recent works in Refs. [190–194]. Sparse regression for the discovery of material models has been much less explored and has only recently gained attention [159,195–197]. See also the review in Ref. [198], where a novel symbolic regression approach based on formal grammars is proposed for hyperelasticity.

In Ref. [159], the authors propose EUCLID (efficient unsupervised constitutive law identification and discovery), a method for discovering symbolic expressions for hyperelastic strain energy density functions using sparse regression starting from displacement and force data. The method was later extended to elastoplasticity [199], viscoelasticity [200], and generalized standard materials [201], see Ref. [202] for an overview. A supervised version (i.e., a version based on stress data) of the EUCLID approach is experimentally validated in Ref. [203] using data stemming from simple mechanical tests on human brain tissue.

The idea behind EUCLID is to construct a set of candidate material models (i.e., a material model library) by introducing a general parametric ansatz for the material behavior, which depends on a large number of parameters (⁠$κ∈Rnκ$ with $nκ≫1$⁠) comprising all material parameters of all candidate material models. Calibrating the parameters using one of the previously discussed calibration methods would result in a highly ill-posed inverse problem. Even assuming that such a problem is solvable, the solution would deliver a parameter vector $κ$ with many nonzero entries and, thus, a highly complex symbolic expression for the material model. Therefore, a regularization term is introduced which penalizes the number of nonzero entries in $κ$⁠, and, consequently, promotes simplicity of the symbolic material model expression while alleviating the ill-posedness of the inverse problem. The so-called $ℓp$-regularization term takes the form $||κ||pp=∑i|κi|p$ with $0<p≤1$⁠, assuming high values for dense vectors $κ$ and smaller values otherwise. In the limit $p→0+$⁠, the $ℓp$-regularization term converges to an operator that counts the number of nonzero entries in $κ$⁠. For p = 1, the approach is called least absolute shrinkage and selection operator (Lasso)-regularization.

where the first term in the objective function quantifies the mismatch between the model and the data. Here, similar to the VFM, see Eq. (60), the model-data mismatch is defined indirectly by the sum of squared residuals of the discretized weak form of the balance of linear momentum. Here, however, we should note that other measures for the model-data mismatch, see Sec. 3.3, could be chosen likewise. The matrix $A(d,u¯)$ and the load vector $pˇ$ in Eq. (67) are constructed as described in Appendix  B for linear elasticity at small strains. If more sophisticated material behavior like hyperelasticity is considered, the matrix $A(d,u¯)$ changes accordingly, see Ref. [159].

The effect of the sparsity-promoting regularization term on the minimization problem is influenced by the weighting factor $λ>0$⁠, which can be chosen to strike a balance between the fitting accuracy and the complexity of the discovered material model, see Refs. [199–203]. Abandoning the deterministic setting, the EUCLID method is investigated from a Bayesian perspective in Ref. [160] to quantify the uncertainty in the discovered models.

In the literature, two different lines of research on model discovery can be distinguished. The first aims to discover symbolic and thus interpretable expressions for the material models, as in the approaches we just described. The second is concerned with the training of noninterpretable black-box machine learning models for encoding the material response. Examples of the latter include machine learning models like splines [204], Gaussian processes [205], neural ordinary differential equations [206], and—the arguably most popular choice—ANNs [207]. Although purely data-driven constitutive models have shown some success [207,208], much effort has been addressed to constrain ANNs such that fundamental physical requirements (such as objectivity, material stability, and thermodynamical consistence) are respected, while maintaining their expressivity. In the context of hyperelasticity, for example, Refs. [209–212] choose specifically designed ANN architectures to ensure the (poly-)convexity of the strain energy density function. In the context of dissipative materials, Refs. [197,213–215] consider ANNs that do not violate thermodynamic requirements (see Ref. [216] for a review). An often disregarded drawback of ANNs is their data-hungriness. Supervised training frameworks for ANNs require a large amount of labeled stress–strain data pairs, which are experimentally not accessible and are therefore typically acquired through computationally expensive microstructure simulations, which require the microstructural properties of the material to be known. To counteract the data scarcity, unsupervised training frameworks that solely rely on experimentally available displacement and reaction force data like the NN-EUCLID proposed by Ref. [211] or the approach presented in Ref. [217] are indispensable. For further details on data-driven material modeling and discovery, the reader is referred to the recent review on the topic in Ref. [218]. While the subsequent analysis focuses on EUCLID, note that in principle any data-driven constitutive model can be implemented within a FEM or VFM code for model discovery, see, for instance, Ref. [217].

where $p(κ)$ refers to the parameter prior and $p(d|κ)$ represents the likelihood function, which can be evaluated by solving the mechanical model.

with the weighted norm $‖x||A2=xTAx$ for any positive definite matrix A, and recover the NLS approach with Tikhonov regularization, see Ref. [220] for details and further links between Bayesian approaches and regularization. One can even obtain information about the posterior based on this optimization-based approach, by perturbing the data d with a randomly drawn noise vector according to the density $πe$⁠. This so-called randomize-then-optimize approach is investigated in a number of papers, see Ref. [221], for instance.