est destinée au dépôt

Nouvelles stratégies numériques pour l'identification robuste, cohérente et efficace de modèles à partir de mesures de champs

Hai Nam Nguyen

To cite this version:

Hai Nam Nguyen. Nouvelles stratégies numériques pour l'identification robuste, cohérente et efficace de modèles à partir de mesures de champs. Solid mechanics [physics.class-ph]. Université Paris-Saclay, 2021. English. NNT : 2021UPAST061. tel-03260879

Albert Einstein

The way to research is not a mystery, yet it is a great deal of work. It comes with lots of ups and downs, lost and found moments. Fortunately, I am blessed to have mentors, family and friends that went through the journey with me.

I am thankful for my guide in this endeavour, Prof. Ludovic CHAMOIN, who showed me the direction and encouraged me to run. I appreciate the fruitful talks with Dr. Cuong HA-MINH who provided me with support not only in academic but also in life aspects. I am filled with gratitude to my supervisors for their continuous support, precious guidance, valuable inputs and patience in all phases of this dissertation. It was truly my honor to receive such a great support.

Most importantly, I am indebted to my mother and father whose love and support brought me so far and encouraged me to go the extra mile.

Finally, I would like to express my very special, heartfelt gratitude to my wife and my son who are always with me on every journey I do. Without their sacrifices, remarkable understanding and endless support with full of love in four years I was away from home, I could not have made this far.

i 1.8 DIC procedure with a global approach (from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]), the measured displacement field being obtained from a FE mesh (that is subsequently used for comparison with FE numerical simulations). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Fitted Poisson ratio versus mesh size (in pixels) in a biaxially loaded elastic specimen: coarse meshes lead to poor kinematic representations whereas fine meshes lead to high noise sensitivity (from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF] 1.12 Tensile test on dog-bone sample with a hole (left), and associated FE mesh with measured and extruded Dirichlet boundary conditions (right), from [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. . . . . . . . . . . . . . . . . 1.13 Global DIC results for the analyzed test: displacement fields expressed in pixels (top and center left) and normalized gray level residual with respect to the dynamic range of the picture in the reference configuration (bottom left); corresponding in-plane strain fields (right), from [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. . . 1.14 Sensitivity fields for some of the final set of material parameters, from [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. The plotted fields are [S i ]p i in order to share the same units and make comparisons. . . . . . . . . . . . . 1.15 I-DIC residuals for the analyzed test: displacement residuals, normalized gray level residuals, and load residuals for the first (left) and converged (right) solutions, from [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. . . . . . .

List of Figures

3.1

Example of adaptive mesh I-DIC procedure, from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]. The element size in the adapted mesh (right) is a decreasing function of sensitivity k@u s /@pk...................89

Evolution of

A ε for the considered 1D non-periodic case (left), and MsFEM basis functions (right). . . 6.4 Comparison between the exact solution and the approximate MsFEM solution: primal fields u ε H and u ε (left), gradients du ε H /dx and du ε /dx (center), and fluxes A ε du ε H /dx and A ε du ε /dx (right). . . . . . . 6.5 Illustration of local problems to be solved on the macro element K with prescribed tractions fK .....147 6.6 Artificial speckle considered (left), and evolution of E(x) in the domain (right). . . . . . . . . . . . . . .

Introduction & motivations

"Essentially, all models are wrong, but some are still useful" (G.E. Box)

In structural mechanics and computational engineering, complex physical systems are now commonly studied by means of modeling and simulation tools, in order to analyze their mechanical behavior or make predictions. In this context, systems are thus represented by physics-based models, described by PDEs, that provide an idealized mathematical abstraction of the physical reality. Out of recent trends that aim at representing and predicting from full data-based information, such mathematical models are the basic ingredient for further numerical analysis which is a daily industrial and research practice; they constitute a fundamental pillar in science and engineering activities. For a given system, the mathematical representation may be picked in a hierarchical list of possible models, with increasing complexity. The selection is traditionally performed from a priori knowledge on the system; merely speaking about the material behavior representation, one may for instance select a simple linear homogeneous elasticity law, but nonlinear multiscale laws may be preferred to represent complex physical phenomena at various scales. The chosen mathematical model is then further numerically processed by means of a discretization scheme such as the finite element method (FEM), leading to a numerical model used as a virtual twin.

Insights of computational approaches depend on the models being a faithful abstraction of the real world, but all models are wrong to some extent. In the framework of computational mechanics based on FE analyses, there are various uncertainties and error sources along the modeling and simulation chain, coming from: (i) partial and noisy experimental information (observations) on the system feeding the virtual representation; (ii) model bias due to an always imperfect representation of the physical system behavior and its environment (related to ignorance, variability in the process, or intentional reduction of complexity); (iii) numerical approximation in terms of discretization error associated with computational solution schemes. For the sake of quantitative computed information, accurate prediction, and safe decision making, there is thus a need for certification of the computed outputs, with consistency between physical, mathematical, and numerical models. This is the matter of verification and validation (V&V) approaches, with procedures for model identification (or updating) from available experimental information, model verification controlling the accuracy of numerical methods (e.g. mesh adaptivity), and model validation (Figure 1); all this should be complemented with uncertainty quantification in order to get quantitative confidence information.

V&V has been an active research topic for more than 30 years [Roache, 1998, Oberkampf et al., 2003]; it has been listed as one of the most challenging applications of simulation-based engineering sciences (SBES) [START_REF] Oden | Simulation Based Engineering Science -Revolutionizing Engineering Science through Simulation[END_REF].

Nevertheless, there is also an increasing need for computational efficiency, with fast simulation and predictions. This is a typical constraint in industry to accommodate engineering times and affordable simulations, and it is becoming a major requirement for online control of systems from modeling and assimilation of in-situ measurements, in the framework of Dynamic Data-Driven Application Systems [Darema, 2004, Darema, 2015]. In this framework, which falls into the broader concept of Cyber-Science & Technology involving a new generation of connected and autonomous systems [START_REF] Russel | Artificial Intelligence: a modern approach[END_REF], computational and physical capabilities are integrated by means of embedded sensors, processors, and controllers in order to perform a numerical feedback loop in which experiments and simulations synergistically and continuously interact (see Figure 1). In modern computational engineering, the goal is thus to compute right at the right cost, with appropriate physics and a smart management of computing resources (trade-off between reliability and computational cost) depending on the objective and on information coming from the physical system. This resorts to model adaptivity in terms of discretization, model upscaling, but also appropriate selection of a reference mathematical model (multi-fidelity approach) with respect to available experimental data. It furthermore requires the use of advanced numerical methods (such as model reduction techniques) in order to address computationally demanding problems with both fast and credible numerical strategies. In this PhD work, we consider all these issues and requirements in the context of model identification or updating from physical observations. This is a fundamental topic as nowadays, and after two last decades that have seen major advances in terms of computational infrastructure and methods to address complex problems, a recurrent limitation in the reliable numerical simulation of material or structural mechanics problems comes from the setting of input model parameter values. We specifically focus on the case of full-field measurements that have also been a major breakthrough of the last two decades, with much increased availability of data. In particular, quantitative imaging by means of Digital Image Correlation (DIC) (or Digital Volume Correlation (DVC)) techniques is an advanced experimental technology that has become a very active research topic over the last two decades (see [Hild and Roux, 2012b] for a review). It is now a common practice for the acquisition of measured displacement fields on materials and structures in research labs and industry; a typical illustration, among others, is the research development of the Correli software [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF], Roux and Hild, 2008, Leclerc et al., 2015, Neggers et al., 2018].

Combined with FE numerical simulations and inverse analysis, DIC is a tool of choice in order to accurately identify the material behavior from full-field measurements. These permit to retrieve several parameters from a single heterogeneous test. Measurements coming from DIC are experimental information which is large but with limited accuracy (as any measurement), due to noise coming from both image acquisition and post-processing. For robustness and effectivity in the V&V process, this uncertainty on observations (which can be quantitatively characterized a priori) needs to be taken into account in a consistent manner. It not only needs to be propagated throughout the identification procedure, but it also needs to be compared and considered with respect to other error sources (modeling, discretization) which are involved in the identification process and which may question identification relevance if not contained or compensated [Mottershead andFriswell, 1993, Simoen et al., 2015].

In recent years, there have been some important works to make the DIC process more efficient in terms of credibility of the obtained experimental information from images (e.g. [START_REF] Pan | Fast, robust and accurate digital image correlation calculation without redundant computations[END_REF],Tomicevic et al., 2013]). Much less progress has been made when it comes to the identification problem itself, i.e. to the way the associated computational model is defined and managed (in terms of ignorance or limited accuracy) from the content of experimental information coming from DIC. On the one hand, an appropriate mathematical (physics-based) model with optimized complexity should be considered as a representation of the observed physics, in order to perform a consistent interpretation and comparison with available data. If the model is too coarse, it will not be able to account for all the richness that is included in experimental observations; if it is too complex, it will lead to expensive computations with no beneficial (and potentially detrimental) impact on identification results, as measurement noise prevents from capturing this extra-level of complexity. It is also intended that the model ignorance be learned from informative DIC data, when the situation occurs. On the other hand, and due to the specific features of the experimental information which is considered here (full-field measurements, not sparse sensing), a special attention should be paid to the mesh size selection (associated with discretization error) in order to exploit measurement information at best with meaningful accuracy; a typical example is given in [START_REF] Hild | Toward 4d mechanical correlation[END_REF]. All these aspects are the topic of the present research work.

Most of applications using DIC measurements do not consider modeling of discretization error; the first source is usually poorly known a priori, while the second source is usually neglected by assuming sufficiently fine meshes (which is questionable). We wish here to propose an alternative vision in which all a priori knowledge and uncertainty sources (coming from both model and experimental data) are considered, informed, connected, and propagated throughout the identification procedure in order to build a consistent and unified methodology. In particular, we wish to inform on bias in the numerical model (coming from an uncertain constitutive law, or inexact fulfillment of equilibrium) in order to better account for the measured data. In this way, the proposed identification procedure appears as a balance between errors in the numerical model and noise corrupting the measurements, with modeling and simulation accuracy governed by that of available experimental information, and providing for the best efficiency for robust, consistent, and effective model updating purposes.

We consider the previously mentioned aspects, referring to predictive data science and smart data assimilation, by designing a suitable approach to combine DIC experimental information and numerical tools. We address this issue by introducing an appropriate identification metric that merges a measurement metric and a metric in the modeling space. This is performed by referring to the deterministic Constitutive Relation Error (CRE) concept associated with duality aspects. This concept, strongly based on physics and material science notions, has been used for a long time in V&V activities (a review can be found in [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF]). It has sound foundations in material and computational mechanics, and provides for quantitative information on the quality of a model. Its modified version, designed for identification/updating and denoted mCRE [Ladevèze et al., 1994a, Deraemaeker et al., 2002], showed performance in many applications and appears as a convenient tool for the PhD objectives. Indeed, the mCRE approach is based on splitting the knowledge on the problem into a reliable part (which is strongly enforced in the inversion process) and an uncertain part (which is relaxed and satisfied at best). The associated hybrid variational formulation thus resorts to reliability of information and provides for a natural metric mixing measurement and modeling uncertainties. Out of various properties (convexity, robustness) and higher performance compared to some alternative inversion techniques, as reported in several works [Ben Azzouna et al., 2015, Waeytens et al., 2016], the mCRE approach directly defines a data-based model enrichment that takes into account features of the physical solution which are not represented a priori in the considered physical model. This is performed by recovering the best state estimation (in terms of so-called admissible fields) from both mathematical model and experimental data, similarly to alternative recent methodologies putting forward the concept of hybrid twins [Maday et al., 2015a, Peherstorfer and Willcox, 2015, Chinesta et al., 2020].

In the context of full-field measurements, the mCRE has been surprisingly employed in quite few works [START_REF] Calloch | Identification de modèles de comportement de matériaux solides : utilisation d'essais et de calculs[END_REF], Ben Azzouna et al., 2015, Barbarella et al., 2016, Hild et al., 2016], despite the fact that DIC (for the processing of images) and mCRE (for model inversion) share the common feature to make benefit of the a priori mechanical knowledge. It is shown in this PhD work that the mCRE concept is a dedicated tool to perform correct (and at correct cost) modeling and simulation for identification purposes, for both linear and nonlinear models and from full-fields measurements. In particular, it naturally provides error indicators on modeling and discretization that identify the computational model inaccuracies and permit to select the optimal one depending on the quality of experimental data. The associated adaptive strategy is in terms of constitutive model class (among a manifold of model classes with various complexity levels), but also in terms of discretization mesh. We mention that the obtained adaptive mesh, dedicated to parameter identification and implicitly considering sensitivity analysis and experimental uncertainties, is naturally different from the one that would be designed for accurate prediction in a direct problem (see [START_REF] Arridge | Approximation errors and model reduction with an application in optical diffusion tomography[END_REF] for a similar adaptive mesh philosophy in another inverse analysis context). In the mCRE context, it is shown that the a posteriori construction of a fully equilibrated admissible stress field, verifying the balance equations in a strong sense, enables one to obtain estimates on both discretization and modeling errors, which can then be compared with measurement noise in order to drive model and mesh adaptivity. Therefore, bias in the numerical model (deviations between model predictions and reality observed through full-field measurements) can be corrected and adjusted so that model complexity and outputs remain consistent with physical observations and noise level. Furthermore, the adaptive procedure can also be seen as a way to filter useless or redundant experimental information for model identification purposes. The proposed approach is first introduced and validated on linear elasticity models, before being extended to more complex nonlinear (elastoplasticity) material behaviors.

In addition, the PhD addresses the use of advanced numerical approaches in combination with mCRE. These aim at complementing the robust inversion method with powerful numerical tools in order to enhance the computational efficiency [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. One aspect which is investigated is the use of reduced order modeling (ROM) which facilitates the solution of the multi-query parametric problems which are involved in the mCRE inversion process (as in any other inversion process). The employed ROM technique is the Proper Generalized Decomposition (PGD) [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulation[END_REF], and the focus is in its use for the minimization of the mCRE functional and for the construction of equilibrated fields, which usually are some computationally intensive procedures. Another aspect is the setting of a specific and efficient solver, based on the LATIN method [Ladevèze, 1989], when considering mCRE with nonlinear material behaviors. A last aspect, which appears as a preliminary study, is the design of a multiscale approach based on MsFEM [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF] when addressing image-based measurements and mCRE with complex material microstructures.

The manuscript is organized as follows:

• a bibliographic review is performed in Chapter 1. It focuses on general aspects of inverse problems, on the (modified) CRE concept applied to address such problems, and on model identification based on DIC;

• the general implementation of the mCRE concept in association with full-field measurements is addressed in Chapter 2, in which a consistent framework taking measurement noise into account is designed. An integrated version of the coupling between image-based experimental techniques and mCRE is also proposed;

• in Chapter 3, a specific interpretation of the modeling error term involved in mCRE is made, in order to inform on the quality of the numerical model (mesh, constitutive relation) with respect to experimental information.

An adaptive strategy is proposed in this context for optimal model and mesh selection;

• advanced numerical tools referring to reduced order modeling and based on the Proper Generalized Decomposition (PGD) technique are developed in Chapter 4. They enable the computational cost to be decreased and the numerical performance of the overall identification process to be enhanced;

• an extension of the proposed inversion procedure to nonlinear models is performed in Chapter 5, and is complemented with a specific solver. Model adaptivity is again illustrated in this context, showing how an initial coarse model can be automatically enriched in order to be consistent with experimental information and uncertainties;

• another extension of the inversion strategy to multiscale problems using the MsFEM framework is proposed in Chapter 6;

• last, conclusions are listed and prospects to the work are drawn.

Chapter 1

Bibliography

In this first chapter, we give a bibliographic review on two topics which are at the heart of the PhD work. The first one is the (modified) CRE concept, which provides for a specific variational formulation with energy-based functional for model identification or updating. The second one deals with the acquisition of full-field measurements from DIC, and their post-processing for identification purposes. We show here the general philosophy and highlight the main features on these two topics.

Modified CRE concept in the context of inverse analysis

Basics on inverse problems

Context

Material and structural mechanics, as other scientific fields, involves models with unknown or at least not perfectly known parameters. A typical case, which is considered here, is constitutive parameters i.e. parameters involved in the modeling of the local material behavior. For credible numerical simulation, analysis, and prediction, a classical procedure is then to infer the missing knowledge on model input parameters by exploiting experimental information on the response of the mechanical system, in terms of indirect measurements from sensors (accelerometers, strain gauges, fiber-optic sensors, cameras. . . ). This refers to model identification or updating, which falls into the larger framework of inverse problems in which the best fit is searched by comparing model outputs and available observations. Inverse analysis is an important topic of mechanical engineering with many applications such as health monitoring and defect identification [START_REF] Mottershead | Model updating in structural dynamics: a survey[END_REF], Doebling et al., 1996, Andrieux and Bui, 2006], geophysical exploration [Menke, 2012], or biomechanical imaging [Arridge, 1999, Barbone andGokhale, 2004], among others. As stated before, it is also of paramount importance for the identification of unknown material parameters (e.g. elastic moduli), possibly heterogeneous, of constitutive laws that feed computational mechanics models.

Mathematically speaking, we are given a parametrized model M(p) chosen in a manifold of possible models, representing a physical system of interest, and generically described by a PDE of the form A M (u; p)=0 . The forward problem consists in computing the state u and related output quantities d(u) by means of the model and from specified input parameters p 2 P (usually chosen dimensionless). Conversely, the inverse problem consists in learning parameters p from overabundant (and usually noisy) observations d obs (Figure 1.1). The idea is to find, by comparison between d(u(p)) and d obs , the parameter set that represents at best experimental data. This ill-posed problem (in the Hadamard sense, with issues on uniqueness and continuity of the solution with observed data) can be solved by means of two families of approaches:

• deterministic approaches in which the inverse problem is formulated as a constrained optimization problem with minimization of the discrepancy between model outputs and observations [Bui, 1993, Cailletaud and Pilvin, 1994, Kern, 2002, Bonnet and Constantinescu, 2005]. Approaches differ one from another in the nature of the imposed constraints and on the choice of the cost function that defines (in a certain metric) the distance between computed and measured data [START_REF] Roux | Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields[END_REF]. Such deterministic approaches will be mainly considered in the PhD work;

• stochastic approaches (referring to statistical inverse analysis) in which the inverse problem is formulated in a probabilistic manner, using the Bayesian framework [START_REF] Kaipio | Statistical and Computational Inverse Problems[END_REF], Tarantola, 2005, Stuart, 2010, Vigliotti et al., 2018]. Several variants can be derived depending on additional assumptions which are made. A major concern, which will be further discussed, is the fact that the inversion process is plagued with several error and uncertainty sources. On the one hand, there is measurement error due to noise and possibly systematic error (caused by imperfections in the measurement setup or during subsequent signal processing). Considering an additive measurement noise ✏ (that may be characterized or not), the observed outputs read d obs = d true + ✏ with d true the physical (undisturbed) output value. On the other hand, the model contains bias as it never exactly describes the behavior of the system (due to errors when translating physics into equations, and when solving these equations numerically); there is thus no value of p such that d(u(p)) = d true . Measurement and model uncertainties combine in the total error d(u(p)) d obs obtained as the difference between model prediction and observed quantities. Consequently, the results of an identification/updating procedure are by definition uncertain, and one should account for this uncertainty by returning not only the identified value of the parameter set but also confidence information (e.g. in terms of covariance matrix).

Remark . The balance between the richness of experimental information and the number of parameters to be identified is an important aspect of effective inversion. A practical issue in inverse problems, which will not be addressed here, is the case of a high-dimensional parametric space as encountered for instance with discretized parameter fields (spatially distributed parameters). In such a case, and with sparse measurements, a strategy to manage numerical difficulties (e.g. ill-conditioning) is to a priori reduce the number of parameters by using projection-based methods [START_REF] Teughels | Damage assessment by fe model updating using damage functions[END_REF], Lieberman et al., 2010] or sensitivity analysis methods which aim at keeping parameters having the largest influence on outputs of the model [START_REF] Saltelli | Sensitivty analysis as an ingredient of modeling[END_REF], Saltelli et al., 2004, Borgonovo and Plischke, 2016, Neggers et al., 2017]. Also, the richness of experimental information may be increased by an optimal placement of sensors (experimental design) which can also been performed from sensitivity analysis approaches [START_REF] Kleiber | Parameter sensitivity in nonlinear mechanics: theory and finite element computations[END_REF], Daescu and Carmichael, 2003, Waeytens et al., 2013], or using information entropy criteria [Papadimitriou, 2004, Huan andMarzouk, 2013].

Deterministic inversion approaches

Deterministic approaches for solving inverse problems are based on the definition of a suitable cost (or misfit) function F(p) that reports on the correlation between the numerical model outputs and the experimental observations. This function is minimized over the parameter space P in order to get the solution p sol of the inverse problem:

p sol = argmin p2P F(p) (1.1)
Typical and most intuitive cost functions are (weighted) least squares functions also referred as squared residuals, based on a measure in the (weighted) Euclidian norm of the mismatch between computed and measured responses [START_REF] Cailletaud | Identification and inverse problems related to material behaviour[END_REF], Mahnken and Stein, 1996, Doyley et al., 2006, Brigham and Aquino, 2009].

They are of the form:

F(p)= 1 2 (d(u(p)) d obs ) T G 1 obs (d(u(p)) d obs ) (1.2)
with G obs a weighting matrix that is assumed diagonal in many practical applications, with coefficient (G obs ) ii proportional to 1/d 2 obs,i for the sake of normalization. Alternative cost functions may be defined, such as these based on auxiliary fields and reciprocity gap [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF], Calderon, 2006], or energy-based functions as defined later with the CRE concept.

A well-known difficulty of deterministic approaches, related to the ill-posed and/or ill-conditioned features of inverse problems, is that the cost function F is usually not convex due to a lack of experimental data (underdetermined minimization problem), with several local minima and high sensitivity to small perturbations in these data. To circumvent this issue and increase robustness, the inversion process is usually complemented with a regularization procedure that consists in replacing the initial inverse problem with an approximate inverse problem such that the committed error is compensated by a better conditioning [Bui, 1993[START_REF] Engl | Regularization of inverse problems[END_REF]. A typical class of regularization methods is Tikhonov regularization [Tikhonov andArsenin, 1977, Tikhonov et al., 1995], in which a priori knowledge or constraints (smoothness, sparsity. . . ) on parameter values is informed by augmenting the cost function with a regularizing quadratic term bringing a convex component [START_REF] Engl | Regularization of inverse problems[END_REF], Hansen, 1998, Neumaier, 1998, Titurus and Friswell, 2008]. As an example, when a prior guess p 0 is known, the regularized cost function may read:

F r (p)=F(p)+ ↵ 2 (p p 0 ) T G 1 0 (p p 0 ) (1.3)
The balance of the two terms that constitute F r is given by the magnitude of the dimensionless regularization parameter ↵ that should be set conveniently: if too small, the regularization effect may be insufficient and the problem may remain ill-posed; if too large, the solved problem has little connexion with the original one. Various methods may be employed to set an optimized value of ↵, such as (i) the L-curve method [Miller, 1970,Hansen, 1992,Hansen, 1998, Vogel, 2002] that heuristically identifies the highest curvature point in a log-log plot of the regularization term versus F (balance of the two terms) with varying ↵; (ii) the Morozov discrepancy principle [Morozov, 1968, Morozov, 1984, Nair et al., 2003] that adapts the regularization parameter to the noise level in the data (assuming this level is known a priori). This latter method will be further discussed in Chapter 3.

As regards the practical solution of the inverse problem (1.1) (or its regularized version), it is usually performed by means of nonlinear numerical optimization schemes with iterative algorithms fed by local Jacobian and/or Hessian information [START_REF] Gill | Algorithms for the solution of the nonlinear least squares problem[END_REF]. As a basic method, the steepest descent (gradient) scheme

p (k+1) = p (k) r (k) p F = p (k) S (k) T p G 1 obs (d(u(p (k) )) d obs )
can be used even though it is associated with slow convergence. S p = @d(u(p)) @p is the Jacobian (or sensitivity) matrix, which may be computed from the adjoint state method or a perturbation method. The Newton-Raphson scheme

p (k+1) = p (k) H (k) 1 r (k) p F = p (k) H (k) 1 S (k) T p G 1 obs (d(u(p (k) )) d obs )
with H = @ 2 F @p 2 the Hessian matrix (assumed to be full column rank), can also be used; it has much faster convergence but also higher computational cost. In most situations, quasi-Newton methods (such as BFGS) or inexact Newton methods are preferred due to their ease of implementation (only the computation of gradients is required) [START_REF] Nocedal | [END_REF]Wright, 2000, Oberai et al., 2003]. We may list in this last category the Gauss-Newton method

p (k+1) = p (k) (S (k) T p G 1 obs S (k) p ) 1 S (k) T p G 1 obs (d(u(p (k) )) d obs )
or the Levenberg-Marquardt method

p (k+1) = p (k) (S (k) T p G 1 obs S (k) p + µI d ) 1 S (k) T p G 1 obs (d(u(p (k) )) d obs )
where µ is a damping parameter. As an alternative to previous local methods, which are sensitive to the initial guess and converge to local minima in the absence of convexity, global methods (e.g. the simulated annealing method) may be used to explore the whole space. Eventually, in extreme cases with non-differentiable functions, the simplex method [START_REF] Nelder | A simplex method for function minimization[END_REF] or genetic algorithms may also be implemented.

Stochastic approaches

Stochastic approaches have been used in several engineering applications of inverse problems, particularly for the monitoring and control of structures [Beck, 2010]. They have also been recently investigated in the context of identification from full-field measurements [START_REF] Zhang | Robust multiscale identification of apparent elastic properties at mesoscale for random heterogeneous materials with multiscale field measurements[END_REF]. They are based on a probabilistic description of uncertain model inputs as random variables with corresponding probability density. Consequently, the result of the inverse problem is a conditional probability density function (pdf) on parameters p, that represents the degree of belief on the value of these parameters given some observations, and that can further be used to derive confidence intervals or for uncertainty quantification on model outputs of interest [START_REF] Rebba | Validation and error estimation of computational models[END_REF]. Stochastic approaches lead to an automatic regularization procedure as the whole space of possibilities for parameter values is explored.

The classical mathematical framework for inverse analysis in a stochastic context is Bayesian inference [START_REF] Beck | Updating models and their uncertainties. i: Bayesian statistical framework[END_REF], Tarantola, 2005, Kaipio and Somersalo, 2004, Stuart, 2010, Rosic et al., 2013]. It proposes a complete probabilistic description which naturally and rigorously considers measurement and modeling (epistemic)

uncertainties, as well as natural irreducible (aleatory) uncertainties. The formalism refers to the Bayes theorem that reads (for a given model class M):

⇡(p|d obs )= 1 C ⇡(d obs |p).⇡(p) (1.4)
where

• ⇡(p) is a prior pdf on p, constructed from a priori knowledge on parameter values. It is chosen independently of measurement results, but rather based on engineering judgment. Several methods have been developed to obtain consistent prior pdfs, such as the conjugate prior method (where the prior pdf is chosen in such a way that the posterior pdf belongs to the same family of distributions), or methods based on information theory e.g.

the principle of maximum entropy [Gull, 1988, Soize, 2008]. Nevertheless, in situations where a large amount of informative observations is at hand (as in DIC), the choice of the prior pdf is not extremely determinative since almost any prior (sufficiently smooth in the region of high likelihood) leads in practice to similar final conclusions;

• ⇡(p|d obs ) is the updated or posterior pdf on p, obtained after taking the observations d obs into account. This conditional probability is the end result that quantifies the learning permitted by measured data, reducing uncertainty and giving the most likely values of p. A solution to the inverse problem may thus be derived from the mean value or the maximal value (maximum a posteriori (MAP) estimator) of ⇡(p|d obs ). Out of specific cases, the description of ⇡(p|d obs ) is usually analytically intractable (e.g. due to complex physical models) and remains implicit. Furthermore, its exploration in potentially high-dimensional spaces needs to be performed [START_REF] Marzouk | Stochastic spectral methods for efficient bayesian solution of inverse problems[END_REF], Rubio et al., 2018];

• C is a normalization constant, which ensures that the posterior pdf integrates to 1. It is usually not computed, but is equal to the probability of having d obs for the given model class M. This is often referred to as the evidence for M that plays a determining role in Bayesian model class selection [START_REF] Beck | Model selection using response measurements: Bayesian probabilistic approach[END_REF].

With additional assumptions (e.g. Gaussian pdfs, uncorrelated variables), a linearized and lighter version of Bayesian inference can be obtained. This is the framework of Kalman filtering [Kalman, 1960] and its extensions to nonlinear models (extended Kalman filtering (EKF), unscented Kalman filtering (UKF) [START_REF] Julier | A new extension of the kalman filter to nonlinear systems[END_REF])

that are well-known in signal processing or robotics, and that have been used in many applications for the solution of inverse problems [START_REF] Hoshiya | Structural identification by extended kalman filter[END_REF], Corigliano and Mariani, 2001, Mariani and Corigliano, 2004, Yang et al., 2006,Mariani and Ghisi, 2007,Liu et al., 2009,Moireau and Chapelle, 2011,Marchand et al., 2016]. Kalman filtering considers low-order statistics (mean, covariance) of parameters alone, so that the computational cost is much lower than in the case where a full probabilistic model is considered.

Alternative strategies for modeling and propagating uncertainties in inverse problems have been developed [Zim-mermann, 1991, Helton et al., 2004]. In particular, in the case of small uncertainties, some techniques avoid a full probabilistic study or make it more affordable. We mention here (i) perturbation methods based on first-or secondorder expansion of random variables; (ii) non-probabilistic methods based on interval analysis (propagation of value range intervals) such as fuzzy set approaches [START_REF] Haag | Model validation and selection based on inverse fuzzy arithmetic[END_REF], Simoen et al., 2015] relying on the fuzzy set theory [Zadeh, 1965], or hybrid interval approaches e.g. the random set theory or the Lack-Of-Knowledge theory [START_REF] Ladevèze | Sur une théorie des méconnaissances en calcul des structures[END_REF], Ladevèze et al., 2006, Louf et al., 2010] in which interval bounds are considered to be random variables. These last interval approaches may be interesting when only very little information concerning the uncertain input parameters is available, as carrying out a complete probabilistic study becomes unreasonable in this specific case.

Eventually, and contrary to parametric approaches [Schueller, 2007] which are considered here, non-parametric approaches [Soize, 2005] based on random matrix theory may be used to construct a probabilistic model and implement it in inverse analysis.

Remark . Besides model parameters, measurement and modeling errors may also be modeled as random variables and identified/updated from observed data in the Bayesian inference process. This avoids incorrect or unsuitable assumptions that unfairly influence the Bayesian updating results. Usually, variance parameters of uncorrelated zero-mean Gaussian representations are included [START_REF] Zhang | A comprehensive bayesian approach for model updating and quantification of modeling errors[END_REF], but correlation parameters such as correlation lengths can be included as well. In the PhD work, we will propose an alternative identification approach that estimates the modeling error along the inversion process.

Remark . When a large amount of experimental data is available, the posterior pdf can be asymptotically approximated by a Gaussian pdf centered at the maximum a posteriori (MAP) point p MAP [START_REF] Papadimitriou | Asymptotic expansions for reliability and moments of uncertain systems[END_REF]. It is characterized by a covariance matrix ΣΣ p such that ΣΣ 1 p = r 2 p F MAP|p M AP with F MAP = log(⇡(p|d obs )). This (Laplace) approximation can be used as a cost-effective alternative to computationally demanding methods (e.g.

MCMC sampling) for uncertainty quantification and propagation.

Link between deterministic and stochastic approaches

Considering the Bayesian inference framework with Gaussian multivariate distributions for prior and likelihood pdfs, the previously described deterministic and stochastic approaches for inverse problems can be linked [Tarantola, 2005]. Indeed, prior and likelihood pdfs are written in this context as:

⇡(p)=C 0 .e 1 2 (p p0) T Σ Σ 1 0 (p p0) ⇡(d obs |p) / ⇡ obs (d obs d true ).⇡ mod (d true d(u(p))|p)=C 1 .e 1 2 (d obs d(u(p))) T [Σ Σ obs +Σ Σ mod ] 1 (d obs d(u(p))) (1.5)
where ΣΣ 0 , ΣΣ obs , and ΣΣ mod are covariance matrices related to a priori knowledge on parameters (with mean p 0 ), additive measurement noise on exact observations d true (with zero mean), and modeling uncertainties at measure-ment points (with mean d(u(p))), respectively. This latter information might be given from estimates on discretization or modeling errors [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF], even though it is most often neglected as very little or no quantitative information is at hand. C 0 and C 1 are normalization constants of the form 1/ p (2⇡) n |ΣΣ | with ΣΣ and n the associated covariance matrix and size of the vector variable, respectively.

Remark . Uncorrelated Gaussian distributions (with diagonal covariance matrices) are usually selected, but this is not valid in some cases when errors show correlations. A typical example for observations is densely populated sensor grids or full-field measurements (see later), or in the presence of a systematic error component (faulty measurement equipment/setup).

Using Bayesian inference, the previous framework leads to a posterior pdf of the form ⇡(p|d obs )=Ce S(p) with:

S(p)= 1 2 (d obs d(u(p))) T [ΣΣ obs + ΣΣ mod ] 1 (d obs d(u(p))) + 1 2 (p p 0 ) T ΣΣ 1 0 (p p 0 ) (1.6)
and the MAP approach indicates that the optimal value p opt of the parameter set is such that

p opt = argmax ⇡(p|d obs ) = argmin[ log(⇡(p|d obs ))]
It thus leads to the following problem:

min p2P ✓ 1 2 (d obs d(u(p))) T [ΣΣ obs + ΣΣ mod ] 1 (d obs d(u(p))) + 1 2 (p p 0 ) T ΣΣ 1 0 (p p 0 ) ◆ (1.7)
which turns out to be a deterministic least squares minimization problem with Tikhonov regularization similar to (1.3), and involving Mahalanobis distances with inverse covariance matrices as weights (giving more weight to reliable information). Therefore, as such, the deterministic counterpart of the Bayesian inference scheme incorporates regularization and choice of norm in a natural way, without any heuristic procedure.

When the model is linear with respect to parameters p, that is d(u(p)) = Dp, it is straightforward to observe from (1.6) that the posterior ⇡(p|d obs ) is a Gaussian pdf with covariance matrix

ΣΣ p = ⇥ D T [ΣΣ obs + ΣΣ mod ] 1 D + ΣΣ 1 0 ⇤ 1 ,
and that the optimal parameter value (mean or MAP estimate) reads:

p opt = ⇥ D T [ΣΣ obs + ΣΣ mod ] 1 D + ΣΣ 1 0 ⇤ 1 D T [ΣΣ obs + ΣΣ mod ] 1 d obs + ΣΣ 1 0 p 0 (1.8)
It is thus obtained as a weighted average of the data set d obs and the prior knowledge p 0 . The expressions obtained in this particular case for the posterior mean value p opt and the covariance matrix ΣΣ p are equivalent to the Kalman filter equations with optimal gain [Kalman, 1960].

When the model is nonlinear, an iterative optimization scheme needs to be used to solve (1.7), in a similar way as in Section 1.1.1. Introducing again the Jacobian or sensitivity matrix S p = @d(u(p)) @p , the use of a Gauss-Newton method (for instance) yields:

p (k+1) opt = p (k) opt H (k) 1 ⇣ S (k) T p [ΣΣ obs + ΣΣ mod ] 1 (d(u(p (k) )) d obs )+ΣΣ 1 0 (p (k) p 0 ) ⌘ (1.9) with H (k) = S (k) T p [ΣΣ obs + ΣΣ mod ] 1 S (k)
p + ΣΣ 1 0 the (inexact) Hessian matrix. At convergence, it is fruitful to observe that a perturbation d obs in measurements leads to a variation p opt in the optimal parameter value, such that:

p opt = ⇥ S T p [ΣΣ obs + ΣΣ mod ] 1 S p + ΣΣ 1 0 ⇤ 1 S T p [ΣΣ obs + ΣΣ mod ] 1 d obs (1.10)
Therefore, the covariance matrix ΣΣ popt = h p opt • p T opt i on parameters can be recovered from that of measurement noise ΣΣ obs = h d obs • d T obs i, from other covariance matrices ΣΣ mod and ΣΣ 0 , and from the sensitivity matrix S p evaluated at p = p opt . This piece of information may be used to get a confidence level on identified parameters, or to perform experimental design with minimization of uncertainties on p opt (e.g. minimization of the largest eigenvalue of ΣΣ popt with respect to sensor positions [Bertin et al., 2016a]).

After these previous general considerations, we introduce in the next section a deterministic inversion method based on a specific cost function which is constructed from the Constitutive Relation Error (CRE) concept.

The CRE concept Definition

The energy-based CRE concept has been used for the robust verification of FEM models, that is the a posteriori estimation of discretization error, for more than thirty years. Pioneering ideas can be found in [Ladevèze, 1975, Ladevèze and Leguillon, 1983,Ladevèze and Rougeot, 1997,Destuynder and Métivet, 1999], and a general overview is given in [Ladevèze andPelle, 2005, Ladevèze andChamoin, 2015]. The CRE concept, based on dual analysis, is actually the only way to compute both strict and effective discretization error bounds for linear or nonlinear models of computational mechanics, that can further be used in mesh adaptivity procedures. We detail here the CRE concept in the context of linear elasticity models.

We consider an open bounded domain Ω ⇢ R d , with boundary @Ω, occupied by a linear elastic material (Figure 1.2). We assume that a displacement field u d is prescribed on part @ 1 Ω of the boundary, and that tractions f s d are prescribed on the complementary part @ 2 Ω such that @ 1 Ω \ @ 2 Ω = ; and @ 1 Ω [ @ 2 Ω = @Ω. A body force field f v d may also be given in Ω. Sufficient regularity is assumed for the prescribed data, that is

u d 2 [H 1/2 (@ 1 Ω)] d , f s d 2 [H 1/2 (@ 2 Ω)] d , and f d 2 [H 1 (Ω)] d .
The associated (well-posed) problem is then classically written by splitting in 3 groups of equations:

• kinematic admissibility (defining the space U ad of compatible displacement fields verifying Dirichlet boundary conditions):

u 2 [H 1 (Ω)] d ; u |∂1Ω = u d (1.11)
• static admissibility (defining the space S ad of H(div, Ω) stress fields satisfying equilibrium equations written here in the weak form of the principle of virtual work):

σ 2 [L 2 (Ω)] d(d+1)/2 s ; r • σ 2 [L 2 (Ω)] d ; Z Ω σ : ε(v)= Z Ω f v d • v + Z ∂2Ω f s d • v 8v 2 U 0 ad
(1.12)

• constitutive relation (Hooke's law):

σ = Kε(u) (1.13)
ε(u)=rr S u denotes the infinitesimal strain associated with the displacement u, K is the symmetric positive definite Hooke tensor, and U 0 ad is the vectorial space associated with U ad . A classical primal FE approximation of the problem yields u h 2 U h ad ⇢ U ad (with associated stress field σ h = Kε(u h ) / 2 S ad ) and leads to a discretization error field

! " Ω $ % Ω & % & % : Sensor location ! ' Ω v )
e h = u u h . A measure ke h k 2 K = R Ω ε(e h ):
Kε(e h ) of this error in the energy norm can be defined, and the objective of FE model verification is to compute an a posteriori error estimate on ke h k K . This may be addressed in two ways:

• a primal variational approach, involving the potential energy

J 1 (v)= 1 2 R Ω Kε(v):ε(v) R Ω f v d • v R ∂2Ω f s d • v
and the search space U ad of compatible displacement fields, leads to:

J 1 (u)= inf v2U ad J 1 (v); ke h k 2 K =2(J 1 (u h ) J 1 (u)) 2(J 1 (u h ) J 1 (v)) 8v 2 U ad (1.14)
so that a computable lower error bound on ke h k K can be obtained from a field u ⇤ 2 U ad at disposal (which should live in a larger space than U h ad in order to get a meaningful bound);

• a dual variational approach, involving the complementary energy

J 2 (⌧ ⌧ )= 1 2 R Ω K 1 ⌧ ⌧ : ⌧ ⌧ R ∂1Ω ⌧ ⌧ n • u d and the
search space S ad of equilibrated stress fields, leads to:

J 2 (σ)= inf τ 2S ad J 2 (⌧ ⌧ ); ke h k 2 K =2(J 1 (u h )+J 2 (σ))  2(J 1 (u h )+J 2 (⌧ ⌧ )) 8⌧ ⌧ 2 S ad (1.15)
so that a fully computable (i.e. without any unknown multiplicative constant) upper error bound on ke h k K is obtained from a field σ 2 S ad at disposal. This bound may be used as a guaranteed error estimate for the assessment of accuracy and as a criterion for mesh adaptivity.

Introducing the energy norm k • k K 1 on stress fields, the previous upper bound on ke h k 2 K /2 is written as:

J 1 (u h )+J 2 (σ)= 1 2 Z Ω (σ Kε(u h )) : K 1 (σ Kε(u h )) = 1 2 kσ Kε(u h )k 2 K 1 = E 2 CRE (u h , σ) (1.16)
It is interpreted as a measure of the residual on the constitutive relation for the admissible pair

(u h , σ) 2 U ad ⇥ S ad ;
this is the definition of the CRE functional E CRE .

The bounding property given by E 2 CRE (u h , σ) is also explained from the Prager-Synge theorem [START_REF] Prager | Approximation in elasticity based on the concept of functions spaces[END_REF], that relates the computable CRE term with distances, in energy norms, to the unknown exact solution (u, σ) of (1.11)-(1.13):

ku u h k 2 K + kσ σk 2 K 1 =2.E 2 CRE (u h , σ) (1.17)
The potential of this theorem in the field of error evaluation, even if not originally applied in the FE context, has been known for a long time [Tottenham, 1970, Aubin andBouchard, 1970].

Remark . For any admissible pair (v, ⌧⌧ ) 2 U ad ⇥S ad , the property E CRE (v, ⌧⌧ ) 0 naturally yields, and

E CRE (v, ⌧⌧ )=0
means that (v, ⌧⌧ ) corresponds to the exact solution (u, σ) of the problem. Using the CRE concept, the reference problem can thus be formulated as:

(u, σ) = argmin (v,τ )2U ad ⇥S ad E CRE (v, ⌧⌧ ) (1.18)

Geometrical interpretation

Two geometrical representations of the CRE philosophy are now given (see Figure 1.3). The first one, classical, is in the space of stress fields with inner product hσ 1 , σ 2 i = R Ω σ 1 K 1 σ 2 and associated energy norm. It illustrates the orthogonality property involved in the Prager-Synge theorem. The distance between σ and σ h , that is

p 2E CRE (u h , σ),
is an upper error bound on the discretization error ku u h k K .

The second one, not classical but that will be later reused, is in the space of strain-stress couples s =( ε, σ).

This space is equipped with the energy inner product hs 1 ,s

2 i = R Ω (ε 1 Kε 2 + σ 1 K 1 σ 2 )
and associated energy norm. We denote (A d ) the space of (kinematically and statically) admissible couples, and Γ the space (linear here) associated with the constitutive law. The exact solution of the well-posed problem (1.11)-(1.13) is defined by the intersection between (Γ) and (A d ). It is then easy to show that the value E CRE exactly corresponds to the distance from the solution ŝ =(ε(u h ), σ) 2 (A d ) at hand to (Γ), with orthogonal projection. The stress field σ m obtained after projection is the average field σ m = 1 2 (σ + σ h ). The Prager-Synge theorem reads in this framework:

hs ŝ, s ŝi =2.E 2 CRE (ŝ)

Construction of an equilibrated stress field

The quality of the upper error bound p 2E CRE (u h , σ) depends on that of the statically admissible field σ. The suitable construction of such a fully equilibrated stress field is the key and technical point of the CRE concept. For this purpose, a first approach may consist in using equilibrated elements in a dual version of FEM [Fraeijs de [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF], Fraeijs de Veubeke and Hogge, 1972, Debongnie et al., 1995[START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF][START_REF] De | Upper bounds of the error in local quantities using equilibrated and compatible finite element solutions for linear elastic problems[END_REF][START_REF] De | Moitinho de Almeida[END_REF]. It is the most effective approach, but also the most technical (use of non-conventional FE spaces in the general case, which are not suited to commercial codes) and expensive (another global problem needs to be solved). Other approaches in the literature are based on the post-processing of the approximate FE field σ h using:

• a hybrid flux (or Element Equilibration Technique -EET) technique [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF], Coorevits et al., 1992, Ladevèze and Maunder, 1996, Florentin et al., 2002, Ladevèze et al., 2010a, Pled et al., 2011];

• a flux-free technique [START_REF] Pares | Subdomain-based flux-free a posteriori error estimators[END_REF], Cottereau et al., 2009[START_REF] Gallimard | A constitutive relation error estimator based on traction-free recovery of the equilibrated stress[END_REF], Pares et al., 2009];

• Raviart-Thomas-Nédélec elements over a dual mesh [START_REF] Ern | An accurate h(div) flux reconstruction for discontinuous galerkin approximations of elliptic problems[END_REF], Vohralik, 2007, Vohralik, 2008, Ern and Vohralik, 2010, Vohralik, 2011].

We briefly describe here the hybrid flux (EET) technique; more details will be given in Chapter 4. It is made of two steps:

1. polynomial tractions fK|Γ , equilibrated with the external loading (f v d , f s d ), are built on edges Γ of each element K. They should satisfy fK|Γ = f s d if Γ ⇢ @ 2 Ω, as well as equilibrium at the element level:

Z K f v d • u ⇤ R + Z ∂K fK • u ⇤ R =0 8u ⇤ R 2 U R (K) (1.19)
where U R (K) denotes the space of rigid body motions on K. In practice, tractions are defined as fK|Γ = ⌘ Γ K fΓ , with ⌘ Γ K = ±1 a signed scalar value that ensures continuity of the stress vector across element boundaries, and they are searched as a linear combination of FE shape functions: fK|Γ (x)= P j2JΓ f j K|Γ j (x). J Γ denotes the set of nodes connected to the edge/face Γ.

2. in each element K, a stress field σh|K that satisfies equilibrium:

r • σh + f v d = 0 in K ;σ h n = fK on @K (1.20)
is constructed. The associated local problems are in practice solved with a quasi-explicit technique and polynomial basis [START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in finite element analysis[END_REF], or with a dual approach with degree enrichment (i.e. using higher-order elements).

The construction of fh in the first step leans on the following prolongation (energy) condition:

Z K (σ h σ h )r i = 0 =) Z ∂K fK i = Z K (σ h r i f v d i ) (1.21)
which is enforced for all elements K and all nodes i connected to K; i is the FE shape function associated with node i. This condition, which automatically ensures the equilibration of fK over K (as P i i|K =1 ), leads to the solution to a system of the form:

Rn X r=1 b r Kn (i)=Q Kn (i) with Q Kn (i)= Z Kn (σ h r i f v d i ) br Kn (i)= Z Γr ⌘ Γr K fΓr i (1.22)
over the set of elements K n connected to each node i (patch). R n is the number of edges of the element K n connected to node i. The existence of a solution for the unknowns br Kn (i) (projections of tractions fΓ on FE shape functions) of the system is ensured by the equilibrium property (in the FE sense) verified by σ h , and uniqueness may be obtained minimizing a cost function [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF].

Remark . It can be shown that using the hybrid-flux technique to construct an admissible stress field σh also yields

E CRE (u h , σh )  Cke h k K ,
where C is a constant independent of the mesh size [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. This proves that the CRE-based error estimate has the same convergence rate as the true discretization error.

Extensions to complex models

The CRE concept, primarily used for linear thermal and elasticity problems, can be extended to more complex problems involving a larger class of constitutive laws, using the duality and convex analysis tools developed in [Moreau, 1966, Nayroles, 1973]. It then enables all numerical error sources of FEM simulations to be controlled, which are space or time discretizations, as well as algebraic errors generated by iterative algorithms.

For nonlinear material behavior, such as hyper-elasticity, the key concept is the use of the convex dual potentials and ⇤ that define the material law as σ = ∂ψ ∂ε or ε = ∂ψ ⇤ ∂σ . Potentials and ⇤ are dual in the Legendre-Fenchel sense, i.e.:

⇤ (σ)=sup ε {σ : ε (ε)}
Linear elasticity corresponds to quadratic potentials (ε)= 1 2 ε : Kε and ⇤ (σ)= 1 2 σ : K 1 σ. The definition of the CRE measure then refers to the previous definition of Legendre-Fenchel duality (related to the symmetrized Bergman divergence used in statistics [START_REF] Chen | Metrics defined by bregman divergences[END_REF]) and reads for an admissible pair

(ε, σ) 2 (A d ): E 2 CRE (ε, σ)= Z Ω ( (ε)+ ⇤ (σ) σ :ε) 0 (1.23)
A geometrical interpretation of this error measure is given in Figure 1.4: for a given point (ε, σ), (ε) is the area in blue, ⇤ (σ) is the area in red, and σ :ε is the area in grey. The residual quantity (ε)+ ⇤ (σ) σ :ε is then the remaining blank area. (ê e , êp , ŝ) such that êe + êp = ê. Terms y • x correspond to the duality product between variables x and y. Residuals ⌘ ψ and ⌘ ϕ are local in space and time quantities, so that the global CRE functional reads [Ladevèze, 1998,Ladevèze et al., 1999, Ladevèze, 2001, Ladevèze and Pelle, 2005]:

E 2 CRE|t = Z Ω ⌘ ψ (ê e , ŝ)+ Z t 0 Z Ω ⌘ ϕ ( ėp , ŝ) (1.24)
More details on this functional, with specific application to the elasto-plastic case, will be given in Chapter 5.

Remark . A variant of the literature is to define the CRE measure from the residual on evolution laws alone, enforcing state equations in the definition of admissibility [Ladevèze, 1998]. This is the concept of dissipation error which has a clear mechanical meaning and emphasizes the dissipation properties of the model.

Moreover, an alternative CRE measure denoted Drucker's error can be defined for dynamics problems [START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF], Ladevèze, 1999]; it is based on the Drucker material stability principle [Drucker, 1964].

The CRE concept has been the topic of many studies and applications for FE model verification. Out of early works, several developments have been proposed over the last two decades for various problems such as stochastics [START_REF] Chamoin | Robust goal-oriented error estimation based on the constitutive relation error for stochastic problems[END_REF], transient dynamics [START_REF] Waeytens | Guaranteed error bounds on pointwise quantities of interest for transient viscodynamics problems[END_REF] and vibratory dynamics [START_REF] Wang | Computable upper and lower bounds on eigenfrequencies[END_REF],

or plasticity [START_REF] Ladevèze | Strict upper bounds of the error in calculated outputs of interest for plasticity problems[END_REF]. Applications to several variants of FEM have also been addressed such as XFEM [START_REF] Panetier | Strict and effective bounds in goaloriented error estimation applied to fracture mechanics problems solved with the xfem[END_REF], domain decomposition [START_REF] Parret-Fréaud | Fast estimation of discretization error for fe problems solved by domain decomposition[END_REF], Rey et al., 2014], model reduction [Ladevèze andChamoin, 2011, Chamoin et al., 2017], non-conforming approximations (e.g. Discontinuous

Galerkin) [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous galerkin, and mixed discretizations[END_REF], isogeometric analysis [START_REF] Thai | A posteriori error estimation for isogeometric analysis using the concept of constitutive relation error[END_REF], or multiscale analysis [START_REF] Chamoin | A posteriori error estimation and adaptive strategy for the control of msfem computations[END_REF]. Eventually, coupled with adjoint-based techniques, the CRE concept was effectively used for goaloriented error estimation [START_REF] Chamoin | A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems[END_REF], Ladevèze, 2008, Ladevèze and Chamoin, 2010, Ladevèze et al., 2013, Wang et al., 2016, Chamoin and Legoll, 2021].

Modified CRE -mCRE Definition

During the 1990s, the CRE functional was extended to the framework of inverse problems with overabundant data and uncertain parameters. This energy-based objective functional, with much physical content, is actually a strong alternative to conventional least square approaches for parameter identification or updating. When unknown parameters are constitutive parameters, it is natural to consider cost functions that measure a residual on constitutive equations (which are less reliable than other equations of the problem) from an admissible solution. These equations are thus relaxed (with primal (displacement) and dual (stress) variables treated independently) so that the inverse problem is regularized by a lower constraint from the model. This is in the same philosophy as regularization coming from the discretization of a continuous problem, even though the present CRE regularization is much stronger.

In this context of inverse analysis, the CRE functional is related to Kohn-Vogelius functionals introduced independently in conductivity imaging (electrical impedance tomography) [START_REF] Kohn | Determining conductivity by measurements at boundary[END_REF], Kohn and Lowe, 1988, Kohn and McKenney, 1990, Knowles, 2001]. It leads to a primal-dual formulation of inverse problems, with physics-guided regularization (see [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF] for a mathematical description on these aspects).

For the sake of simplicity, we consider the previous linear elasticity problem with prescribed displacement u d on @ 1 Ω, tractions f s d on @ 2 Ω, and body force f v d in Ω. In addition, observations d obs are made in order to compensate the lack of information on material parameters p defining the Hooke tensor K(p); these observations may be local displacements, strains, load resultants. . . A first approach for using the CRE concept in the context of this ill-posed problem, in order to identify or update parameters p, consists in selecting the classical CRE functional (1.16) as a cost function and integrating experimental data as additional constraints to be fulfilled in admissibility [START_REF] Kohn | A variational method for parameter identification[END_REF], Ladevèze and Reynier, 1989[START_REF] Ladevèze | Error on the constitutive relation in dynamics: theory and application for model updating[END_REF], Rota, 1994, Constantinescu, 1995, Geymonat and Pagano, 2003, Latourte et al., 2008, Florentin and Lubineau, 2010]. In other words, the definition of the admissibility space is revisited, prescribing observations in addition to standard (reliable) kinematic and static information such as equilibrium. Denoting by (A + d ) this enriched admissibility space, the resulting optimization problem reads:

p sol = argmin p2P " min (û,σ)2(A + d ) E 2 CRE (û, σ; p) # with E 2 CRE (û, σ; p)= 1 2 kσ K(p)ε(û)k 2 K 1 (1.25)
It is in practice solved in an iterative process.

The optimal admissible solution (û opt , σopt ) 2 (A + d ) which is obtained at the end of this process satisfies a relation which (in general) is different from the constitutive relation of the model; this mismatch is directly measured by the CRE functional, with E CRE (û opt , σopt ; p sol ) > 0. In the idealistic case where all admissibility conditions, including these associated with measurements, are compatible with the constitutive relation defined by K(p sol ), then E CRE (û opt , σopt ; p sol )=0.

The CRE strategy for inverse problems can be geometrically interpreted as shown in Figure 1.5, as an extension of Figure 1.3; the minimal distance is searched between an admissible solution (û, σ) 2 (A + d ) and the manifold (Γ p ) generated by the parametrized constitutive model. In most situations, (A + d ) \ (Γ p )=;.

Remark . Duality arguments indicate that the optimal admissible field σ, for given p, derives from a displacement field v, i.e. σ = K(p)ε(v). This property may be advantageously used when solving the problem numerically (see later).

Remark . The CRE functional still leads two dual formulations in the case of ill-posed problems. Introducing any .

solution (û 0 , σ0 ) 2 (A + d ) and noticing that R Ω (σ σ0 ):ε(u u 0 )=0, the CRE functional is written as E 2 CRE (û, σ; p)= J 1 (û; p)+J 2 (σ; p) with J 1 (û; p)= 1 2 Z Ω K(p)ε(û):ε(û) Z Ω σ0 : ε(û) 1 2 Z Ω σ0 : ε(û 0 ) J 2 (σ; p)= 1 2 Z Ω K 1 (p)σ :σ Z Ω σ : ε(û 0 ) 1 2 Z Ω σ0 : ε(û 0 ) (1.26)
Remark . In the specific case where both Dirichlet (ũ d ) and Neumann ( f s d ) conditions are prescribed on the whole boundary @Ω, through the enriched admissibility space (A + d ) with over-specified boundary data, one defines ûD and ûN that satisfy the well-posed problems:

r • (K(p)ε(û D )) = r • (K(p)ε(û N )) = f v d ; ûD |∂Ω = ũd ; K(p)ε(û N ) |∂Ω = f s d (1.27)
Then, it is easy to show that:

E 2 CRE (û D , K(p)ε(û N ); p)= 1 2 Z Ω (ε(û N ) ε(û D )) : K(p)(ε(û N ) ε(û D )) = 1 2 Z ∂Ω (K(p)ε(û D )n f s d ) • (ũ d ûN ) (1.28)
which again illustrates that the CRE functional informs on incompatibilities between prescribed conditions and the constitutive relation.

In the original procedure (described above) that uses the CRE functional for inverse problems, strongly enforcing measured data as admissibility constraints is often inappropriate as these data are polluted with measurement noise, if not corrupted. A more flexible and effective approach is obtained with a modified version of the CRE functional, called modified CRE (mCRE) and initially studied for model updating in vibration problems [Ladevèze et al., 1994a, Chouaki et al., 1996, Ladevèze and Chouaki, 1999]. This version is associated with a general framework in which a distinction is made between:

• reliable information on the inverse problem, such as equilibrium equations, location of sensors, or known boundary conditions;

• less reliable information, such as constitutive relations, measurement values, or imperfectly known boundary conditions (when applicable, e.g. see [Feissel andAllix, 2007, Diaz et al., 2015]).

Only reliable information is enforced through an admissibility space denoted (A d ), e.g. using Lagrange multipliers, while other information is relaxed and verified at best when minimizing an appropriate residual functional E 2 mCRE . The inverse procedure thus reads:

p sol = argmin p2P " min (û,σ)2(A d ) E 2 mCRE (û, σ; p) # (1.29)
where the mCRE functional is defined as (considering here that only the constitutive relation and measurement values are relaxed):

E 2 mCRE (û, σ; p)=E 2 CRE (û, σ; p)+ ↵ 2 (d(û) d obs ) T G 1 obs (d(û) d obs ) (1.30)
G obs is a scaling matrix similar to that introduced in Section 1.1.1, while ↵ is a positive scalar factor that has here the dimension of an energy. Therefore, the mCRE cost function appears as a weighted additive combination between the CRE functional (that corresponds to a distance to the constitutive model) and a quadratic error term (that corresponds to a distance to experimental data). As seen later, modeling error and measurement noise can naturally be taken into account in this framework.

Through its formulation, the mCRE approach thus prescribes measurements in a weak manner through penalization. It can also be seen as a classical least squares formulation regularized by physical information on the model (with regularization parameter 1/↵).

Remark . In [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF], it is highlighted that the state constraint in a deterministic least-squares formulation can be included through a penalized potential energy J 1 , thus minimizing

1 2 (d(û) d obs ) T G 1 obs (d(û) d obs )+ 1 α J 1 (û).
As the minimal value of J 1 is problem-dependent, making this approach somewhat impractical, it is proposed to define the penalty term as the sum J 1 + J 2 of potential and complementary energies (with problemindependent minimal value). This comes down to the definition of the mCRE functional.

The optimal admissible solution (û opt , σopt ) obtained at the end of the process represents a compromise between modeling and experimental information, which may also be used for data completion in inverse Cauchy problems [START_REF] Andrieux | An energy error-based method for the resolution of the cauchy problem in 3d linear elasticity[END_REF]. This compromise between satisfying the less reliable modeling information (particularly the constitutive relation) and matching the measured data is performed through the value of the scalar weight ↵. Limiting values correspond to the solution of a classical least-squares minimization (when ↵ ! 0, with emphasis put on satisfying the constitutive relation) and of a pure CRE minimization (when ↵ !1 , with emphasis put on reducing the discrepancy with data). The influence of ↵ on the sensitivity with respect to measurement uncertainties, and therefore on the quality of the inversion performed using mCRE, was illustrated in [START_REF] Deraemaeker | Model validation in the presence of uncertain experimental data[END_REF], Feissel and Allix, 2007, Banerjee et al., 2013]. In any case, the value ↵ should generally be set in regards to the a priori knowledge on model and measurements (e.g. using the Morozov principle or the L-curve method) which is a fundamental aspect of inverse problems participating in their regularization; we will come to this point in

Properties and applications

Several advantages and attractive properties have been reported in the literature when using the mCRE functional.

We may list:

• improved convexity (empirically observed) over its least-squares counterpart, and much less sensitivity to initial guess, with experimental evidence [START_REF] Bonnet | Inverse problems in elasticity[END_REF], Feissel and Allix, 2007, Hadj-Sassi, 2007] or mathematical proof [START_REF] Aquino | Analysis of the error in constitutive equation approach for time-harmonic elasticity imaging[END_REF];

• excellent capability for spatial localization of errors (to localize defects or erroneous sensors for instance) [START_REF] Bonnet | Structural parameter identification using nonlinear gaussian inversion[END_REF],Bui and Constantinescu, 2000,Banerjee et al., 2013,Ben Azzouna et al., 2015,Bonnet and Aquino, 2015, Barbarella et al., 2016];

• robustness to measurement noise [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF], Feissel and Allix, 2007, Nguyen et al., 2008, Ben Azzouna et al., 2015];

• tolerance to incompletely specified boundary conditions [START_REF] Bonnet | Three-dimensional transient elastodynamic inversion using an error in constitutive relation functional[END_REF], Diaz et al., 2015, Aquino and Bonnet, 2019].

Inheriting from the flexibility of the CRE concept, the mCRE strategy was successfully applied in many situations of inverse problems, involving:

• forced vibrations dynamics [Ladevèze et al., 1994a,Chouaki et al., 1996,Ladevèze and Chouaki, 1999,Deraemaeker et al., 2002, Barthe et al., 2004, Deraemaeker et al., 2004, Ladevèze et al., 2006, Faverjon and Sinou, 2008, Banerjee et al., 2013, Charbonnel et al., 2013, Diaz et al., 2015, Guchhait and Banerjee, 2016, Silva and Maia, 2017, Guchhait and Banerjee, 2018];

• transient dynamics [START_REF] Allix | A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues[END_REF], Allix et al., 2005, Feissel and Allix, 2007, Nguyen et al., 2008, Bonnet and Aquino, 2015];

• acoustics [START_REF] Decouvreur | Updating 2d acoustic models with the constitutive relation error method[END_REF], Decouvreur et al., 2007, Decouvreur et al., 2008, Warner et al., 2014];

• nonlinear material behaviors [START_REF] Chouaki | Updating structural dynamic models with emphasis on the damping properties[END_REF], Chouaki et al., 2000, Nguyen et al., 2008, Marchand et al., 2019];

• probabilistic models (when considering a family of similar structures) [START_REF] Ladevèze | Validation of structural dynamics models containing uncertainties[END_REF], Faverjon et al., 2009];

• joint or connexion parameters [START_REF] Gant | Definition and updating of simplified models of joint stiffness[END_REF], Oliveira et al., 2016, Oliveira et al., 2020];

• corrupted measurements [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF]Allix, 2007];

• measurements from imaging [Ben Azzouna et al., 2015, Huang et al., 2016, Ghosh et al., 2017];

• in-situ measurements [START_REF] Charbonnel | A robust cre-based approach for model updating using in situ measurements[END_REF], Bouclier et al., 2013].

In recent years, a goal-oriented variant of the mCRE approach was also proposed [START_REF] Chamoin | Goal-oriented updating of mechanical models using the adjoint framework[END_REF], Djatouti et al., 2020, Djatouti et al., 2021], and the mCRE was used with Kalman filtering for sequential data assimilation [START_REF] Alarcon | A coupled unscented kalman filter and modified error in constitutive relation technique for structural dynamics identification[END_REF], Alarcon, 2012, Marchand et al., 2016] with evolving model parameters such as structural damage.

Simple 1D illustration

As a simple illustration of the mCRE philosophy, let us consider a bar of length L and cross section S. It is clamped at x =0 , with prescribed traction loading F d (assumed to be known) and measured displacement u obs at end point x = L (Figure 1.7). The goal is to identify Young's modulus E. In this simple case, the mCRE functional reads:

E 2 mCRE (û, N ; E)= 1 2 Z L 0 1 ES ( N ES dû dx ) 2 + ↵ 2 (û(L) u obs ) 2 2 u (1.31)
with 2 u a normalization factor (e.g. variance of the measurement noise), and (û, N ) an admissible axial displacementforce pair. The admissibility constraints are:

û(0) = 0 ; d N dx =0on ]0,L[; N (L)=F d (1.32)
so that the admissible force field reads N (x)=F d and is unique.

Introducing = ↵/ 2
u , the minimization of the mCRE functional with respect to û (for fixed E) reads:

Z L 0 (ES dû dx F d ) dû ⇤ dx + (û(L) u obs )û ⇤ (L)=0 8û ⇤ such that û⇤ (0) = 0 (1.33)
and integration by parts leads to the following conditions:

d 2 û dx 2 =0on ]0,L[; ES dû dx (L) F d + (û(L) u obs )=0 (1.34)
Consequently, the optimal admissible displacement field reads:

û(x)= F d + u obs ES + L x (1.35)
We recover that this field corresponds to (F d /ES)x when ↵ ! 0 (classical least-squares formulation) and (u obs /L)x when ↵ !1(pure CRE formulation). In any case, the solution of the inverse problem is

E sol = F d L/Su obs ,
independently of ↵ and with E mCRE (û opt , Nopt ; E sol )=0in this simple case.

Practical implementation

The minimization (1.29) of the mCRE functional is in practice performed from an iterative two-steps algorithm (alternating direction strategy of block Gauss-Seidel type), in which optimal admissible fields are first computed (for fixed p), before minimizing the obtained mCRE cost function with respect to model parameters (for fixed admissible fields). This appears as a relaxation method. Also, a convergence criterion is defined such as a threshold ✏ on the value of the mCRE cost function. The algorithm is summarized as: 0. Initialize the parameter set p (0) and define the stopping criterion threshold ✏ Iteration loop (iteration n +1)

1. Compute optimal admissible fields (û (n+1) , σ(n+1) ):

(û (n+1) , σ(n+1) ) = argmin (û,σ)2(A d ) E 2 mCRE (û, σ; p (n) ) (1.36)
2. Update model parameters:

p (n+1) = argmin p2P F mCRE (p) with F mCRE (p)=E 2 mCRE (û (n+1) , σ(n+1) ; p) (1.37) 3. Stop if F mCRE (p (n+1) )  ✏.
Otherwise, increment n and go to Step 1.

Remark . With the previous iterative scheme, the two partial and complementary minimization steps are made independent. Consequently, the cost function used to update model parameters in Step 2 may be defined in alternative ways, such as from the CRE term E 2 CRE (û (n+1) , σ(n+1) ; p) alone [Feissel andAllix, 2007,Nguyen et al., 2008] or from the Fröbenius norm on constitutive discrepancy [START_REF] Guchhait | Anisotropic linear elastic parameter estimation using error in the constitutive equation functional[END_REF]. Nevertheless, such approaches do not allow for the direct evaluation of gradients of the cost function any longer (see below).

When p contains a large number of parameters (e.g. when describing a parameter field [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF]), the spatial distribution of the cost function F mCRE (p) may be used to select and update only parameters that contribute most to the mismatch [START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF], in addition to detecting corrupted sensors. This is the so-called localization step, performed at the end of Step 1. The associated hierarchical updating, with correction of parameters in zones with high local error alone (similarly to mesh adaptation procedures), is associated with a minimization problem in a lower-dimension parameter space. It thus reduces the computational cost and participates in the regularization process by guiding the identification/updating (since naturally favoring an optimal configuration close to the initial one). Zones may correspond to sub-structures in engineering applications, or to finite elements when identifying a parameter field.

We now describe the numerical implementation of the two minimization steps. We still consider an elasticity problem with constitutive materials to be identified or updated. We also assume that out of overabundant measured data d obs , the prescribed loading and boundary conditions are reliable and correspond to these of the well-posed problem (1.11)-(1.13), so that (A d )=(A d ); alternative situations will be considered in Chapter 2.

In Step 1 of the iterative algorithm, the problem to be solved reads:

min (û,σ)2(A d ) E 2 mCRE (û, σ; p) (1.38)
and thus consists in finding optimal admissible fields for given p (= p (n) at iteration n +1). The associated constrained minimization is conducted by introducing a Lagrangian functional

L(û, σ, ; p)=E 2 mCRE (û, σ; p) Z Ω σ : ε( ) Z Ω f v d • Z ∂2Ω f s d • (1.39)
in order to enforce the static admissibility constraint on σ (while the kinematic admissibility constraint on û is directly included in the search space here). is a Lagrange multiplier field chosen in U 0 ad in order not to reveal unknown reaction forces on @ 1 Ω. As such, the approach is a mixed one; however it is possible from dualization to derive σ from a displacement field v, under the form σ = K(p)ε(v) [START_REF] Ladevèze | Error on the constitutive relation in dynamics: theory and application for model updating[END_REF].

After discretization in space using the FEM (neglecting discretization error for the moment), and using the definition (1.30) of the mCRE, the Lagrangian functional reads:

L h ( Û, V, Λ; p)= 1 2 ( Û V) T K(p)( Û V)+ ↵ 2 (Π Û U obs ) T G obs (Π Û U obs ) Λ T (K(p) V F) (1.40)
with Û, V, and Λ the vectors of nodal values for û, v, and , respectively. K and F are the global stiffness matrix and global load vector, respectively, while Π is a Boolean operator that extracts observation outputs from the displacement field. As û should be kinematically admissible, we introduce the notation Û =[ U d , Ûa ] T with U d the vector of prescribed dofs (obtained from the discretization of u d on @ 1 Ω), and Ûa the complementary vector of active (or free) dofs. Similarly, as 2 U 0 ad , we write Λ =[0, Λ a ] T .

The first-order Karush-Kuhn-Tucker necessary optimality conditions, by searching the saddle-point of L h , read:

ÛT ⇣ K(p)( Û V)+↵Π T G obs (Π Û U obs ) ⌘ =0 8 Û =[0, Ûa ] T VT ⇣ K(p)( V Û) K(p)Λ ⌘ =0 8 V Λ T (K(p) V F)=0 8 Λ =[0, Λ a ] T (1.41)
They lead to V = Û + Λ (+ a rigid body displacement which is taken as 0 here) and to the solution of a coupled linear system of the form:

2 6 4 ↵(Π T G obs Π) aa K aa (p) K aa (p) K aa (p) 3 7 5 0 B @ Ûa Λ a 1 C A = 0 B @ ↵(Π T G obs ) ao U obs ↵(Π T G obs Π) ad U d F a K ad (p)U d 1 C A (1.42)
where we used the block notation

K = 2 6 4 K dd K da K ad K aa 3 7 5 = 2 6 4 K do K ao 3 7 5 ; F = 2 6 4 F d F a 3 7 5 (1.43)
Re-organizing this system yields the following equation for Ûa :

[K aa (p)+↵(Π T G obs Π) aa ] Ûa = F a + ↵(Π T G obs ) ao U obs K ad (p)U d ↵(Π T G obs Π) ad U d (1.44)
This again shows that the mCRE leads to a predicted state, in terms of optimal admissible field Û, that is recovered from both model and measurements. It thus defines a hybrid construction with data-based model enrichment (i.e. correction of the model ignorance) similarly to other variational data assimilation approaches such as the 3D-VAR [START_REF] Law | Data Assimilation[END_REF], the PBDW method [Maday et al., 2015a[START_REF] Maday | Pbdw state estimation: Noisy observations; configuration-adaptive background spaces; physical interpretations[END_REF], the sparse-PGD [START_REF] Ibanez | A multi-dimensional data-driven sparse identification technique: the sparse proper generalized decomposition[END_REF], Ibanez et al., 2019, Chinesta et al., 2020], or other related approaches [Peherstorfer andWillcox, 2015, Rubio et al., 2019a].

Remark . Properties of the coupled system (1.42), which were mathematically studied [START_REF] Aquino | Analysis of the error in constitutive equation approach for time-harmonic elasticity imaging[END_REF],

play a fundamental role in the qualitative and computational aspects of the mCRE minimization. It was shown that this specific system leads to a unique and stable solution (provided data is abundant enough) even in the case of missing information on boundary conditions for the forward problem.

In addition, optimized numerical methods may be used to solve the coupled system. As an example, a (block)

successive over-relaxation (SOR) technique was used in [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF] in the case of large-scale inverse identification, enabling for the use of existing parallel FE codes with minimal modifications. Also, the Sherman-Morrison-Woodbury formula was used in [START_REF] Marchand | Parameter identification and model updating in the context of nonlinear mechanical behaviors using a unified formulation of the modified constitutive relation error concept[END_REF]; it is of the form

(K+↵Π T Π) 1 = K 1 ↵K 1 Π T (I d + ↵ΠK 1 Π T )
1 ΠK 1 and enables the system to be solved by the factorization of K and of (I d + ↵ΠK 1 Π T ) which is usually a low-rank matrix.

Remark . In case the loading F is uncertain and has to be updated, the mCRE functional should read:

1 2 ( Û V) T K(p)( Û V)+ ↵ 1 2 (Π Û U obs ) T G U obs (Π Û U obs )+ ↵ 2 2 (F F obs ) T G F obs (F F obs ) (1.45)
and stationarity on fields Û, V, Λ and F should be written to solve the saddle-point problem and find optimal fields.

In Step 2 of the iterative algorithm, the problem to be solved reads:

min p2P E 2 mCRE (û, σ; p) or min p2P E h 2 mCRE ( Û, V; p) (discretized version) (1.46)
for fixed admissible fields (û, σ) or vectors ( Û, V) obtained at Step 1. This nonlinear minimization may be advanta-geously conducted using a steepest descent approach, of the form:

p (n+1) = p (n) B n r p E h 2 mCRE|p (n)
This takes advantage of good convexity properties of the mCRE functional, and of the fact that the gradient of this functional is easily computed using the adjoint-state method and available fields obtained at Step 1:

d dp E h 2 mCRE ( Û, V; p) |p (n) = d dp L h ( Û, V, Λ; p) |p (n) = @ Û @p |p (n) @L h @ Û + @ V @p |p (n) @L h @ V + @Λ @p |p (n) @L h @Λ + @L h @p |p (n) = 1 2 ( Û V) T @K(p) @p |p (n) ( Û V) Λ T @K(p) @p |p (n) V = 1 2 ( Û V) T @K(p) @p |p (n) ( Û + V) (1.47)
In practice, a backtracking line search (Armijo-Goldstein rule [Armijo, 1966]) may be used to set the step length n , with B n = n I d , in the gradient algorithm, but other methods (e.g. BFGS [Shanno, 1970] or Gauss-Newton) may also be used to define B n .

Remark . The iterative procedure used to minimize the mCRE functional is multi-query, as a series of similar systems (with different values for p) needs to be solved to define the saddle-point and compute gradients of the cost function.

Therefore, to highly reduce the computational cost, ROM techniques may be used [START_REF] Bouclier | Real-time validation of mechanical models coupling pgd and constitutive relation error[END_REF], Marchand et al., 2016] in order to find explicit representations of the optimal admissible fields with respect to p (and possibly to other quantities ↵, d obs . . . ) over the range of input values. This will be the topic of Chapter 4.

The global minimization process is stopped when

E h 2 mCRE  ✏ = ✏ r E h 2 0 with E h 2 0 = 1 2 U T 0 KU 0 (
reference value obtained from an initial U 0 ) and ✏ r a relative error threshold. Alternatively, the stopping criterion may be based on the stagnation of the mCRE functional along successive iterations; this latter criterion is more appropriate in some cases, particularly when the model is biased (see Chapter 3).

Confidence information

In [START_REF] Charbonnel | A robust cre-based approach for model updating using in situ measurements[END_REF], a confidence interval is given for the obtained identified/updated parameter set p (N ) obtained after stopping the iterative process at iteration N , depending on the threshold ✏ r . Performing a secondorder Taylor expansion of the mCRE functional around p (N ) , under the form

F mCRE (p (N ) + p)=F mCRE (p (N ) )+ p T r p F mCRE|p (N ) + p T r 2 p F mCRE|p (N ) p the inequality F mCRE (p (N ) + p)  ✏ r E h 2 0 is solved.
This provides for possible parameter values inside a confidence interval associated with the chosen threshold isovalue.

At convergence (n = 1) of the mCRE procedure with solution p sol , it is again fruitful to observe how a perturbation d obs in measurements leads to a variation p in the optimal parameter value. Using a steepest descent method, we get:

p = B 1 [ @ @d obs r p E h 2 mCRE|p sol ] d obs = B 1 [ @ @d obs r p E h 2 CRE|p sol + S T p|p sol G obs ] d obs (1.48)
with S T p = @d(û)/@p. In the case of a Gauss-Newton algorithm, it yields

B 1 = H 1 mCRE|p sol =[ H CRE|p sol + S T p|p sol G obs S T p|p sol ]
1 with H CRE and H mCRE the (approximate) Hessian matrices of CRE and mCRE functionals.

From these pieces of information, the covariance matrix ΣΣ p = h p • p T i on identified parameters can be recovered from that ΣΣ obs = h d obs • d T obs i of measurement noise and other computable sensitivity quantities.

Stochastic interpretation

A first attempt to draw a parallel between the CRE concept and the general probabilistic inverse problem theory can be found in [START_REF] Bonnet | Structural parameter identification using nonlinear gaussian inversion[END_REF], in which the local value of the CRE functional is used to evaluate the relative variances on model parameters. A more general stochastic interpretation of mCRE was proposed in [START_REF] Deraemaeker | Model validation in the presence of uncertain experimental data[END_REF], referring to the Bayesian inference framework introduced in Section 1.1.1. This interpretation, which we give below, will be useful in further developments of the PhD work, particularly to define the overall quality of the model and perform model selection with respect to information contained in experimental data (see Chapters 3 and5).

Since covariance on the modeling error is usually not known, the idea is to integrate modeling error in a different manner into Bayesian inference, in a more global and less strict framework that allows for more flexibility in the model structure. The value of the CRE functional is thus used to globally quantify the confidence on the less reliable parts of the model (constitutive relation in particular), into a pdf

⇡ mod / e E 2 CRE (û,σ;p) α
accounting for modeling error.

The confidence on modeling exponentially decreases when the CRE value E 2 CRE (û, σ; p) increases, with a rate specified by scalar value ↵.

The likelihood pdf defined in (1.5) is then revisited and is now written as:

⇡(d obs |p)=C 1 .e 1 2 (d obs d(û(p))) T Σ Σ 1 obs (d obs d(û(p))) .e E 2 CRE (û,σ;p) α (1.49)
for any admissible solution (û, σ). The MAP principle (with uniform prior ⇡(p)) or the maximum likelihood principle then yield the minimization of a mCRE-type functional:

min p2P ✓ 1 ↵ E 2 CRE (û, σ; p)+ 1 2 (d(û) d obs ) T ΣΣ 1 obs (d(û) d obs ) ◆ (1.50)
We thus clearly observe how the regularization process is performed with mCRE; this is a regularization from physics, imposing strong constraints on admissible fields while relaxing unreliable information. Furthermore, the mCRE metric is here naturally derived from statistics concepts.

Inverse analysis in the context of DIC 1.2.1 Basics on the DIC method for full-field measurement Introduction

The definition and exploitation of image-based full-field measurement methods has been one of the main breakthroughs of the last decades in experimental and computational mechanics. Since the late 20th century, the rapid development of these contactless methods (in which pixels are sensors) has offered rich input information for model identification and updating from inverse analysis. They permit the design of rich heterogeneous tests (in terms of material, specimen geometry or loading) in order to characterize behaviors of complex material such as composites [START_REF] Périé | Digital image correlation and biaxial test on composite material for anisotropic damage law identification[END_REF]. They thus progressively replaced traditional time-consuming procedures in with multiple homogeneous tests were considered. Among all optical techniques that are available to retrieve full-field experimental data (photoelasticity, interferometry, holography, topography. . . ), digital image correlation (DIC) is an attractive technique as it provides large quantitative experimental information in the form of displacement/strain fields. Typically 1,000 to 10,000 independent measurement points are obtained in 2D from 16-bit images of 1 to 100-Mpixel definition (2-200 MB data sets).

Initiated in the early 80s [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF], Burt et al., 1982, Peters and Ranson, 1982, Sutton et al., 1983, Chu et al., 1985, Bruck et al., 1989], DIC has been widely studied, developed, and applied during the last two decades using modern devices such as high-speed CCD cameras [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF],Sutton et al., 2009,Hild and Roux, 2012b, Sutton, 2013, Sutton and Hild, 2015]. It has now become a common tool in mechanical engineering and materials science. It has also been extended to stereo correlation (SC) that use several cameras [START_REF] Bay | Digital volume correlation: three-dimensional strain mapping using x-ray tomography[END_REF], Garcia et al., 2002, Dufour et al., 2015, Pierré et al., 2016], as well as digital volume correlation (DVC) using X-ray tomography or magnetic resonance imaging (MRI) [Bay, 2008, Leclerc et al., 2011].

The DIC procedure

The purpose of DIC is to extract the displacement field u DIC (x) that best matches between two images: (i) a reference image f (x); (ii) a deformed image g(x). Each image is a set of pixels (x denotes the pixel coordinates)

with distribution of light intensity or gray level values. These values are generated by a random texture (speckle), are encoded over few bits (typically 8 to 16), and are stored as matrices. Invoking the local gray-level conservation f (x)=g(x + u DIC (x)), and relaxing it in order to tackle ill-posedness due to noise, the field u DIC (x) is determined by minimizing a nonlinear correlation residual. Introducing the covariance matrix

ΣΣ pix =[Σ pix (x 1 , x 2 )
] of gray-level noise, this residual functional is defined as

2 (u)= 1 N pix X x12ω X x22ω (f (x 1 ) g(x 1 + u(x 1 ))) Σ 1 pix (x 1 , x 2 )(f (x 2 ) g(x 2 + u(x 2 ))) (1.51)
with ! the domain occupied by the image which is analyzed, and N pix the associated number of pixels. Assuming that the acquisition noise field is Gaussian and white (i.e. spatially uncorrelated), with zero mean and covariance

Σ pix (x 1 , x 2 )= 2 f . (x 1 , x 2 ) ( f
being the standard deviation of noise on each pixel), and that this noise affects each image f and g independently, the previous residual functional reduces to:

2 (u)= 1 2 2 f N pix X x2ω (f (x) g(x + u(x))) 2 (1.52)
and thus involves the L 2 -norm which is optimal for white noise.

Remark . Since x + u(x) may be non-integer, a gray level interpolation is required to evaluate g(x + u(x)). This may be performed using linear or spline interpolation for instance.

Minimizing 2 (u) is still an ill-posed problem as it attempts to find a vector u(x) per scalar quantity f (x) and at every pixel position x. Computational approaches addressing this issue can be split in two categories, depending on the domains on which the gray level residual is defined and the way the regularization is performed.

On the one hand, local subset-based approaches (widely employed in commercial DIC codes) use a local matching procedure on sub-images or zones of interest (ZOI) which are treated independently [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF], Chu et al., 1985,Bay et al., 1999,Bay, 2008,Sutton et al., 2009,Sutton, 2013]. The associated correlation residual reads:

2 L (u)= 1 2 2 f N pix X x2ZOI (f (x) g(x + u(x))) 2 (1.53)
Assuming that u(x) is constant over the ZOI (even though linear or quadratic displacement fields may also be considered), the minimization of (1.53) is equivalent to maximizing the cross-correlation (f ⇤ g)(u)= P x2ZOI f (x)g(x + u) which may be performed using the Fourier transform.

In such local approaches, the finite extension of the ZOI naturally introduces a filtering of the displacement field, whereas the size of the ZOI drives the measurement uncertainties. However, due to the solution of separate problems with no interconnectivity, a drawback of local approaches is the non-continuity of the recovered displacement field.

On the other hand, global approaches minimize a correlation residual defined over the whole measurement zone (region of interest or ROI) [START_REF] Besnard | Finite-element displacement fields analysis from digital images: application to portevin-le châtelier bands[END_REF]:

2 G (u)= 1 2 2 f N pix X x2ROI (f (x) g(x + u(x))) 2 (1.54)
This global formulation results in lower measurement uncertainties compared to local approaches [Hild and Roux, 2012a] as the continuity of the searched displacement field u is imposed throughout the ROI.

In order to solve the optimization problem (1.54), a first-order Taylor expansion g(x+

u(x)) ⇡ g(x)+rf (x)• u(x)
associated with a Gauss-Newton scheme is commonly employed (assuming a small perturbations regime and taking gradients at theoretical convergence) [START_REF] Bruck | Digital image correlation using newton-raphson method of partial differential correction[END_REF]. This leads to an iterative construction of u DIC (x) with successive minimization problems of the form:

u (k) DIC = argmin δu (k) X x2ROI ⇣ f (x) g(x) rf (x) • u (k) (x) ⌘ 2 (1.55)
where g(x)=g(x + u As the minimization of (1.55) is ill-posed and can not be solved as such (the displacement being only detectable along the direction of rf ), global approaches enforce additional regularization [Hild and Roux, 2012b]. This may be performed with global constraint of smoothness on the searched displacement field u DIC (equivalent to a Tikhonov regularization) [START_REF] Horn | Determining optical flow[END_REF], or using mechanics-based regularization techniques that restrict the search of u DIC to a kinematic subspace, with decomposition on a meaningful and a priori chosen kinematic basis.

A typical choice is to use a Rayleigh-Ritz formulation over the ROI. Another typical choice, that we use in the PhD (see Figure 1.8), is to use a Galerkin formulation over a discretized ROI with FE basis [START_REF] Besnard | Finite-element displacement fields analysis from digital images: application to portevin-le châtelier bands[END_REF], i.e. defining u DIC (x)=N(x)U DIC with N gathering the FEM basis functions and U DIC the vector of unknown nodal values. Regular elements are usually considered (i.e. quadrilateral 4-node (Q4) elements for 2D images [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF], and hexahedral 8-node (C8) elements for 3D images) but irregular grids can also be used (e.g. see [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF],Tomicevic et al., 2013] where unstructured meshes with triangular elements (T3) are used).

The successive minimization problems (1.55) thus lead to solutions of linear systems of the form: Figure 1.8: DIC procedure with a global approach (from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]), the measured displacement field being obtained from a FE mesh (that is subsequently used for comparison with FE numerical simulations).

M DIC U (k) DIC = b
Remark . The regularization from a FE basis is a natural approach as it permits a direct comparison with FE simulations for identification purposes. If the same mesh is used for DIC and simulation, no interpolation/projection error is introduced and a seamless coupling is obtained between experimental data and FE codes [Avril et al., 2008a, Leclerc et al., 2009].

Remark . The correlation residual 2 G (u) provides for a good indication on the quality of the recovered displacement u DIC . It is ideally expected to be of the order of 1. Otherwise, 2 G (u) and its spatial distribution may be used to enrich the interpolation kinematic basis with mesh refinement or specific basis functions when a priori relevant mechanical information is available. A typical case is cracked domains for which the enriched XFEM basis is a suitable choice [START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF], Roux et al., 2009, Rannou et al., 2010] as an alternative to mesh refinement [START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF].

In FE-based DIC methods, the definition of the approximation subspace has a direct impact on the accuracy of the estimated field u DIC and subsequent parameter identification (see an illustration in Figure 1.9). A very fine mesh is not necessarily optimal; the number of dofs should indeed be far lower than the number N pix of pixels in the ROI due to the ill-posedness nature of the correlation problem. The smaller the number of pixels per element (defining the spatial resolution expressed in pixels), the larger measurement uncertainties that may be a drawback for further inverse parameter identification. Actually, the image contrast should be sufficiently fine-scaled with respect to the element size in order to obtain a positive definite matrix M DIC ; if not, M DIC is ill-conditioned if not rank-defficient. However, the FE space should be rich enough to accurately represent the unknown kinematics (e.g. strain gradients) [START_REF] Lindner | On the evaluation of stress triaxiality fields in a notched titanium alloy sample via integrated dic[END_REF] and complex geometries [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF].

Consequently, a compromise has to be found (with optimized DIC mesh) in order to manage both measurement uncertainties and accuracy of the displacement interpolation with small correlation residuals, which may not be an easy task [START_REF] Bornert | Assessment of digital image correlation measurement errors: methodology and results[END_REF]. As long as the approximation error due to the mesh remains well below the uncertainty of DIC, a coarse mesh is beneficial. Conversely, or when the picture contrast is poor, finer meshes should be used. We will show in Chapter 3 how this can be automatically addressed from the CRE concept.

Figure 1.9: Fitted Poisson ratio versus mesh size (in pixels) in a biaxially loaded elastic specimen: coarse meshes lead to poor kinematic representations whereas fine meshes lead to high noise sensitivity (from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]).

Remark . The global approach may be associated with a multiscale strategy in order to overcome this compromise between spatial resolution and uncertainty in measurements [Passieux et al., 2015a]. In this strategy, information coming from cameras with different resolutions (farfield/nearfield) is used, as well as DIC meshes at different scales.

In any case, the covariance matrix ΣΣ DIC of the measured dofs due to image noise g can be explicitly computed, so that quantitative information on the signal-to-noise ratio is available. Using (1.56), we directly get U DIC = M 1 DIC P x2ROI gN T rf and consequently:

h U DIC i =0(no systematic bias) ; ΣΣ DIC = h U DIC • U T DIC i =2 2 f M 1 DIC (1.57)
As M 1 DIC is not diagonal, this shows that the measurement noise stemming from FE-DIC is a spatially correlated (not white) displacement noise, with higher noise in zones with low contrast, and for edge nodes due to a reduced connectivity. The spatial correlations thus depend on the image texture and the chosen kinematic basis.

As a result, the output of a FE-based DIC analysis is not only a displacement field over a kinematic basis, but also correlation residuals map and measurement uncertainty quantification in terms of covariance matrix ΣΣ DIC .

These enable to adjust the DIC parameters to make the best out of available images, and to weigh the measured information as it should. This is an essential aspect of the DIC methodology, since measured displacements are almost never the ultimate quantity of interest, but rather intermediate data that are to be further processed to estimate sought mechanical parameters.

Additional regularization

Out of "hard" regularization using a selected kinematic basis and reduced number of dofs, additional mechanicsbased regularization from a priori information on the mechanical behavior may be introduced in the DIC procedure.

In [START_REF] Claire | A finite element formulation to identify damage fields: the equilibrium gap method[END_REF],Roux and Hild, 2008[START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF],Ben Azzouna et al., 2011], elastic regularization involving the balance of momentum equations is used into the DIC algorithm, and it is shown that this acts as a filter removing the unphysical components of the displacement. Mechanical regularization by forcing equilibrium with penalization terms is also used in [START_REF] Tomicevic | Mechanics-aided digital image correlation[END_REF], Naylor et al., 2019] to limit ill-conditioning and therefore mitigate noise sensitivity. In this context, a specific treatment of boundary dofs was proposed [START_REF] Mendoza | Complete mechanical regularization applied to digital image and volume correlation[END_REF], distinguishing the different roles played by Neumann or Dirichlet boundaries in mechanical tests, and proposing a low-pass filtering (with specific characteristic lengths) of surface tractions and boundary layers. Also, a specific filtering of the DIC displacement was analyzed in [START_REF]Estimation of the strain field from full-field displacement noisy data: comparing finite elements global least squares and polynomial diffuse approximation[END_REF],Ben Azzouna et al., 2013] to decrease measurement noise and improve identification. It is based on a diffuse approximation algorithm (moving least-squares approximation), with reconstruction of a smooth field from a cloud of data points.

Identification from DIC displacements

The FEMU method

The measured displacement field u DIC previously obtained from (global) DIC can be used as input data for indirect identification of a set p of constitutive parameters. Among all possible identification methods available with fullfield measurements (see [Avril et al., 2008a] for an overview), the Finite Element Model Updating (FEMU) method is the most popular and widely employed in practical applications [Grédiac, 2004, Kajberg and Lindkvist, 2004, Anghileri et al., 2005,Molimard et al., 2005,Giton et al., 2006,Lecompte et al., 2007,Pagnacco et al., 2013,Rahmani et al., 2013]. In its FEMU-U version, it is a versatile incremental least-squares method that aims at minimizing, by varying the parameter set p, the gap between the displacement field u DIC and its simulated counterpart u s (p) (see Figure 1.10). The field u s (p) is in practice obtained from a FE computation, using a constitutive model parametrized by p and some Dirichlet boundary conditions extracted from DIC. The same mesh is usually used for the numerical construction of u s and the representation of the DIC displacement u DIC , even though this is not compulsory.

Denoting N DIC the total number of kinematic dofs in the DIC mesh, the following weighted least squares function involving the appropriate Mahalanobis distance (originating from the Bayesian framework, see Section 1.1.1) is introduced:

F U (p)= 1 N DIC (U s (p) U DIC ) T ΣΣ 1 DIC (U s (p) U DIC ) (1.58)
U s and U DIC are nodal value vectors associated with fields u s and u DIC , respectively, while

ΣΣ DIC =2 2 f M 1

DIC

is the available covariance matrix of DIC measurement noise defined in the previous section (see (1.57)). The use of ΣΣ DIC enables to transfer the impact of noise coming from images onto the uncertainty of the identified parameters [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF],Gras et al., 2013]. Furthermore, the prefactor 1/N DIC is such that, at convergence, noise in the measured displacement should by itself endow F U (p)=1. The value of F U is thus an intensive quantity that may be used to evaluate the quality of the identification and select an appropriate mechanical model [START_REF] Neggers | Improving full-field identification using progressive model enrichments[END_REF]; any deviation from noise induces a level of the cost function that is on average greater than 1.

Remark . When image information is not available (e.g. when considering synthetic data), a simplified weighting

matrix Σ Σ 1 DIC = 1 2γ 2 f R ROI N T N
may be used in the previous cost function F U (mean field assumption).

Measured reaction forces or applied loads F obs may also be added in the FEMU identification process in order to complement kinematic experimental data. Out of providing additional data, different in nature from the kinematics, they are crucial information when identifying stress-related parameters (as Young's moduli) since purely kinematic information (with no data set on scales) does not permit to provide an absolute characterization of such constitutive parameters; dimensionless quantities such as the ratio of Young's moduli or Poisson's ratio can yet be extracted from purely kinematic data [START_REF] Gras | Identification of the out-of-plane shear modulus of a 3d woven composite[END_REF]. In this framework, it is valuable to compare measured reaction forces F obs with simulated counterparts F s , these being obtained from the parametrized constitutive law and Dirichlet boundary conditions that are provided by DIC measurements. An additional term may thus be added to the FEMU cost function [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF], Bertin et al., 2016a]; it reads:

F F (p)= 1 N F (F s (p) F obs ) T ΣΣ 1 F (F s (p) F obs ) (1.59)
where N F is the number of load cells, and ΣΣ F is the covariance matrix of the measured loads, usually chosen as

2 F I d with 2 F = ⇢ 2 1 |F| 2 + ⇢ 2 0
(sum of load uncertainty, proportional to the magnitude of the load, and measurement uncertainty). The resulting FEMU identification based on the combined displacement fields and reaction forces, referred to as FEMU-UF, then involves the global cost function:

F UF (p)= N DIC N DIC + N F F U (p)+ N F N DIC + N F F F (p) (1.60)
The optimized weighting in the additivity of functionals F U and F F stems from a Bayesian foundation. It is the counterpart of the statistical independence between image and load measurements (which implies that probabilities are multiplied), and it leads to the smallest variance for the estimated parameter set p.

Remark . When F F is minimized alone, it corresponds to a load-based FEMU procedure called FEMU-F [START_REF] Pagnacco | Inverse strategy from displacement field measurement and distributed forces using fea[END_REF].

The solution of the FEMU optimization problem can then be performed with iterative correction on the simulated FE quantities. Considering for instance the FEMU-U functional F U (p) defined in (1.58), a first-order Taylor expansion is introduced:

U s (p (n) + p (n) )=U s (p (n) )+S p p (n) (1.61)
with p (n) the correction of p at iteration n +1and S p = ∂Us ∂p the matrix of sensitivity fields. This matrix may be obtained analytically [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF], even though it is usually approximated in an offline phase using (expensive) finite differences or more elaborated strategies based on reduced order modeling [START_REF] Neggers | Reducing full-field identification cost by using quasi-newton methods[END_REF], Neggers et al., 2018]. From an initial guess p (0) , the Gauss-Newton algorithm then leads to the solution of the following linear system at iteration n +1:

H FEMU p (n) = b (n) FEMU (1.62) with H FEMU = S T p ΣΣ 1 DIC S p = 1 2 2 f S T p M DIC S p (inexact Hessian matrix) b (n) FEMU = S T p ΣΣ 1 DIC (U DIC U s (p (n) ) (1.63)
It should be emphasized that H FEMU , which can be seen as the restriction of M DIC to a subspace generated by sensitivity fields, is thus better conditioned than M DIC .

Remark . The previous iterative algorithm, detailed for the FEMU-U version, can be defined in a similar way for FEMU-F and FEMU-UF versions [Grédiac andHild, 2013, Mathieu et al., 2015]. It would involve the Hessian matrix

H FEMU F = S F T p ΣΣ 1 F S F p with S F p = ∂Fs ∂p .
At convergence, the system (1.62) indicates that the fluctuating part p of the identified parameters due to noise is p = H 1 FEMU S T p ΣΣ 1 DIC U DIC . Therefore, h pi =0(unbiased identification) and the covariance matrix ΣΣ p of the identified parameters due to image noise reads:

ΣΣ p = h p • p T i = H 1 FEMU (1.64)
It indicates the identification uncertainty (associated with a chosen constitutive law), after tracking down the uncertainties from the measurement step (pixel level) down to the identification step (parameter level) [START_REF] Roux | Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields[END_REF]. Several scalar indicators qualifying the robustness of parameter identification and enabling optimal experimental design (by maximizing these indicators) may then be defined, such as: (i) the smallest eigenvalue (or the ratio between smallest and highest eigenvalues) of H FEMU [Bertin et al., 2016a, Chamoin et al., 2020]; (ii) the quantity det(H FEMU ) (inverse of the uncertainty volume). In addition, parameters for which uncertainty is below the measurement noise can be discarded from the identification procedure.

Alternative identification methods

A large overview of available identification techniques using full-field measurements obtained from imaging techniques is given in [Avril et al., 2008a]. Out of the FEMU method, we mention:

• the Equilibrium Gap Method (EGM) [START_REF] Claire | Identification of damage fields using kinematic measurements[END_REF]Hild, 2008] based on the minimization of equilibrium residuals after discretizing the balance equations. It was used for the identification of a distribution of elastic properties and its evolution during damage tests [START_REF] Claire | A finite element formulation to identify damage fields: the equilibrium gap method[END_REF], Claire et al., 2007], for the identification of damage models [START_REF] Périé | Digital image correlation and biaxial test on composite material for anisotropic damage law identification[END_REF], Ben Azzouna et al., 2011] or parameters of orthotropic composite materials [START_REF] Crouzeix | An orthotropic variant of the equilibrium gap method applied to the analysis of a biaxial test on a composite material[END_REF].

• the Virtual Fields Method (VFM) [Grédiac, 1989, Grédiac et al., 2002a, Grédiac et al., 2006, Grédiac and Pierron, 2006, Promma et al., 2009] using virtual works principle and consisting in minimizing the difference between internal and external virtual works by optimizing material parameters. The key point of the method, which requires full-field measurements over the whole domain, is the appropriate choice of virtual fields in order to extract relevant information for identification purposes. The VFM was applied in various contexts, in particular for the characterization of composite materials [START_REF] Pierron | Identification of the through-thickness properties of thick laminates using the virtual fields method[END_REF], Pierron et al., 2007] or damage models [START_REF] Chalal | Experimental identification of a nonlinear model for composites using the grid technique coupled to the virtual fields method[END_REF]. It was also shown [Avril and Pierron, 2007] that the stationarity conditions for cost functionals associated with alternative identification methods (including FEMU or EGM) could be interpreted in terms of the virtual fields method for suitably chosen virtual fields.

• the Reciprocity Gap Method (RGM) [Bui, 1993, Andrieux et al., 1999] relying on the Maxwell-Betti reciprocity relation between the measured field and a fictitious solution field in the absence of the unknown elements to be identified. Valid for linear constitutive laws, it was mainly used for source or crack identification [START_REF] Andrieux | Identification of planar cracks by complete overdetermined data: inversion formulae[END_REF], Bui et al., 2004, Andrieux and Bui, 2006, Andrieux et al., 2013] by minimizing the reciprocity gap. It is a suitable method in situations where mechanical field measurements are available on the full outer boundary of a body.

• the Constitutive Equation Gap Method (CEGM), using similar ideas as the CRE concept. As detailed in Section 1.1.3, a specificity of this method is that it relaxes the constitutive equation (which is an approximate mathematical representation) and uses the associated residual as a guiding norm for the identification. Conversely, reliable information such as equilibrium is strongly enforced. A pure CRE version (with experimental data imposed as constraints) was used for identification from full-field measurements in [START_REF] Geymonat | Identification of elastic parameters by displacement field measurement[END_REF], Geymonat and Pagano, 2003, Florentin and Lubineau, 2010, Florentin and Lubineau, 2011, Pagano and Bonnet, 2013, Moussawi et al., 2013, Lubineau et al., 2015, Babaniyia et al., 2017] for linear problems and in [START_REF] Latourte | Elastoplastic behavior identification for heterogeneous loadings and materials[END_REF],Blaysat et al., 2012] for nonlinear elastoplasticity problems. The modified CRE version was investigated in [START_REF] Calloch | Identification de modèles de comportement de matériaux solides : utilisation d'essais et de calculs[END_REF], Ben Azzouna et al., 2015, Barbarella et al., 2016, Barbarella et al., 2017, Ghosh et al., 2017].

We mention that several of these methods consider strain information, obtained by numerical differentiation of DIC data, which often causes significant amplification of the original measurement errors.

In the next chapters, we will investigate a robust identification procedure based on the mCRE concept.

Integrated approach: I-DIC

Numerical procedure

In the previously described DIC-FEMU two-step procedure, displacement data are only an intermediate quantity, whose computation may imply constraints (e.g. coarse mesh for DIC to be well-conditioned) that are not ideal to FEMU. The idea in integrated DIC (I-DIC) is to merge measurement and identification procedures in a single step with integration of the identification step into the DIC algorithm [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF]. I-DIC thus directly uses the images (rather than a measured displacement field) to determine material parameters, as well as bases dedicated to the mechanical problem which is considered. A specificity of I-DIC is that it allows any mesh size and structure (avoiding discretization errors), as the mesh size does not prevent the convergence of computations; in [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF] for instance, the element size was chosen as a decreasing function of the sensitivity with respect to the studied parameters. I-DIC also allows standard FE codes (e.g. Abaqus) to be used in a non-intrusive way, and in many contexts [START_REF] Roux | Stress intensity factor measurements from digital image correlation : post-processing and integrated approaches[END_REF], Leclerc et al., 2009[START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF], Mathieu et al., 2012, Lindner et al., 2015, Mathieu et al., 2015[START_REF] Bertin | Integrated digital image correlation applied to elastoplastic identification in a biaxial experiment[END_REF], Neggers et al., 2017].

In the I-DIC framework, dofs in image correlation are directly related to the parameters p to be identified (elasticity parameters, stress intensity factors. . . ). A restricted number of dofs is then involved, which is beneficial for the conditioning of the inverse procedure. I-DIC thus consists in minimizing the gap between image frames f and g using a specific kinematics u s (p) driven by a parametrized model, rather than simple standard FE representations.

The parametrized closed-form solutions u s (p), satisfying the constitutive equation and equilibrium, may be obtained analytically, from classical FEM, or using ROM.

The I-DIC functional to be minimized reads:

F IDIC (p)= 1 2 2 f N pix X x2ROI (f (x) g(x + u s (x, p))) 2
(1.65)

and leads (after discretization with u s = NU s ) to an iterative strategy with successive linear systems of the form (at iteration n +1):

H IDIC p (n) = b (n) IDIC
(1.66)

with

H IDIC = X x2ROI S T p N T rf r T f NS p (symmetric semi-positive Hessian matrix) b (n) IDIC = X x2ROI (f g)S T p N T rf (1.67)
These two quantities come down to H IDIC = S T p M DIC S p and b

(n) IDIC = S T p b (n)
DIC when similar meshes are used in DIC-FEMU and I-DIC. This shows that I-DIC merely consists of a projection onto the parameter space through sensitivity fields S p = @U s /@p [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF].

At convergence, and if only random acquisition noise is present, the expectation value of F IDIC should approach unity. Also, the system (1.66) implies that the fluctuating part p of the identified parameters due to image noise is p = H 1 IDIC S T p P x2ROI gN T rf , so that [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF]:

h pi =0 ; ΣΣ p = h p • p T i = H 1 IDIC (1.68)
This result indicates that the procedure is unbiased and provides an estimate on the reliability of the identified parameters.

In theory, if the same mesh is used in DIC and in FEMU (for the computation of sensitivity fields), the two-step DIC-FEMU procedure and the integrated I-DIC procedure should be equivalent. However, this is valid only for small noise level [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF], Ruybalid et al., 2016], when an appropriate metric with ΣΣ DIC is used in FEMU (weighted FEMU), and when the smallest element size chosen in I-DIC (addressing complex specimen geometry, strain concentrations. . . ) is not critical for displacement uncertainty with classical DIC.

Remark . As for FEMU-F and FEMU-UF versions, a load-based term of the form (1.59) may be added in the I-DIC procedure [START_REF] Bertin | Integrated digital image correlation applied to elastoplastic identification in a biaxial experiment[END_REF], Neggers et al., 2017].

Remark . In the case of insignificant or limited sensitivity, the system (1.66) may remain ill-conditioned. In order to circumvent this difficulty, a Tikhonov-type regularization can be used [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite : Uncertainty analysis and regularization[END_REF] with penalization from a set p 0 of reference parameters obtained from other sources (other experiments, expert knowledge).

A flowchart showing the different steps of weighted FEMU and I-DIC (U versions) is given in Figure 1.11.

Figure 1.11: Different steps of weighted FEMU-U and IDIC-U.

Illustrative example

As an illustration of parameter identification from I-DIC, we report here a test-case considered in [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. It deals with the calibration of elastoplastic material parameters from a tensile experiment (cyclic loading) on an aluminum alloy (AA2219) sample. This sample has a dog-bone geometry with a hole that increases the sensitivity to the sought material parameters (Figure 1.12). A 4-Mpixel camera (with 16-bit digitization) is used to get the surface images, with a pixel physical size of 11.14 µm.

In order to drive I-DIC computations, Dirichlet boundary conditions are measured via global DIC with a triangular mesh, before being extruded along the thickness of the sample in order to be prescribed on the 3D computation mesh. Global DIC results are shown in Figure 1.13. Components of the displacement and in-plane strain are plotted in two cross-sections. They enable to observe four strained bands at the end of the test, emanating from the hole and developing over the last 200 loading increments.

Ten material parameters are identified (Young's modulus E, Poisson's ratio ⌫, yield stress y , and seven additional stress increments ∆ i in a non-parametric hardening law). The identification is performed from images and measured load data F obs . Figure 1.14 shows the corresponding sensitivity fields for some of the parameters, computed using finite differences. Figure 1.15 shows various residuals for the initial step of the I-DIC procedure and Figure 1.14: Sensitivity fields for some of the final set of material parameters, from [START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]. The plotted fields are [S i ]p i in order to share the same units and make comparisons.

mCRE framework with full-field measurements

Motivation and preliminary aspects

The previously described DIC-FEMU or I-DIC methods are committed to propose a convenient and robust processing of measurement perturbations which are inherent to experimental tests and can make the identification problem difficult to handle (especially for low-strain elasticity where the signal-to-noise ratio is small). This is a sound concern that has to be kept in any new developments around identification procedures. Nevertheless, and out of discretization error which is assumed to be negligible, the error in the mathematical model itself is set aside in the direct use of FE computations in association with experimental information. Consequently, uncertainty in the model class (particularly in the structure of the constitutive law which is supposed to fit within a given parametric representation) is usually not considered, except in some few works e.g. [START_REF] Hild | Toward 4d mechanical correlation[END_REF]. However, quite often, the a priori choice of the constitutive law is empirical and does not offer sufficient flexibility with regard to available experimental data, leading to model errors that may hinder a proper identification. It thus seems valuable to question the reliability of the constitutive model, relaxing it as well as other uncertain parts of the model.

On the other hand, a feature in DIC-FEMU or IDIC methods is the limitation of the simulation domain to the spatial subregion that is analyzed with cameras, thereby reducing computational costs. Nevertheless, this is also associated with the drawback that boundary conditions, prescribed to compute simulated fields, are extracted from full-field measurements so that measurement uncertainty and high-frequency fluctuations due to noise are directly transferred to the identification algorithm, even though the associated spurious effects may be attenuated with empirical interpolation techniques. A sounder procedure would be to consider available information on external boundary conditions (that may be far for the measurement zone) alone, without adding supplementary assumptions to recover unavailable information for simulation purposes.

We wish here to address these issues by using the mCRE framework that is based on the reliability of information, thus proposing an appropriate metric with natural consideration of all error and uncertainty sources. In this framework, and as described in Chapter 1, available modeling and experimental information on the identification problem is partitioned into reliable and less reliable sets. Reliable information (such as local equilibrium or known conditions on the specimen boundary) is enforced in an admissibility space (A d ), while uncertain information on the model and measurements (constitutive relation, part of boundary conditions, measurement values) is relaxed during identification. In addition, the mCRE philosophy is perfectly suited to situations where interior data is abundant over a subset of the medium alone, with a possible lack of information outside the measurement zone or on the boundary of the specimen (incomplete data).

In the context of identification from full-field measurements, the (modified) CRE concept with metric in a space of constitutive laws has already been investigated (see [START_REF] Pagano | Constitutive equation gap. Full-Field Measurements and Identification in Solid Mechanics[END_REF] for a review). It was mainly used in its pure CRE version (referred to as Constitutive Equation Gap method), with compatibility with measurements strongly enforced in the admissibility space [START_REF] Geymonat | Identification of elastic parameters by displacement field measurement[END_REF], Geymonat and Pagano, 2003, Latourte et al., 2008,Florentin and Lubineau, 2010,Florentin and Lubineau, 2011,Blaysat et al., 2012,Moussawi et al., 2013]. Such a version may not be desirable for the PhD objectives as it presents several drawbacks. On the one hand, and as mentioned in Section 1.1.3, prescribing displacement data exactly implies that information on measurement noise is not taken into account. On the other hand, a pure CRE-based inversion method is fed by strain measurements, obtained in practice by differentiating the displacement measurements; this leads to an increase of the perturbation level and thus requires the use of specific non-mechanical filters [START_REF] Avril | Comparison of two approaches for differentiating full-field data in solid mechanics[END_REF] that may introduce artificial bias in the identification results [Ben Azzouna et al., 2013].

The modified CRE version (mCRE), with physical regularization from the model, was considered in association with full-field measurements in quite few works; we may list [START_REF] Calloch | Identification de modèles de comportement de matériaux solides : utilisation d'essais et de calculs[END_REF], Ben Azzouna et al., 2015, Huang et al., 2016] for elastostatics models, or [Bonnet and[START_REF] Diaz | A modified error in constitutive equation approach for frequency-domain viscoelasticity imaging using interior data[END_REF] for (visco-)elastodynamics models. In many of these works, the attractive properties of the mCRE functional (in terms of improved convexity and robustness with respect to perturbations compared to least-squares approaches) were observed and analyzed.

Here, we go one step further by defining a mCRE procedure which is fully consistent with experimental information and measurement noise; it will be denoted DIC-mCRE when used in combination with DIC measurements, and mI-DIC in an integrated version that is also proposed.

DIC-mCRE method

Reference identification problem

We consider an elastic medium in a domain Ω with boundary @Ω (Figure 2.1). On this boundary @Ω, we distinguish 3 zones:

• @ 1 Ω on which reliable Dirichlet boundary conditions u d are prescribed (with potential measured load);

• @ 2 Ω on which Neumann boundary conditions f s d are prescribed. These may be reliable (e.g. free edge) or not;

• @ ; Ω on which there is no modeling or experimental information.

A body force field f v d may also be applied in Ω, even though it will be zero in the following numerical applications. The objective is to recover material parameters p from overabundant and noisy full-field displacement measurements. We assume they are available in a restricted measurement zone Ω m ⇢ Ω, and that they are obtained from global DIC (for the moment). Additional experimental information may be available such as global loads on

@ 1 Ω [ @ 2 Ω.
The inverse problem is here solved using the mCRE framework, that reads (cf Section 1.1.3):

p sol = argmin p2P " min (û,σ)2(A d ) E 2 mCRE (û, σ; p) # (2.1)
From the above ill-posed problem, appropriate admissibility space (A d ) and mCRE functional E mCRE need to be defined depending on the configuration. The definition of the mCRE functional (with associated admissibility space) proposed below is based on reliability of information, and also takes advantage of a key feature of the DIC methodology, which is the explicit characterization of measured displacement uncertainty through the covariance

matrix ΣΣ DIC =2 2 f M 1 DIC .

Definition of the cost function

Managing the overall uncertainty in the identification process appears as a compromise between empirical modeling (in particular with regards to the constitutive law that describes the material behavior) and noise corrupting the fullfield measurements. When phrased this way, it appears that the FEMU metric between measured and computed displacement fields (see (1.58)) has to be balanced by a metric in the constitutive law space.

In the simple case where information on @ 2 Ω is reliable, with no load measurement, and in agreement with the stochastic framework of Section 1.1.3, we thus define the following mCRE functional as:

E 2 mCRE (û, σ; p)=E 2 CRE (û, σ; p)+ ↵ 2 . 1 N DIC (d(û) U DIC ) T ΣΣ 1 DIC (d(û) U DIC ) with E 2 CRE (û, σ; p)= 1 2 Z Ω (σ K(p)ε(û))K(p) 1 (σ K(p)ε(û)) (2.2)
E CRE (û, σ; p) is the classical CRE functional, and d is an operator that evaluates the admissible displacement field û on DIC measurement points. The distance to observations is defined from the appropriate norm already used in weighted FEMU (Mahalanobis distance with inverse covariance matrix ΣΣ 1 DIC ), which enables to track the impact of measurement noise starting from image acquisition.

In this context, the pair (û, σ) lives in an admissibility space (A d ) in which reliable information is enforced. This space exhibits here the following constraints:

û|∂1Ω = u d r • σ = f v d in Ω ;σn |∂2Ω = f s d or Z Ω σ : ε(v)= Z Ω f v d • v + Z ∂2Ω f s d • v 8v 2 U 0 ad (2.3)
Remark . In order to circumvent mesh projection issues when evaluating û on nodes of the DIC mesh, a possible strategy is to project both û and u DIC = NU DIC on a shared regular fine grid before comparing nodal values on this fine grid. Denoting by P the transfer matrix to go from U DIC to the new vector ŨDIC of nodal measurements on the fine grid, the weighting matrix in the second term of the mCRE functional then reads P T ΣΣ 1 DIC P 1 .

In its FEM discretized version (neglecting discretization error for the moment), the previous mCRE functional reads:

E h 2 mCRE ( Û, V; p)= 1 2 ( Û V) T K(p)( Û V)+ ↵ 2 . 1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) (2.4)
where Π extracts local dofs from the vector Û (defined over the whole domain Ω) to fit the measurement zone Ω m , and projects these extracted dofs from the simulation mesh to the DIC mesh thanks to FE shape functions; Π is a restriction of the identity matrix when the simulation mesh and the DIC mesh coincide locally. Let us remind that

V is statically admissible in the sense that is verifies the discretized FEM equilibrium V T (K(p) V F)=0for any virtual vector V 2 U h,0 ad .

Consequently, the mCRE functional is composed of two parts with quadratic terms associated with each uncertainty source: one part is related to the constitutive relation error, and the other part represents the gap between outputs derived from admissible fields and observed data.

Remark . In the case where part of the loading force is measured (in terms of resultant usually), the previous functional should be generalized as:

E h 2 mCRE ( Û, V; p)= 1 2 ( Û V) T K(p)( Û V) + ↵ 2 . 1 N DIC + N F h (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC )+(Π F F F obs ) T ΣΣ 1 F (Π F F F obs ) i (2.5)
with Π F the operator that extracts measured components or resultants from the global loading vector F, and ΣΣ F the covariance matrix of load measurements. The second part of the mCRE (which may merge heterogeneous data sets) is properly scaled with dimensionless terms.

Also, when part of Neumann boundary conditions is uncertain and not measured, it should not be included in the admissibility space (A d ).

Implementation details

As described in Chapter 1, the minimization (2.1) of the mCRE functional for solving the inverse problem is in practice conducted by using an iterative two-steps algorithm that is summarized below.

0. Initialize the parameter set p (0) and set the stopping criterion threshold Iteration loop (iteration n +1)

1. Compute optimal admissible fields ( Û(n+1) , V(n+1) ) for given p (n) :

( Û(n+1) , V(n+1) ) = argmin ( Û, V)2(A d ) E h 2 mCRE ( Û, V; p (n) ) (2.6)
2. Update model parameters for fixed admissible fields:

p (n+1) = argmin p2P F h mCRE (p) with F h mCRE (p)=E h 2 mCRE ( Û(n+1) , V(n+1) ; p) (2.7) 3. Stop if the stopping criterion is reached ( F h mCRE (p (n+1) )  ✏ r E h 2 0
, or stagnation criterion). Otherwise, increment n and go to Step 1. The optimization problem (minimization with linear constraint) associated with Step 1 is addressed with a Lagrangian functional. It reads (for the case (2.4)):

Remark . The criterion

F h mCRE (p (n+1) )  ✏ r E h
L h ( Û, V, Λ; p (n) )= 1 2 ( Û V) T K(p (n) )( Û V) + ↵ 2 . 1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) Λ T (K(p (n) ) V F) (2.8)
Referring to Section 1.1.3, the search of the saddle-point leads to V(n+1) = Û(n+1) + Λ (n+1) and to the linear system:

K (n) mCRE 0 B @ Û(n+1) a Λ (n+1) a 1 C A = F (n) mCRE (2.9) with K (n) mCRE = 2 6 4 ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa K aa (p (n) ) K aa (p (n) ) K aa (p (n) ) 3 7 5
F (n) mCRE = 0 B @ ↵. 1 N DIC (Π T ΣΣ 1 DIC ) ao U DIC ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) ad U d F a K ad (p (n) )U d 1 C A (2.10)
We again used the notation Û =[ U d , Ûa ] T and Λ =[ 0, Λ a ] T to take prescribed dofs into account. Re-organizing this system yields the following equation for

Û(n+1) a : K (n) aa Û(n+1) a = F (n) a (2.11) with K (n) aa = K aa (p (n) )+↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa F (n) a = F a + ↵. 1 N DIC (Π T ΣΣ 1 DIC ) ao U obs K ad (p (n) )U d ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) ad U d (2.
12)

It is important to notice again that in Step 1, optimal admissible fields are defined from both modeling and experimental information available, with a compromise on the reliability of each.

Remark . The symmetric semi-positive matrix of the linear system (2.9) may be a rank-deficient matrix when Ω m is much smaller than Ω. In such a case, the system is solved through a QR decomposition with fixed truncation level (based on rank-deficiency), as proposed in [Ben Azzouna et al., 2015]. We thus get a particular solution to the system, which is sufficient as any solution in the matrix kernel does not impact the magnitude of the mCRE functional.

Also, the system to be solved is twice the size of a classical FE system, but it yields gradients of the cost function at once (see below).

As regards the nonlinear minimization in Step 2, it is performed using a steepest descent approach:

p (n+1) = p (n) B n r p E h 2 mCRE|p (n)
with gradient computed using the adjoint-state method and available admissible fields obtained at Step 1:

r p E h 2 mCRE|p (n) = 1 2 ( Û(n+1) V(n+1) ) T @K(p) @p |p (n) ( Û(n+1) + V(n+1) ) (2.13)
The overall DIC-mCRE algorithm is shown in Figure 2.2.

Setting of the weight parameter

The scalar ↵, which has here the dimension of an energy, is a parameter of the mCRE method that can be tuned to play on the relative weight of the two parts (CRE term, distance to measurements) of the functional. It achieves a balance between minimizing the CRE and the measurement misfit. Out of the specific continuation scheme which was studied in the context of mCRE [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], several approaches related to regularization techniques can be employed to define an optimized value of ↵. In particular, the impact of ↵ on the robustness with respect to measurement noise was thoroughly studied [START_REF] Barthe | Validation and updating of industrial models based on the constitutive relation error[END_REF], Deraemaeker et al., 2004, Warner et al., 2014]; it was empirically found that an adequate weight, minimizing the effect of noise on the computation of the CRE, could be obtained by considering the same confidence in the model and in the measurements. After normalization, this comes down to choosing ↵ such that the two terms of the mCRE functional have the same order of magnitude. This strategy is closely related to the L-curve technique [Miller, 1970, Calvetti et al., 2000], which was investigated for mCRE [START_REF] Huang | Modified constitutive relation error: An identification framework dealing with the reliability of information[END_REF].

Another technique, which we propose to use here and which takes advantage of the a priori knowledge on the noise level, is the Morozov discrepancy principle [Morozov, 1968, Morozov, 1984, Nair et al., 2003, Isakov, 2006].

This principle has a clear physical sense and is actually in the philosophy of reliability of information. It comes down to setting ↵ such that the distance to measurements falls within the noise level. It is indeed meaningless trying to decrease the distance to measurements below this level, as this corresponds to the experimental uncertainty range.

With the covariance-based weighted norm which is used and the employed prefactors, the term of distance to measurements in the new definition (2.4) or (2.5) of the mCRE functional appears as an intensive quantity; satisfying at convergence the Morozov discrepancy principle thus reads:

• when considering (2.4):

1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) ⇡ 1 (2.14)
• when considering (2.5):

1 N DIC + N F h (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC )+(Π F F F obs ) T ΣΣ 1 F (Π F F F obs ) i ⇡ 2 (2.15)
Consequently, in the following, the coefficient ↵ is adjusted as the smallest positive value that enables to reach the previous criterion, ensuring a final discrepancy between computed and measured quantities which is at the level of the noise. Starting from small values of ↵, a bisection method or a progressive increase is implemented to catch this value. The procedure for setting ↵ may be performed after reaching convergence of the iterative inversion process, or along the iterations (during Step 1).

Remark . For small values of ↵, it is easier (i.e. faster convergence) to find a corresponding minimizer for the mCRE optimization problem as there is less constraint from data. Also, the minimizer ( Û(j) , V(j) , p

sol ) obtained with ↵ (j) can be used as an improved initial guess for the optimization problem with the next scalar weight ↵ (j+1) .

Confidence on identified parameters

As mentioned in Chapter 1, at convergence (n = 1) of the mCRE procedure with solution p sol , it is useful to observe how a perturbation U DIC in full-field measurements leads to a variation p in the optimal value of the parameter set. From the steepest descent algorithm, we get:

p = B 1 [ @ @U DIC r p E h 2 mCRE|p sol ] U DIC = B 1 [ @ @U DIC r p E h 2 CRE|p sol + S T p|p sol Π T ΣΣ 1 DIC ] U DIC (2.16)
with S T p = @ Û/@p. From these pieces of information, the covariance matrix ΣΣ p = h p• p T i on identified parameters can be recovered from that ΣΣ DIC = h U DIC • U T DIC i on measurement noise and other computable sensitivity quantities. Therefore, the propagation of uncertainties from image acquisition (i.e. from pixels and measured FE displacement fields) to the sought material parameters can again be performed in the mCRE-based inversion technique.

Integrated version -mI-DIC

In a similar manner as for the (FEMU-integrated) I-DIC method or integrated versions of other identification methods (e.g. EGM [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis of fractured samples[END_REF]), an integrated version of the mCRE functional (denoted mI-DIC) is proposed here when associated with image-based information. This again circumvents finding a compromise to manage both measurement uncertainties and kinematic accuracy when designing a DIC mesh. Following the philosophy of I-DIC (see Section 1.2.3), the corresponding mI-DIC functional is naturally defined as (in its discretized version):

E h 2 mIDIC ( Û, V; p)= 1 2 ( Û V) T K(p)( Û V)+ ↵ 2 . 1 2 2 f N pix X x2ROI ⇣ f (x) g(x + N(x) Û) ⌘ 2 (2.17)
It is again minimized in a sequential two-steps procedure. At iteration n +1 , optimal admissible fields are computed in a first step by solving:

( Û(n+1) , V(n+1) ) = argmin ( Û, V)2(A d ) E h 2 mIDIC ( Û, V; p (n) ) (2.18)
Introducing the associated Lagrangian functional:

L h ( Û, V, Λ; p (n) )= 1 2 ( Û V) T K(p (n) )( Û V) + ↵ 2 . 1 2 2 f N pix X x2ROI ⇣ f (x) g(x + N(x) Û) ⌘ 2 Λ T (K(p (n) )V F) (2.19)
and linearizing g(x + N(x) Û), the search of the saddle-point consists in finding Û(n+1) = Û(n) + Û(n) and Λ (n+1)

(and V(n+1) = Û(n+1) + Λ (n+1) ) solution of the linear system:

K (n) mIDIC 0 B @ Û(n+1) a Λ (n+1) a 1 C A = F (n) mIDIC (2.20) with K (n) mIDIC = 2 6 4 ↵. 1 2γ 2 f Npix [ P x2ROI N T rf r T f N] aa K aa (p (n) ) K aa (p (n) ) K aa (p (n) ) 3 7 5
F (n) mIDIC = 0 B @ ↵. 1 2γ 2 f Npix P x2ROI (N T rf ) a (f g) ↵. 1 2γ 2 f Npix [ P x2ROI N T rf r T f N] ad U d F a [K(p (n 1) ) Û(n 1) ] a K ad (p (n) )U d 1 C A (2.21) and g(x)=g(x + N(x) Û(n) ).
The updating of model parameters (for fixed admissible fields) in the second step reads:

p (n+1) = argmin p2P F h mIDIC (p) with F h mIDIC (p)=E h 2 mCRE ( Û(n+1) , V(n+1) ; p) (2.22)
and is still performed by means of a steepest descent algorithm (e.g. BFGS method [Shanno, 1970]) with gradients obtained from the adjoint-state method.

As regards the setting of the weighting parameter, it is again performed from the Morozov principle; the coefficient ↵ is adjusted as the smallest positive value that enables to reach, at the end of the identification process or along iterations:

1 2 2 f N pix X x2ROI ⇣ f (x) g(x + N(x) Û) ⌘ 2 ⇡ 1 (2.23)
The overall mI-DIC algorithm is shown in Figure 2.3.

Numerical results

In this section, we perform experiments in order to check the performance of the DIC-mCRE or mI-DIC approaches.

We particularly study their robustness with respect to noise and reliability of information. All numerical examples consider 2D elasticity problems with isotropic behavior and plane stress assumption.

Example 1: simple tensile test

This first example is similar to the one considered in [Ben Azzouna et al., 2015]. It consists in a squared plate simply supported on its bottom side and subjected to a uniformly distributed traction on its top side (Figure 2.4). Other sides are free. The material is homogeneous and the objective is to identify the Young modulus E. We use here synthetic displacement measurements. Kinematic data are built numerically, by considering an analytical reference solution U 0 for a given reference set of material coefficients (E 0 , ⌫ 0 ). It is projected on a regular measurement grid of size 50 ⇥ 50, and measurement noise is generated by adding an uncorrelated Gaussian white perturbation to the reference field, under the form U i DIC = U i 0 + i ⌘ with ⌘ a zero-mean and unit variance Gaussian distributed variable. Consequently, there is no modeling error for this example. The mesh used for simulations is made of T3 elements.

Robustness with noise level

We evaluate the robustness of the DIC-mCRE approach by studying the dependency of the cost function with the noise level. For a value of the weight factor ↵ set using the Morozov principle, the values of the mCRE functional E h 2 mCRE ( Û, V; E) and its two components

E h 2 CRE = 1 2 ( Û V) T K(E)( Û V) and kΠ Û U DIC k 2 Σ Σ 1 DIC = 1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC )
are represented in Figure 2.5 as a function of E/E 0 . Various noise levels are considered, with moderate perturbations from 0% to 5% of the mean of the displacement field; the solid line corresponds to the case where no perturbation is added. For all tested noise levels, we note that the minimum of the cost function perfectly corresponds to the reference value. Besides, the component on mismatch with measurements is insensitive to the noise level in this example (the six curves are superposed); this is mainly due to the strategy used to select ↵. mCRE and its two components as a function of E/E 0 , for various noise levels (from 0% to 5%).

For more significant noise levels (20% and 50% of the mean of the displacement field), evolutions are given in Figures 2.6 and 2.7. We observe that the reached minima are a bit shifted but remain close to the reference value. mCRE and its two components as a function of E/E 0 , for a 20% noise level.

mI-DIC strategy

We now investigate the integrated version of mCRE (mI-DIC). Using the same reference solution as before with additional noise to get synthetic data, we now work at the fine grid level that would represented pixels of an image in practice. The synthetic displacement field over this grid is represented in Figure 2.8. A sub-pixel discretization is used to obtain the synthetic image g and to compute the updated deformed images g at each iteration.

Using similar simulation meshes, the same identification results as for DIC-mCRE are obtained for a range of low noise levels (up to 5%). In Figure 2.9, we display the evolution of the obtained confidence interval on the identified parameter, over a larger range of noise levels. We observe, as expected, that the accuracy on the identified parameter degrades with the noise level, but it remains acceptable even for large noise levels.

Example 2: plate with a hole

The second example considers a tensile test on a plate with a hole (Figure 2.10). Using symmetry properties, only one quarter of the plate is conserved. The simulation mesh, also shown in Figure 2.10, is made of 368 T3 elements.

The material is homogeneous and the objective is to identify the Young modulus E. Synthetic data are again used, considering a reference solution U 0 obtained from a simulation with the reference set of material coefficients E 0 = 200 GPa and ⌫ 0 =0.26. It is represented in Figure 2.11. Consequently, here again there is no modeling error.

Analysis with moderate noise level

We first analyze the performance of the mCRE-based identification method when choosing a 10% additive measurement noise. The distribution of noise and resulting measured displacements are given in Figure 2.12.

In order to analyze the effect of the weighting factor ↵, it is defined as ↵ = 10 β .U T 0 K 0 U 0 where U 0 and K 0 are displacement field and stiffness matrix associated with the reference solution. The coefficient is tuned in the range [ 2,[START_REF]1 Left: a coarse mesh of the domain Ω is introduced, with elements of diameter H much larger than the small characteristic size " of the heterogeneities (here the diameter of the inclusions). Right: on each coarse element, we solve a local problem, using in practice a discretization at the size h adapted to the heterogeneities[END_REF], and the curve associated with the Morozov principle is represented in Figure 2.13. We observe that for small , and as the model is here compatible with experimental information, the weighted term on model-data discrepancy stagnates to the desired target 1. Conversely, choosing a too large makes this term go to 0, which is undesirable in the case of noisy data. In this specific example, we select the value of that satisfies the criterion on Morozov principle (with measurement data approached up to noise level) while not over-smoothing too much the solution from the model. The curve thus indicates here to choose =0. For this optimal value =0 , the evolution of the mCRE functional and its two components is represented in Figure 2.14. We observe that the value of E is correctly identified.

We represent in Figure 2.15 the admissible fields û and v (and their difference, which is actually the field ) obtained at the end of the identification process. We observe that they differ very slightly, which is due to the absence of modeling error in this example. We also show the comparison between the kinematically admissible field and the measurement field in Figure 2.16.

Analysis with large noise level

We now consider a larger 30% additive measurement noise. The distribution of noise and resulting measured displacements are given in Figure 2.17.

When tuning in the range [ 2,[START_REF]1 Left: a coarse mesh of the domain Ω is introduced, with elements of diameter H much larger than the small characteristic size " of the heterogeneities (here the diameter of the inclusions). Right: on each coarse element, we solve a local problem, using in practice a discretization at the size h adapted to the heterogeneities[END_REF], the curve associated with the Morozov principle is represented in Figure 2.18. For the optimal value =0 , the evolution of the mCRE functional and its two components is represented in Figure 2.19. We observe that the value of E is again correctly identified.

Last, we show in Figure 2.20 the admissible fields û and v (and their difference) obtained at the end of the identification process. Again, they differ very slightly due to the absence of modeling error in this example. We also show the comparison between the kinematically admissible field and the measurement field in Figure 2.21.

Example 3: three-point bending test

We now consider an isotropic beam subjected to three-point bending (Figure 2.22), in a similar way as a test-case presented in [Ben Azzouna et al., 2015]. Boundary conditions are assumed to be perfectly known. The objective is again to identify the homogeneous Young modulus E, using the DIC-FEMU method. Displacement measurements are created from a reference FE calculation with specified boundary conditions and reference parameter E 0 (Fig-

Identification results

In Figure 2.24, we show the evolution of the mCRE cost function with respect to the ratio E/E 0 , for various perturbation levels. We again observe the good convexity and estimation of the minimum up to quite high perturbation levels (up to 20%), as well as the convexity of the cost function.

In order to show the practical interest of mCRE, we solve the identification problem by using a restricted measurement zone Ω m ⇢ Ω. Reliable information available outside the measurement zone is then included in admissibility spaces. We analyze for cases (see Figure 2.25) corresponding to ratios |Ω m |/Ω =0 .48, 0.64, 0.81, and 1, with a measurement noise of 5%. Results are reported in Figure 2.26, where we show the identified value and confidence interval. It appears that identification error increases but remains reasonable for smaller measurements zones.

Comparison with FEMU

In this section, we compare identification results obtained from the proposed mCRE method and the classical FEMU method, using similar data. We assume Ω m = Ω. Available data is the measured displacement field, as well as reliable information on the boundary @Ω. Boundary conditions for FEMU are thus made of available displacement measurements on the boundary and free edge conditions which are strongly prescribed (Figure 2.27). In order to perform sound comparison, the mCRE method only uses free-edge boundary conditions to define admissibility spaces.

Identification results are shown in Figure 2.28. We observe that mCRE performs better than FEMU. Bad results obtained with FEMU are due to perturbed boundary conditions which are prescribed to very few nodes (supported points of the beam), hence avoiding any averaging of the noise and leading to a non-representative loading in the FEMU simulation; the mCRE method better manages such a case.

Example 4: identification from real images

As a last example in this chapter, we consider the biaxial test described in [Avril et al., 2008a, Leclerc et al., 2009].

In this test, a cross-shaped specimen (see Figure 2.29) is loaded in a multiaxial testing machine, in an equibiaxial manner F x = F y (same loading magnitude in both directions). The specimen material is made of vinylester matrix reinforced by glass fibers. The quasi-uniform distribution of fiber orientations leads to a quasi-isotropic elastic behavior priori to matrix cracking and fiber breakage.

Here, we analyze experimental data obtained for the load level F x = F x =5kN, which still corresponds to the elastic regime of the composite material. The displacement field is measured by DIC, with a measurement zone covering the whole specimen geometry. Images are taken with a 8-bit CCD camera (resolution 1008 ⇥ 1016 pixels).

The measurement noise is evaluated from a first image with almost zero load. Two DIC meshes are considered (Figure 2.29): a coarser mesh with 781 nodes, and a finer mesh with 2865 nodes. These also correspond to the simulation meshes in the mCRE strategy.

Reference values E 0 = 10 GPa and ⌫ 0 =0.28 are considered for the Young modulus and Poisson ratio, respectively. They are taken from identification results in [Avril et al., 2008a]. Here, and as we only consider experimental information coming from measured displacements (with Dirichlet boundary conditions), we identify the Poisson ratio ⌫ alone. This identification is made difficult due to the fact that F x = F y .

Starting with the coarser mesh (and therefore reduced measurement uncertainty), the displacement field ob- In order to implement the Morozov principle, we again define the weighting factor ↵ as ↵ = 10 β .U T 0 K 0 U 0 where U 0 and K 0 are displacement field and stiffness matrix associated with a reference solution. The obtained Morozov curve is shown in Figure 2.31. When is too small, there is loss of information contained in the measurements and over-smoothing of the solution, resulting to values of the model-data discrepancy term larger than 1. Nevertheless, the curve does not go too high for small , as admissible fields which are compared to measurements are already constructed from both model and measurements, so that model correction from measurement data is already performed. The curve indicates that the optimal weight value is =0 . For this value, the evolution of the mCRE cost function and its two components is displayed in Figure 2.32. We observe that the identified value is ⌫ ⇡ 1.3⌫ 0 , and that the cost function has good convexity properties.

We represent in Figure 2.33 the admissible fields û and v (and their difference) obtained at the end of the identification process. The difference is not so small now, as there is modeling error. We also show the comparison between the kinematically admissible field and the measurement field in Figure 2.34.

We now perform the same analysis but with the second (finer) DIC mesh, and consequently with a larger measurement uncertainty. The associated displacement field obtained from DIC is shown in Figure 2.35 (we removed black lines of element edges for more readability).

The new Morozov curve, shown in Figure 2.36, indicates the optimal value = 0.5. For this value, the evolution of the mCRE cost function and its two components is displayed in Figure 2.37. We observe that the identified value remains ⌫ ⇡ 1.3⌫ 0 , and that the cost function still has good convexity properties.

We represent in Figure 2.38 the admissible fields û and v (and their difference) obtained at the end of the identification process. Similar trends as for the coarser mesh are observed. We also show the comparison between the kinematically admissible field and the measurement field in Figure 2.39. Chapter 3

Modeling error estimation with adaptive strategy

In this chapter, we introduce an error estimation method, and an associated adaptive strategy, in order to define models which are compatible with noisy data information coming from full-field measurements. For this purpose, we show that the analysis of the CRE term in the context of DIC-mCRE or mI-DIC methods appears as a natural tool to define appropriate models with optimal accuracy with regards to experimental information. Consequently, we highlight that the mCRE framework enables the whole structure of a reference model, not only its parameters, to be selected from available measurements. We also show that parts of the modeling error coming from the simulation mesh (discretization error) and from the mathematical model itself can be dissociated in order to drive the adaptive process in an optimal manner. The proposed strategy is here illustrated with several numerical experiments involving linear models.

Motivation and objectives

Out of model input uncertainty, that is uncertainty regarding the values of input parameters of a selected model, modeling uncertainty arises as a result of assumptions and simplifications made in answer to a lack of knowledge on the physical system. Modeling errors are caused by all hypotheses and shortcuts taken while translating physical phenomena into equations (mathematical model), and translating these equations into a form that can be solved numerically (discretization). They may have a large impact in inverse problems due to high sensitivity of the solution;

a global optimum of the inverse problem may for instance be found with unrealistic physical values for the material properties as a consequence of such modeling errors. Conversely, if the goal of modeling is not further prediction, but merely the characterization of physical properties from experimental data, a complex model (with large number of parameters) may yield large uncertainty intervals on the parameter estimates. Furthermore, the mesh size should be designed in association with the spatial accuracy of experimental observations. Adaptive modeling (in terms of model class and discretization mesh) thus appears as an important issue for effective inverse analysis; in particular, it can be seen as a suitable manner to save computational cost and perform parameter reduction [START_REF] Lieberman | Parameter and state model reduction for large-scale statistical inverse problems[END_REF], with a number of input parameters that progressively increases with local enrichment of the model.

Questioning and assessing the quality of the model, in terms of model class and discretization selection, is not a common practice in inverse problems. Quite few works addressed this issue [START_REF] Mottershead | Model updating in structural dynamics: a survey[END_REF],Becker and Vexler, 2004, Becker and Vexler, 2005, Arridge et al., 2006, Kaipio and Somersalo, 2007, Johansson et al., 2007, Oden et al., 2010, Simoen et al., 2015, Calvetti et al., 2018], with procedures being usually computationally intensive.

In Bayesian statistical inversion theory, the computational model inaccuracy may be represented as a random variable (treated as noise) with statistics that are updated from data. Moreover, information on the quality of the model can be recovered through the normalization constant C (appearing in (1.4)). This constant corresponds to the probability of making observations d obs for a given model class M; it is called model evidence and informs on the relevance of a model with respect to observed data. It can thus be used to discriminate between concurrent model classes (keeping the most likely to have generated the data), and to perform model class selection [Beck andYuen, 2004, Mthembu et al., 2011]. Considering for instance N models {M 1 ,...,M N }, each model M k depending on a parameter set p k , the pdf informing on the probability that model M k generated data d obs reads:

⇡(M k |d obs )= ⇡(d obs |M k ).⇡(M k ) P N k=1 ⇡(d obs |M k ).⇡(M k ) / ⇡(d obs |M k ).⇡(M k ) (3.1)
where ⇡(M k ) is the prior pdf (a priori knowledge) on model M k , while ⇡(d obs |M k ) is the integrated likelihood function of M k , computed by marginalization:

⇡(d obs |M k )= Z ⇡(d obs , p k |M k )dp k = Z ⇡(d obs |p k , M k ).⇡(p k |M k )dp k = C k (3.2)
where C k is the model evidence associated with M k . Therefore, when prior pdfs ⇡(M k ) (k 2 {1,...,N}) are equal, we may define the Bayes factor F ij between two models M i and M j :

F ij = ⇡(M i |d obs ) ⇡(M j |d obs ) = C i C j (3.3) It indicates when F ij > 1 (resp. F ij < 1
), that the model M i is more (resp. less) likely than the model M j to be associated with the considered data d obs .

In a deterministic framework, cross validation may be used to a posteriori assess the quality of a model [START_REF] Hastie | Elements of Statistical Learning: data mining, inference, and prediction[END_REF]. Typically, the original dataset is split in two sets (holdout method): (i) a training set used to iden-tify/update the model; (ii) a test set for its validation with error estimation. Goal-oriented estimates on the discretization error have also been designed [START_REF] Becker | A posteriori error estimation for finite element discretization of parameter identification problems[END_REF], Johansson et al., 2007, Johansson et al., 2011] for identification problems.

In the context of full-field measurements, an important feature is that data are numerous with noise that can be quantitatively characterized. This rich experimental information not only needs to be propagated throughout the identification procedure, it also needs be compared and considered with respect to other error sources (modeling, discretization) involved in the identification process. Most applications using full-field measurements assume that there is no modeling of discretization error; the first source is usually poorly known a priori, while the second source is usually neglected assuming that the mesh is fine enough. This last assumption seams sensible for localized sensors (e.g. strain gauges) where experimental richness is usually far below the kinematics offered by the mesh, so that discretization error is not an issue, but it is questionable when considering full-field measurements as discretization error may then become significant compared to the measurement error. Consequently, there is a need to adapt the numerical model (model class and mesh) used with full-field measurements, in order to obtain a numerical complexity which is consistent with the available experimental information, and therefore conduct an effective post-processing of full-field measurement data (right computation at the right cost for the objective of identification).

On the one hand, a suitable mathematical model with appropriate physics should be considered as a representation of image-based measurements in order to perform a relevant interpretation and comparison with available noisy data. Model selection with full-field measurements was addressed in [START_REF] Neggers | Improving full-field identification using progressive model enrichments[END_REF]] by using information from correlation residuals toward successive enrichment of the constitutive model in order to progressively reduce the experiment-model gap; we propose here another numerical strategy.

On the other hand, special attention should be paid to the mesh size selection (associated with discretization error) in order to exploit measurement information at best with meaningful accuracy. In [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF], mesh adaptivity was performed in the context of I-DIC by using an empirical (no error estimate) procedure that is based on sensitivity fields to define the mesh size; an example of obtained adapted mesh is given in Figure 3.1.

In [START_REF] Wittevrongel | A self adaptive global digital image correlation algorithm[END_REF], p-adaptive FE analysis was implemented in DIC to obtain a self-adapting higher order mesh capable of describing high-gradient displacement fields, here again without clear indicator on discretization error.

The objective in this chapter is to develop error measures capable of quantifying the quality of the approximate numerical model with respect to full-field measurements. These should consider simplified modeling and discretization error, and they should be related to the limited experimental information with noise. Therefore, an associated question is to define threshold values for these error measures which permit to certify that the numerical model is valid. We thus propose a framework, based on the mCRE functional, in which all a priori knowledge and uncertainty Figure 3.1: Example of adaptive mesh I-DIC procedure, from [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]. The element size in the adapted mesh (right) is a decreasing function of sensitivity k@u s /@pk.

sources (coming from both model and experimental data) are considered, informed, connected, and propagated throughout the identification procedure in order to build a consistent and unified methodology for guiding model adaptation.

The use of the CRE concept in model verification showed that the value of the CRE functional, when computed with appropriate admissible fields, is a valuable, robust, and quantitative information on the quality of the model (see Section 1.1.2). The goal is to transfer this powerful tool to image-based identification procedures. A preliminary study was conducted [START_REF] Hild | Toward 4d mechanical correlation[END_REF], where the CRE functional was a posteriori used with DIC data in order to control the quality of the mesh, for a given mathematical model and depending on the way boundary conditions are prescribed from kinematic measurements in I-DIC simulations (e.g. depending on the polynomial order used to interpolate measured fields).

We propose to take advantage of specific features of parameter identification with mCRE. Actually, the mCRE functional contains a mCRE term that can be interpreted as a distance to a given model. We show that this provides a posteriori error indicators on modeling and discretization that identify the computational model inaccuracies with regards to experimental data. These indicators are further used in a greedy manner for adaptive modeling purposes, in order to perform mesh adaptivity and select an appropriate model (inside a manifold of model classes with increasing complexity and number of parameters). Therefore, bias in the numerical model (deviations between model predictions and reality observed through full-field measurements) can be corrected and adjusted so that model complexity and outputs remain consistent with physical observations. The idea is to have a confidence in the model that is driven by confidence in measurements.

The framework proposed in this chapter is geometrically represented (Figure 3.2), as an extension to the representation of mCRE proposed in Figure 1.6. We now play with a manifold (Γ +obs M ) of possible models, with choice on the model class M (in practice among a list of parametrized constitutive models with hierarchic complexity) and on the mesh T h used to represent the state of the system.

Global estimate and specific indicators on modeling errors

For the sake of clarity, we develop the strategy proposed in this chapter by considering the DIC-mCRE version of the identification method; a similar strategy can be developed for mI-DIC. We thus consider the following mCRE functional:

E 2 mCRE (û, σ; p)=E 2 CRE (û, σ; p)+ ↵ 2 . 1 N DIC (d(û) U DIC ) T ΣΣ 1 DIC (d(û) U DIC ) with E 2 CRE (û, σ; p)= 1 2 Z Ω (σ K(p)ε(û))K(p) 1 (σ K(p)ε(û)) (3.4)
or, in its discretized version:

E h 2 mCRE ( Û, V; p)= 1 2 ( Û V) T K(p)( Û V)+ ↵ 2 . 1 N DIC (ΠU U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) (3.5)
with ↵ a scalar value that has the dimension of an energy.

A dimensionless mCRE functional could also be defined, under the form:

E 2 mCRE (û, σ; p)=E 2 CRE (û, σ; p)+ ↵ 2 . 1 N DIC (d(û) U DIC ) T ΣΣ 1 DIC (d(û) U DIC ) (3.6) with ↵ =↵.E 2 0 and E 2 mCRE = E 2 mCRE /E 2 0 . The quantity E 2 0 is a reference energy value; a typical choice is E 2 0 = 1 2 Û(0) T K(p (0) ) Û(0)
where initial values of the parameter set (p (0) ) and of admissible fields ( Û(0) ) are used. In the following, an alternative definition of the reference energy value is proposed to conduct an adaptive modeling algorithm.

Estimate on the mathematical model error

In order to define a quantitative criterion for adaptive modeling, we first analyze the idealistic scenario when the outputs of the updated model perfectly reproduce the real physical observations, that is ΠU = U true obs where U is the (admissible) solution of the direct problem K(p sol )U = F. In such a scenario, the updated model is compatible with physics so that there is no modeling error, while U DIC ΠU is a random vector ⇠ with zero mean and covariance ΣΣ DIC . We also note that in the discretized mCRE framework, the optimal admissible field V (for given p) corresponds to a kinematically admissible field of the form [U d , Va ] (as V = Û + Λ + a rigid body displacement which is considered 0) which further satisfies K(p) V = F. Consequently, for p = p sol , V = U.

In the hypothetic case of no measurement noise, then DIC measurements exactly correspond to the outputs of the updated model K(p sol )U = F, i.e. U DIC = ΠU (= Π V when p = p sol ). Therefore, at convergence, the system describing the optimal admissible fields (Step 1 of the mCRE minimization) reads:

[K aa (p sol )+↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa ] Ûa = F a + ↵. 1 N DIC (Π T ΣΣ 1 DIC ) ao U DIC K ad (p sol )U d ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) ad U d = K ao (p sol ) V + ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) ao V K ad (p sol )U d ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) ad U d =[K aa (p sol )+↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa ] Va (3.7)
so that Ûa = Va , and consequently Û = V (= U).

Therefore, E h 2 mCRE ( Û, V; p sol )=0with its two terms (including E h 2 CRE ( Û, V; p sol ) vanishing for any value of ↵.

When now considering DIC measurement noise, still with an exact model, then U DIC = Π V +⇠ (when p = p sol ).

Step 1 of mCRE minimization thus yields:

[K aa (p sol )+↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa ] Ûa = K aa (p sol ) Va + ↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa Va + ↵. 1 N DIC (Π T ΣΣ 1 DIC ) ao ⇠ (3.8) so that Ûa = Va +[K aa (p sol )+↵. 1 N DIC (Π T ΣΣ 1 DIC Π) aa ] 1 ↵. 1 N DIC (Π T ΣΣ 1 DIC ) ao ⇠ (3.9)
Consequently, the value of the CRE term is:

E h 2 CRE ( Û, V; p sol )= 1 2 ( Û V) T K(p sol )( Û V)= 1 2 ( Ûa Va ) T K aa (p sol )( Ûa Va ) = ↵ 2 2N 2 DIC .⇠ T (Π T ΣΣ 1 DIC ) T ao K T aa (p sol , ↵)K aa (p sol )K 1 aa (p sol , ↵)(Π T ΣΣ 1 DIC ) ao ⇠ (3.10)
where we again used the notation K aa (p sol , ↵)=K aa (p sol )+↵. Therefore, at the end of the identification process with p = p sol , the value of the CRE term E h 2 CRE ( Û, V; p sol ) naturally informs on a mismatch in the mathematical model:

• when the model is compatible with (noisy) observations, we should get E h CRE ( Û, V; p sol )/E ref (↵) ⇡ 1 (=1 on average); this is the configuration shown in Figure 3.3;

• when E h CRE ( Û, V; p sol )/E ref (↵)
1, this indicates bias in the employed mathematical model. This model is then too poor and deviates from a model that should represent the partial physical reality seen from observations;

• when E h CRE ( Û, V; p sol )/E ref (↵) ⌧ 1, the employed mathematical model appears to be too rich with regards to available experimental information and associated measurement uncertainty. Remark . We remind that in the strategy developed in the PhD, the value of ↵ is systematically calibrated such that the measurement error term in the mCRE functional is made consistent with measurement noise (see Section 2.1.2).

This term should satisfy (Morozov principle):

1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) ⇡ 1 (3.11)
indicating that DIC data are approached up to noise level (i.e. that the residual vector on measurements is noise only).

The value E h 2 CRE ( Û, V; p sol )= 1 2 ( Û V)K(p sol )( Û V) of the CRE term in the mCRE functional can thus be interpreted, at convergence of the identification algorithm, as a modeling error estimator ⌘ 2 mod . It informs on the quality of the updated model in view of noisy measurements, when compared with the reference value E 2 ref (↵). The estimate ⌘ mod can be further used as a criterion in an adaptive strategy, in order to select a model that suitably explains observations (see later).

Remark . The absolute estimate ⌘ mod implicitly takes into account model and measurements as admissible fields Û and V depend on them.

Error indicators and adaptive strategy

The modeling error estimate

⌘ 2 mod = 1 2 ( Û V)K(p sol )( Û V)
, derived from a discretized version of the mCRE functional, does not consider discretization error as it refers to a reference model which is already discretized in space. The stress field associated with V, that is σ( V)=Kε(N V), is indeed merely equilibrated in a FE sense. It satisfies:

V T (K(p) V F)=0 8 V 2 U h,0 ad (3.12)
The estimate ⌘ 2 mod thus informs on the error due to the mathematical model alone.

In order to capture the discretization error source, and therefore evaluate the strict quality of the mesh, one needs to refer to a continuous reference model. It is thus necessary to recover a stress field σ which is fully equilibrated (as performed in FE model verification with CRE, see Section 1.1.2). It should satisfy the property (assuming that boundary conditions on @ 2 Ω are reliable):

Z Ω σ : ε(v)= Z Ω f v d • v + Z ∂2Ω f s d • v 8v 2 U 0 ad (3.13)
which is related to static admissibility.

For this purpose, we propose to use the classical hybrid-flux (EET) technique (see Section 1.1.2) with direct postprocessing of the FE stress field σ( V) at hand after conducting the identification process with mCRE. The hybrid-flux can be employed without any change as σ( V) satisfies the required properties (in terms of FE equilibration). This a posteriori construction of a fully equilibrated admissible stress field σ thus enables one to obtain a modeling error estimate ⌘ 2 tot that integrates both discretization and mathematical model error sources. It is constructed from the CRE term, and reads:

⌘ 2 tot = E 2 CRE (û, σ; p sol ) (3.14)
From orthogonality properties, an error indicator on the discretization error alone is then recovered as:

⌘ 2 dis = ⌘ 2 tot ⌘ 2 mod (3.15)
Remark . When it is envisioned to generate synthetic measurement data using a FE mesh, and if this same mesh is further used for simulation of admissible fields in mCRE, there is no discretization error (with respect to measurements) so that the indicator ⌘ dis on discretization error should be close to 0.

Remark . In [START_REF] Geymonat | Identification of elastic parameters by displacement field measurement[END_REF], Geymonat and Pagano, 2003, Florentin and Lubineau, 2010, Moussawi et al., 2013], inverse analysis with the CRE concept is performed by using stress fields which are made statically admissible in a strong sense from the start. The strict equilibration is obtained using Airy functions in [START_REF] Geymonat | Identification of elastic parameters by displacement field measurement[END_REF], Geymonat and Pagano, 2003, Latourte et al., 2008]. In [Florentin andLubineau, 2010, Moussawi et al., 2013],

a global finite-dimensional optimization problem is solved by considering a regular simulation grid; unknowns are then equilibrated polynomial tractions over element edges, after defining a direct link between these tractions and the inside local stress field. Such procedures may yield an estimate on the global modeling error, but not indicators on individual error sources.

Adaptive process

From the above modeling error estimator ⌘ tot and indicators (⌘ mod , ⌘ dis ), and by comparison with the reference value E ref , an adaptive strategy is proposed to drive optimal model and mesh selection. We highlight that optimality has to be understood with the objective of identification; adapted model and mesh, dedicated to parameter identification and implicitly considering sensitivity analysis and experimental information, are strongly coupled with the amount and richness of measurements. They differ from adapted model and mesh that should be used for prediction in a direct problem.

Starting from an initial coarse configuration for the model class (M (0) ) and mesh (T

(0)
h ), adaptivity is performed in a greedy manner at iteration k. When ⌘

(k) tot E (k)
ref , indicating that the numerical model is not consistent with experimental information, adaptivity is conducted by comparing relative values of M (0) and mesh T (0) h , as well as their local spatial contributions to perform local adaptation. A scheme of the adaptive algorithm is given in Figure 3.4.

Remark . In practice, a change in model class is here performed over the whole domain Ω, but coupling between concurrent models could also be envisioned in order to permit local change in the model class. This would be a relevant strategy for keeping coarse model and mesh in regions which are not sensitive to identified parameters.

At the end of the adaptive process, the algorithm should tend to a simulated model that is consistent with data, yielding savings in modeling complexity without sacrificing accuracy. When writing a normalized mCRE functional (3.6), with normalization from E 2 ref , modeling and measurement terms should be of the same unit order (see Figure 3.3). Also, a balance between modeling and discretization error is intended.

Remark . We define here an adaptive process which is driven by the measurement noise level. We indicate that the reverse procedure could also be possible, that is with given reference model and mesh, we could define the richness of experimental information which is consistent to use. In the case of full-field measurements, this would come down to selecting a consistent image resolution (e.g. compressing data images with SVD) or a suitable correlation mesh (e.g. to get iso-measurement noise) within a multiscale storage of images.

Start k =0 , M (0) , T (0) h mCRE solver ⌘ (k) tot /E (k) ref 1 End ⌘ (k) mod > ⌘ (k) dis k = k +1 , M (k) , T (k)
Remark . When the same mesh is used between simulations and digital image correlation in DIC-mCRE, or when mI-DIC is used, the proposed mesh adaptation procedure can also be seen as an automatic manner to filter useless or redundant experimental information for model identification purposes. It is a pruning process, implicitly based on sensitivity analysis and measurement reliability, that enables saving in the amount of data [START_REF] Neggers | Improving full-field identification using progressive model enrichments[END_REF].

Numerical results

Application 1: heterogeneous material

As a first illustration, we consider the identification of elastic properties of a plate with 3 inclusions (with known positions) inside a matrix, from a tensile test (Figure 3.5). Each material is assumed to be isotropic, with coefficients andinclusions 1, 2, and3). We consider that ⌫ i =0 .3 (i =0 , 1, 2, 3) and that

(E i , ⌫ i ) in each phase (matrix 0
E 0 = 200 GPa.
The objective is to identify the other Young moduli E 1 , E 2 , and E 3 . tot and adapted mesh are also given in Figure 3.7. As could be foreseen, the model class is adapted during the first iterations in order to progressively incorporate the 3 inclusions which are not represented in the initial model, then the mesh is adapted around inclusions and in the vicinity of the edge with Dirichlet boundary conditions; the corresponding adapted mesh is 4770 dofs and 7712 dofs, respectively. After iteration 6, the modeling error has a level which is consistent with experimental information and the algorithm is stopped.

Application 2: anisotropic material

As a second illustration, we consider the identification of homogeneous orthotropic elastic properties. We consider again the example of Section 2.2.2 (tensile test on a plate with a hole), but now synthetic data are generated from An additive 5% noise is added.

Starting from a priori knowledge on homogeneity and on the value of the Poisson ratio alone, an isotropic elasticity model class M (0) with regular coarse mesh is first used to identified the associated Young modulus E.

The identification process indicates E ⇡ 80 GPa as well as a very large modeling error mainly coming from the choice of the mathematical model (too poor model).

Consequently, the model class is then changed into an orthotropic elasticity model, with a larger set of associated parameters (E x , E y , G xy ) to be identified. The mesh size is kept unchanged. At the end of the identification procedure, the error indicators then show that the modeling error has decreased but the discretization error source remains high (see Figure 3.8).

Mesh adaptation is thus performed using subdivision of existing elements (with management of hanging nodes)

until error indicators have converged to acceptable values. The values obtained along the adaptive process are given in Table 3.1.

Remark . Identification of anisotropic material parameters from full-field measurements was investigated in several works [START_REF] Bruno | A full-field approach for the elastic characterization of anisotropic materials[END_REF][START_REF] Grédiac | Special virtual fields for the direct determination of material parameters with the virtual fields method. 2-application to in-plane properties[END_REF],Molimard et al., 2005,Lecompte et al., 2007,Rahmani et al., 2014,Bal et al., 2014, Bal et al., 2015]. It was addressed using the mCRE concept in [START_REF] Guchhait | Anisotropic linear elastic parameter estimation using error in the constitutive equation functional[END_REF]. 3.2 1.9 2.6

Table 3.1: Relative estimate and indicators on modeling error obtained along the adaptive process.

ROM techniques, unlike meta-modeling techniques, do not simplify physics models but rather decrease their computational complexity by using specific numerical tools that generate an adequate approximate solution from a low-dimensional basis (manifold), facilitating the map from the input space to the set of outputs. They exploit the fact that the response of complex models can often be approximated with a reasonable accuracy using a reduced basis. Consequently, ROM techniques have the potential to dramatically reduce CPU costs and memory resources without sacrificing too much of the solution accuracy. They have been in rapid expansion over the last decade and their performance, in terms of savings in computational time and memory storage, are impressive (several orders of magnitude).

Most ROM procedures consist of the generation, in an intensive offline (learning) stage, of a relevant reducedorder basis that captures the dominant dynamics of the physical model. This basis is then operated in an online phase to obtain approximate solutions at low cost. We may list here:

• the POD method [Chatterjee, 2000, Kunisch and Volkwein, 2001, Kunisch and Xie, 2005, Gunzburger et al., 2007], which is similar to the Singular Value Decomposition (SVD), the Principal Component analysis (PCA),

or the Karhunen-Loeve Decomposition (KLD);

• the Reduced Basis (RB) method [START_REF] Maday | A reduced-basis element method[END_REF], Barrault et al., 2004, Maday, 2006, Rozza et al., 2008, Drohmann et al., 2012];

• the PGD method [START_REF] Chinesta | A short review on model reduction based on proper generalized decomposition[END_REF], Chinesta et al., 2014], which will be specifically detailed below.

In the case of nonlinear problems, a second reduction procedure aiming at reducing the evaluation step over a lower dimensional space is needed; this may be performed with several methods such as the Empirical Integration Method (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF], Maday and Mula, 2013, Radermacher and Reese, 2016] or the hyper-reduction method [Ryckelynck, 2009], to name a few.

ROM is effective to address multi-query procedures and parametrized problems encountered in many computational engineering activities such as optimization (sensitivity analysis), inverse analysis, uncertainty propagation, or optimal control [START_REF] Grepl | Certified rapid solution of partial differential equations for real-time parameter estimation and optimization[END_REF], Nguyen et al., 2010, Ghnatios et al., 2012, Maday et al., 2015a, Cui et al., 2015, Nadal et al., 2015, Yu and Chakravorty, 2015, Manzoni et al., 2016, Chen et al., 2017, Karcher et al., 2018].

PGD technique

In contrast to the POD or RB methods in which the reduced-order basis is extracted from pre-computed solutions of the system (learning phase), a priori methods follow a different path by progressively building an approximate representation of the solution, without assuming any prior basis or knowledge on the problem dynamics. The Proper Generalized Decomposition (PGD) [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF], Nouy, 2010, Chinesta et al., 2011, Chinesta et al., 2014] belongs to this family. PGD origins go back to the early 1980s with radial loading approximation as a key point of the LATIN method framework [Ladevèze, 1985, Ladevèze, 1989, Ladevèze, 1999]. Developed in this context for the simulation of evolution problems, PGD consists in seeking the multidimensional solution in an offline phase by using a low-rank modal approximation, i.e. a finite sum of separated-variable functions (modes). The solution to an evolution problem is, for instance, searched as the sum of products of space by time functions (Figure 4.1).

Additional model parameters (related to material properties, boundary conditions, geometry. . . ) can also be inserted as extra-coordinates of the problem. With such a modal representation, the complexity scales linearly with the number of dimensions. Then, a progressive construction of successive best rank-one approximations is performed.

PGD does not require any knowledge on the separated functions; they are computed from scratch by the iterative solver, using a variational formulation and a greedy algorithm.

The obtained PGD approximation explicitly depends on all model parameters; it constitutes a handbook of solutions that are further particularized in the online phase with cheap and fast evaluation of the modal representation for any value of the parameters. It thus permits real-time parametric analysis computations, using portable computing platforms, for optimization, inverse identification, uncertainty propagation, or optimal control purposes. During the last decade, PGD was extensively used to solve multidimensional problems and perform efficient simulation. We may cite:

• stochastic problems [Nouy, 2008];

• multiphysics problems [START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF];

• data assimilation and inverse analysis [START_REF] Ghnatios | Proper generalized decomposition based dynamic data-driven control of thermal processes[END_REF], Gonzalez et al., 2012, Louf and Champaney, 2013, Beringhier and Gigliotti, 2015, Nadal et al., 2015, Berger et al., 2016, Chamoin et al., 2016, Marchand et al., 2016, Berger et al., 2017, Signorini et al., 2017, Badias et al., 2018, Rubio et al., 2018[START_REF] Rubio | Transport map sampling with pgd model reduction for fast dynamical bayesian data assimilation[END_REF], Rubio et al., 2019a];

• problems with varying geometry [START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF], Modesto et al., 2015, Zlotnik et al., 2015, Courard et al., 2016, Chamoin and Thai, 2019, Sevilla et al., 2020];

• DVC [START_REF] Jailin | Fast 4d tensile test monitored via x-ct: Single projection based digital volume correlation dedicated to slender samples[END_REF].

A review on PGD applications can be found in [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulation[END_REF].

Also, some approaches have also been proposed to apply the PGD to nonlinear models, using Newton-type algorithms [START_REF] Chinesta | A short review on model reduction based on proper generalized decomposition[END_REF], the LATIN-PGD method [Ladevèze, 1989,Ladevèze, 1999[START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF], Vitse et al., 2014, Néron et al., 2015, Ladevèze, 2016, Vitse et al., 2019], or alternative methods [Ryckelynck, 2009].

In addition, the certification of PGD has been investigated in several works [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF], Ladevèze and Chamoin, 2011[START_REF] Moitinho De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF], Alfaro et al., 2015, Chamoin et al., 2017, Chamoin and Thai, 2019, Reis et al., 2020].

The standard approach to compute PGD modes is the so-called progressive Galerkin approach [Nouy, 2010],

which starts from a global weak formulation of the multi-dimensional problem. Consider a general linear Ddimensional problem of the form:

Lu = g, u 2 X = X 1 ⌦ X 2 ⌦ •••⌦ X D (4.1)
where L is an operator defined on the tensor space X . PGD consists in searching an approximation u m of u in a low-dimensional tensor subspace of X made of canonical format tensors of rank m:

u m = m X i=1 w 1 i ⌦ w 2 i •••⌦ w D i ,w µ i 2 X µ (4.2)
Introducing the global weak formulation of the problem:

Find u 2 X such that B(u, v)=F (v) 8v 2 X (4.3) with B(u, v)= Z Ω1 Z Ω2 ... Z Ω D b(u, v); F (v)= Z Ω1 Z Ω2 ... Z Ω D f (v) (4.4)
and assuming the rank m 1 decomposition u m 1 is known, the rank

m decomposition u m = u m 1 +w 1 ⌦w 2 •••⌦w D is searched such that: B(u m , v)=F ( v) 8 v = w 1 ⌦ w 2 •••⌦ w D + w 1 ⌦ w 2 •••⌦ w D + •••+ w 1 ⌦ w 2 •••⌦ w D (4.5)
with w µ 2 X µ . This formulation naturally leads to a nonlinear system where a set of coupled low-dimensional problems have to be solved:

B(w 1 ⌦ w 2 •••⌦ w D , w 1 ⌦ w 2 •••⌦ w D )=R m 1 ( w 1 ⌦ w 2 •••⌦ w D ) 8 w 1 2 X 1 B(w 1 ⌦ w 2 •••⌦ w D ,w 1 ⌦ w 2 •••⌦ w D )=R m 1 (w 1 ⌦ w 2 •••⌦ w D ) 8 w 2 2 X 2 . . . = . . . B(w 1 ⌦ w 2 •••⌦ w D ,w 1 ⌦ w 2 •••⌦ w D )=R m 1 (w 1 ⌦ w 2 •••⌦ w D ) 8 w D 2 X D (4.6) with R m 1 (v)=F (v) B(u m 1 ,v).
This system is in practice solved with an iterative fixed-point (or alternated directions) strategy. Additional ingredients may be added in the modal construction in order to optimize numerical performance, such as the orthogonalization of PGD modes.

Minimization of mCRE using PGD

Presentation

In the iterative inversion technique based on the mCRE functional (see Chapter 2), Step 1 of the procedure requires to solve a potentially large linear system for many values of the parameter set p (sequential minimization). This enables to recover optimal admissible fields ( Û(p), V(p)) at each iteration (i.e. each time p is updated), to derive the cost function F mCRE (p) (or F mIDIC (p)) to be minimized, and to compute its gradient. Also, we detailed in Chapter 2 how the scaling factor ↵ could be tuned in order to satisfy the Morozov principle; this again requires to solve the linear system for several values of ↵. In this multi-query context, PGD appears as a natural tool to decrease the computational cost and accelerate the inversion method. The idea is to compute, in an offline phase, a multiparametrized admissible solution that explicitly depends on p and ↵. This solution can then be easily evaluated in the online stage when conducting the inversion process.

Pioneering works on the use of PGD into the mCRE context can be found in [START_REF] Bouclier | Real-time validation of mechanical models coupling pgd and constitutive relation error[END_REF], Chamoin et al., 2016, Marchand et al., 2016]. An important feature is that PGD modes are computed from both model and experimental data, so that they integrate sensing information; this is in opposition to approaches where the PGD solution is constructed from the direct problem [START_REF] Gonzalez | Proper generalized decomposition based dynamic data-driven inverse identification[END_REF] before being employed for inversion. We also indicate that reduced basis approaches were also used in [START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF], Oliveira et al., 2020] when minimizing the mCRE functional.

We explain the PGD strategy by considering the functional associated with the DIC-mCRE method (Section 2.1.2), that reads:

E h 2 mCRE ( Û, V; p)= 1 2 ( Û V) T K(p)( Û V)+ ↵ 2 . 1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) (4.7)
A similar strategy can be implemented when considering the functional associated with the mI-DIC method (Section 2.1.3).

Proposed strategy

The constrained minimization performed in Step 1 of the mCRE minimization is based on the following Lagrangian functional (written in its FEM discretized form):

L h ( Û, V, Λ; p)= 1 2 ( Û V) T K(p)( Û V) + ↵ 2 . 1 N DIC (Π Û U DIC ) T ΣΣ 1 DIC (Π Û U DIC ) Λ T (K(p) V F) (4.8)
The first-order Karush-Kuhn-Tucker necessary optimality conditions, by searching the saddle-point of L h , read:

ÛT ✓ K(p)( Û V)+↵. 1 N DIC Π T ΣΣ 1 DIC (Π Û U DIC ) ◆ =0 8 Û =[0, Ûa ] T VT ⇣ K(p)( V Û) K(p)Λ ⌘ =0 8 V Λ T (K(p) V F)=0 8 Λ =[0, Λ a ] T (4.9)
or, after substituting V (= Û + Λ):

ÛT ✓ K(p)Λ + ↵. 1 N DIC Π T ΣΣ 1 DIC (Π Û U DIC ) ◆ =0 8 Û =[0, Ûa ] T Λ T ⇣ K(p)( Û + Λ) F ⌘ =0 8 Λ =[0, Λ a ] T (4.10)
This Galerkin formulation (in space) may be written in a more condensed form:

b([ Û, Λ], [ Û, Λ]) = f ([ Û, Λ]) 8[ Û, Λ] (4.11) with b([ Û, Λ], [ Û, Λ]) = ÛT K(p)Λ + ↵. 1 N DIC ÛT Π T ΣΣ 1 DIC Π Û + Λ T K(p)( Û + Λ) f ([ Û, Λ]) = ↵. 1 N DIC ÛT Π T ΣΣ 1 DIC U DIC + Λ T F (4.12)
We implement a PGD reduced model by finding, in an offline phase, parametrized solutions ( Û, Λ) with ↵ and p as extra-parameters. They are searched under the form:

Ûm (↵, p)= m X i=1 2 4 Ψ U i  U i (↵) P Y j=1 U j,i (p i ) 3 5 ; Λ m (↵, p)= m X i=1 2 4 Ψ Λ i  Λ i (↵) P Y j=1 Λ j,i (p i ) 3 5 (4.13)
Remark . We assume that measurement values in U DIC are known upstream to the updating procedure, and that this procedure is conducted for a single set of measurement values. In contrary cases such as data assimilation on time-dependent problems, measurement values could also be considered as extra-parameters in the PGD decomposition [START_REF] Bouclier | Real-time validation of mechanical models coupling pgd and constitutive relation error[END_REF], Marchand et al., 2016]. To limit the number of extra-parameters, full-field measurements should then be decomposed over a reduced basis (using POD for instance).

The conventional progressive Galerkin approach described before is used with the global bilinear form B and linear form F , defined as:

B([ Û, Λ], [ Û, Λ]) = Z Ωα Z P b([ Û, Λ], [ Û, Λ]) ; F ([ Û, Λ]) = Z Ωα Z P f ([ Û, Λ]) (4.14)
with Ω α = R ⇤+ and P the spaces of parameters ↵ and p, respectively.

In the fixed point algorithm used to compute each mode, functions of ↵ and p are normalized. Therefore, the magnitude of mode i is in the space functions Ψ U i and Ψ Λ i .

The PGD solutions Ûm (↵, p) and Vm (↵, p)= Ûm (↵, p)+Λ m (↵, p) are then used in the online model updating phase with mCRE. This has several advantages:

• the explicit dependency on parameters p enables: (i) to evaluate very fast and for any value of p the optimal admissible fields arising from the first constrained minimization; (ii) to compute gradients of the cost function

F mCRE (p)
in a straightforward manner. The second minimization (correction) step of the mCRE approach is thus performed very easily, by means of simple evaluations;

• the explicit dependency on parameter ↵ makes the definition of its optimal value (with respect to measurement noise using the Morozov principle) straightforward as well.

Implementation details

For illustration, we consider only one material parameter p, that is p = {p}.

Assuming the order m 1 PGD decomposition ( Ûm 1 , Λ m 1 ) is known, the order m decomposition of the form

Ûm (↵,p)= Ûm 1 (↵,p)+Ψ U  U (↵) U (p) and Λ m (↵,p)=Λ m 1 (↵,p)+Ψ Λ  Λ (↵) Λ (p) is searched by solving: B([ Ûm , Λ m ], [ Û, Λ]) = F ([ Û, Λ]) 8[ Û, Λ] (4.15)
Test functions are chosen in the tangent space and read:

Û = Ψ U . U . U + Ψ U .  U . U + Ψ U . U . U Λ = Ψ Λ . Λ . Λ + Ψ Λ .  Λ . Λ + Ψ Λ . Λ . Λ (4.16)
Consequently, we get the following system of coupled equations:

• problem in the physical space Ω (P1):

([ U ,  Λ ], [ U , Λ ]) 7 ! [Ψ U , Ψ Λ ]=g 1 ([ U ,  Λ ], [ U , Λ ]) B([Ψ U . U . U , Ψ Λ . Λ . Λ ], [ Ψ U . U . U , Ψ Λ . Λ . Λ ]) = R m 1 ([ Ψ U . U . U , Ψ Λ . Λ . Λ ]) 8 Ψ U , Ψ Λ (4.17)
• problem in the scalar weight space Ω α (P2):

([Ψ U , Ψ Λ ], [ U , Λ ]) 7 ! [ U ,  Λ ]=g 2 ([Ψ U , Ψ Λ ], [ U , Λ ]) B([Ψ U . U . U , Ψ Λ . Λ . Λ ], [Ψ U .  U . U , Ψ Λ .  Λ . Λ ]) = R m 1 ([Ψ U .  U . U , Ψ Λ .  Λ . Λ ]) 8  U ,  Λ (4.18)
• problem in the material parameter space P (P3):

([Ψ U , Ψ Λ ], [ U ,  Λ ]) 7 ! [ U , Λ ]=g 3 ([Ψ U , Ψ Λ ], [ U ,  Λ ]) B([Ψ U . U . U , Ψ Λ . Λ . Λ ], [Ψ U . U . U , Ψ Λ . Λ . Λ ]) = R m 1 ([Ψ U . U . U , Ψ Λ . Λ . Λ ]) 8 U , Λ (4.19) with R m 1 (•)=B([ Ûm 1 , Λ m 1 ], •) F (•).
Remark . In order to take into account the possibly non-homogeneous Dirichlet boundary conditions associated with Û, an initial mode Û0 = Ψ U 0 satisfying these conditions is inserted in the PGD decomposition Ûm . The next modes in this PGD decomposition, as well as those in the PGD decomposition Λ m , then contain space functions that satisfy homogeneous Dirichlet boundary conditions.

The coupled system is solved with a fixed-point strategy; after initializing with ( U (0) ,  Λ (0) , U (0) , Λ (0) ) (constant functions in practice), the following computations are iteratively performed:

1. Solve problem P1 to compute space functions [Ψ U (k+1) , Ψ Λ (k+1) ] from [ U (k) ,  Λ (k) ] and [ U (k) , Λ (k) ] 2. Solve problem P2 to compute scalar weight functions [ U (k+1) ,  Λ (k+1) ] from [Ψ U (k+1) , Ψ Λ (k+1) ] and [ U (k) , Λ (k) ],
and normalize them 3. Solve problem P3 to compute material parameter functions k+1) ] and k+1) ,  Λ (k+1) ], and normalize them

[ U (k+1) , Λ (k+1) ] from [Ψ U (k+1) , Ψ Λ (
[ U (
The (discretized) space problem P1 leads to a linear system, while problems P2 and P3 correspond to algebraic equations, for which a discretization on a fine grid may be performed.

The iterative process to compute mode m may be stopped when it has converged (with stagnation of modal functions along iterations, or criterion on the global residual in the tangent space). Alternatively, the process may be stopped after a fixed number k max of iterations. In the numerical studies, we implemented this second choice. FE sense) in order to derive a stress field σ which is fully equilibrated. As indicated in Chapter 1, the hybrid flux technique is made of two steps. In the first one, some polynomial tractions fK (defined as linear combinations fK|Γ (x)= P j2JΓ f j K|Γ j (x) on FE shape functions i ) are constructed on edges of each element K; they should be balanced with the external loading (f v d , f s d ). In a second step, the stress field σ is locally computed on each element K of the mesh. It should satisfy:

r • σ + f v d = 0 in K ;σn = fK on @K (4.20)
It can be shown that the best stress field σ (i.e. locally minimizing the complementary energy or, equivalently, the CRE functional) is derived from a displacement field, which is σ = K(p)ε(⇢). Consequently, the independent problems that are solved in each element K read:

Z K K(p)ε(⇢):ε(v)= Z K f v d • v + Z ∂K fK • v 8v 2 U (K) (4.21) with U (K)=[H 1 (K)] d .
The solution ⇢ of these problems, which is defined up to rigid body motions in U R (K), is in practice numerically solved with a single FE element with higher-order shape functions (i.e. degree p + k with p the polynomial degree of the initial FE computation and k the extra-degrees). Numerical studies [START_REF] Babuska | Validation of a posteriori error estimators by numerical approach[END_REF] showed that accurate error estimates are obtained choosing k 3, even though the stress field is not rigorously balanced in each element

K.

The previous procedure again leads to a multi-query context, as the above problem (4.21) (illustrated in Figure 4.4) needs to be solved for various p, for various elements (i.e. computational domain geometries), and for various boundary loadings (tractions). Even though it involves local computations, it is in practice the costliest part in the hybrid-flux method (in particular for 3D applications). The PGD technique can thus be advantageously used to compute, in an offline phase and at the element level, an approximation of ⇢ which is valid for any p, any geometry of K, and any boundary loading (using superposition). This is performed after parametrization of local problems, and it leads to locally equilibrated fields (at the element level), valid for any configuration, that can further be used in the online error estimation phase. Pioneering works on this topic can be found in [START_REF] Chamoin | Synergies between the constitutive relation error concept and pgd model reduction for simplified v&v procedures[END_REF], Allier et al., 2018, Chamoin and Thai, 2019] for scalar linear elliptic (thermal) problems; they are here extended to elasticity with equilibrium on stress fields.

Proposed strategy

We detail here the PGD procedure in 2D cases, considering 3-node triangle elements; it can be extended to other cases in a straightforward manner. On each edge Γ of any element K, tractions are linear combinations of FE shape functions j and thus read:

fK|Γ (x)= X j2JΓ [ f jx K|Γ x j (x)+ f jy K|Γ y j (x)] with ( f jx K|Γ , f jy K|Γ ) T 2 R 2 (4.22)
where we used the notation x j (x)=( j (x), 0) T and y j (x)=( 0 , j (x)) T . J Γ corresponds to the set of two nodes on Γ. Consequently, the solution ⇢ to (4.21) is written from linearity as:

⇢(x)= X Γ⇢∂K X j2JΓ [ f jx K|Γ ⇢ jx Γ (x)+ f jy K|Γ ⇢ jy Γ (x)] (4.23)
where ⇢ jx/y Γ is the solution, up to rigid body motions, to the following elementary problem:

Z K K(p)ε(⇢ jx/y Γ ):ε(v)= Z Γ x/y j • v Z K (a 1 ^X + a 2 ) • v 8v 2 U (K) (4.24)
X denotes the barycentric coordinates in element K, and a 1 and a 2 are defined as:

a 1 = R Γ X ^ x/y j R K X • X ; a 2 = 1 |K| Z Γ x/y j (4.25)
They ensure equilibrium at the element level, and therefore existence of a solution to elementary problems (4.24).

Remark . The number of elementary problems (4.24) depends on the space dimension and the type of element.

For 3-node triangle elements, there are 12 elementary problems (2 for each of the 2 nodes on each of the 3 edges).

There would be 18 elementary problems for 6-node triangle elements (2 for each of the 3 nodes on each of the 3 edges), and 36 elementary problems for 4-node tetrahedron elements (3 for each of the 3 nodes on each of the 4 faces).

The solution ⇢ jx/y Γ to each problem (4.24) can be computed with the PGD technique, for any element K,b y parametrizing the geometry of K with a set of parameters p geo 2 P geo . Following the approach described in [START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF], Zlotnik et al., 2015], problem (4.24) is reformulated by introducing a parameter-dependent mapping

M (p geo ): K ref ! K(p geo ) from a reference fixed element K ref to the geometrically parametrized element K(p geo ).
Such a geometrical transformation then allows the problem to be recast in a tensor product space and PGD to be applied, in order to compute generic parametrized solutions ⇢ jx/y Γ (x ref , p geo , p), which can be used for any element geometry and material properties.

In the present case, the mapping is defined from two transformations:

• a first scaling mapping M 1 : K ! K maps a homothetic element K with diameter 1 to the actual element K with diameter h. This mapping involves a transformation matrix T 1 = hI d such that x = T 1 x;

• a second mapping M 2 : K ref ! K maps a reference element K ref to element K. This mapping uses an isoparametric formulation and involves a transformation matrix T 2 such that x = T 2 x ref .

The global mapping is therefore M (p geo )=M 1 M 2 with transformation matrix T(p geo )=hT 2 , Jacobian matrix J(p geo )=T(p geo ), and Jacobian J(p geo )=det(J(p geo )).

Considering 3-node triangle elements (Figure 4.5), mapping M 2 reads: Remark . The number of parameters involved in mapping M (i.e. the size of p geo ) again depends on the space dimension and the type of element. There are 3 geometrical parameters (6 degrees of freedom with 3 rigid body motions) for 3-node triangle elements, 9 geometrical parameters (12 degrees of freedom with 3 rigid body motions) for 6-node triangle elements, and 6 geometrical parameters (12 degrees of freedom with 6 rigid body motions) for 4-node tetrahedron elements. In this last case of a 3D problem with 4-node tetrahedron elements, the geometry parameter vector would read p geo =( h, x 3 , y 3 , x 4 , y 4 , z 4 ) (6 scalar parameters) where local coordinates of nodes 3 and 4 in the coordinate system associated with element K are involved.

0 B @ x y 1 C A = 0 B @ P 3 i=1 x i i (⌘, ⇠) P 3 i=1 y i i (⌘, ⇠) 1 C A = T 2 0 B @ ⌘ ⇠ 1 C A with T 2 =
Consequently, problem (4.24) is parametrized and defined in a fixed reference element K ref under the form:

Z K ref JK(p)[J T rr ] s ⇢ jx/y Γ :[J T rr ] s v = Z Γ ref J s x/y j • v Z K ref J(a 1 ^X + a 2 ) • v 8v 2 U (K ref ) (4.28)
with J s the Jacobian associated with the transformation of element edge Γ:

J s = h 8 > > > > < > > > > : 1 on Γ 12 p (1 x 3 ) 2 + y 2 3 on Γ 23 p x 2 3 + y 2 3 on Γ 13 (4.29)
A separated-variable approximation of J s (using HO-SVD) is defined in practice for computational purposes.

We now introduce a reduced model, based on PGD, to obtain parametrized approximations of the solutions to (4.28). We illustrate the procedure for the approximation of ⇢ jx Γ corresponding to the boundary loading x j (x)= ( j (x), 0) T (associated with node j) along the edge Γ. This approximation with variable-separated modal decomposition reads:

⇢ jx Γ,m (x ref , p geo , p)= m X i=1 0 @ i (x ref ) Y 1kPgeo  k,i (p geo,k ) Y 1lP l,i (p l ) 1 A (4.30)
It is computed in an offline phase using the previously described progressive Galerkin approach and introducing the global forms:

B(⇢ jx Γ , v)= Z Pgeo Z P Z K ref JK(p)[J T rr ] s ⇢ jx/y Γ :[J T rr ] s v F (v)= Z Pgeo Z P " Z Γ ref J s x/y j • v Z K ref J(a 1 ^X + a 2 ) • v # (4.31)
Remark . For extremal cases (elements K with bad shape), the number of PGD modes should be chosen carefully in order to ensure accuracy. However, too much distorted elements are usually avoided in practical meshes, and a few number of PGD modes leads to a PGD error which is negligible compared to the discretization error that is assessed.

Combining with (4.23), a PGD representation of the solution to (4.21) is finally obtained as: which can be directly used online in the a posteriori error estimation procedure of Chapter 3.

⇢ m (x ref , p geo , p, { f jx/y K|Γ })= X Γ⇢∂K X j2JΓ h f jx K|Γ ⇢ jx Γ,
A flowchart of the use of PGD in Step 2 of the hybrid-flux approach is given in Figure 4.6.

Implementation details

For illustration purposes, we consider the PGD decomposition of the solution ⇢ jx Γ . We also assume that there is only one geometry parameter p geo (that is p geo = {p geo }) and one material parameter p (that is p = {p}). Consequently, the decomposition reads:

⇢ jx Γ,m (x ref ,p geo ,p)= m X i=1 i (x ref ) i (p geo ) i (p) (4.34)
B(⇢ jx Γ,m , ⇢ jx Γ )=F ( ⇢ jx Γ ) 8 ⇢ jx Γ (4.35)
Test functions ⇢ jx Γ are chosen in the tangent space and read:

⇢ jx Γ = .. + . . + .. (4.36)
Consequently, we get the following system of coupled equations:

• problem in the physical space K ref (P1):

(, ) 7 ! = g 1 (, )

B( .. , .. )=R m 1 ( .. ) 8 (4.37)

• problem in the geometry parameter space P geo (P2):

( , ) 7 !  = g 2 ( , ) B( .. , . . )=R m 1 ( . . ) 8  (4.38)

• problem in the material parameter space P (P3):

( , ) 7 ! = g 3 ( , ) B( .. , .. )=R m 1 ( .. ) 8 (4.39) with R m 1 (•)=B(⇢ jx Γ,m 1 , •) F (•).
The coupled system is solved with a fixed-point strategy; after initializing with ( (0) , (0) ) (constant functions in practice), the following computations are iteratively performed:

1. Solve problem P1 to compute the space function (k+1) from  (k) and (k)

2. Solve problem P2 to compute the geometry parameter function  (k+1) from (k+1) and (k) , and normalize it 3. Solve problem P3 to compute the material parameter function (k+1) from (k+1) and  (k+1) , and normalize it

The space problem P1 is solved using FEM with a single element of degree p + k. Problems P2 and P3 correspond to algebraic equations, for which a discretization on a fine grid may be performed. The iterative process is stopped when it has converged or after a fixed number of iterations.

Numerical results

We illustrate the approach by considering again the test-case of Section 3.4.1 with heterogeneous material. When implementing PGD for this 2D problem to obtain equilibrated stress fields at the T3-element level, decompositions of the form:

⇢ jx/y Γ,m (x ref ,h,x 3 , y 3 ,E)= m X i=1 i (x ref ) h,i (h) x,i (x 3 ) y,i (y 3 ) i (E) (4.40)
are computed (E is the local value of the heterogeneous Young modulus to be identified).

We show in Figure 4.7 the first 3 modes (functions i ,  x,i and  y,i ) corresponding to the decomposition ⇢ 1x Γ12,m , associated with the elementary boundary loading x 1 (x)=( 1 (x), 0) T (thus oriented in the x-direction and evolving as the shape function of node 1) along the edge Γ 12 . Evaluations of the PGD solution in the element K, for various values of x 3 and y 3 , are displayed in Figure 4.8.

When implementing the modeling adaptation process with the PGD strategy, an average speed-up of 16 is obtained in the online phase. This speed-up level takes into account the CPU time of the offline PGD simulations. Chapter 5

Extension to nonlinear behavior with model selection

In this chapter, we extend to nonlinear constitutive (elasto-plasticity) models the mCRE framework in association with full-field measurements. We first generalize the mCRE functional by using sound thermodynamics concepts, before introducing a specific solver in order to solve the optimization problem that drives the minimization of this functional. Model selection and adaptivity are again illustrated in this context.

Modified CRE functional in the nonlinear framework

Motivation

There is a constant development of new material models with increasing constitutive complexity. As they need to be identified, considering complex nonlinear behaviors has now become usual when dealing with inverse problems.

We refer to [START_REF] Constantinescu | On the identification of elastoviscoplastic constitutive laws from indentation tests[END_REF], Corigliano and Mariani, 2001, Claire et al., 2004, Azam, 2014] for general works on this topic, and to [START_REF] Latourte | Elastoplastic behavior identification for heterogeneous loadings and materials[END_REF],Ben Azzouna et al., 2011,Blaysat et al., 2012,Bouterf et al., 2015,Mathieu et al., 2015,Barbarella et al., 2016,Neggers et al., 2017] for applications with full-field measurements.

From an experimental point of view, parameter identification of such complex nonlinear models appears to be timeconsuming and expensive. It is thus important to master the identification procedure in terms of management of error sources; in particular, selecting the correct constitutive model among a large manifold with hierarchical complexity (see [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]) is a fundamental issue. In this context, we wish here to extend the DIC-mCRE (or mI-DIC) framework that we developed in previous chapters.

Actually, the CRE concept is strongly based on constitutive models, with a sound background on continuum thermodynamics. Even though originally defined for linear constitutive models, a generalization to complex non-linear constitutive models (viscoplasticity, damage) was proposed [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF], Ladevèze et al., 1999,Ladevèze, 2001] still in the context of model verification. It was derived from a standard formulation of the material behavior involving internal variables [START_REF] Halphen | On generalized standard materials (in french)[END_REF], and using Legendre-conjugate dual convex potentials associated with constitutive laws within this formulation.

We propose to use this generalized CRE measure into the mCRE functional, thus performing a natural transition between CRE and mCRE for nonlinear models, and providing for a unified formulation of mCRE-based inverse analysis. Preliminary works on the use of nonlinear models with mCRE can be found in [START_REF] Marchand | Parameter identification and model updating in the context of nonlinear mechanical behaviors using a unified formulation of the modified constitutive relation error concept[END_REF],

and in [Hadj-Sassi, 2007, Blaysat et al., 2012] to a lesser extent. We also mention that other works proposed methodologies based on CRE-type identification tools for nonlinear behaviors [Allix andVidal, 2002, Barbarella et al., 2016], but these are specific and functionals which are introduced are tuned to the considered problem (e.g.

using an incremental version of the linear elastic case); they do not consider a pure definition of the CRE contrary to what we propose below.

In this section, we develop and analyze the methodology by considering (elasto-)plasticity models, and using the local state method (which postulates that the thermomechanical state of the medium, at given point and time, is fully defined by knowing local values of a given number of variables).

Thermodynamical framework

We consider the quasi-static evolution of an open bounded body Ω over the time interval [0,T]. Assuming small deformations and isothermal conditions, the (nonlinear) constitutive relation reads:

σ |t = A(ε |τ , ⌧  t) 8(x,t) 2 Ω ⇥ [0,T] (5.1)
where A is a given operator. This is the functional description of the behavior, which will not be kept in the following.

We rather focus on the wide class of nonlinear material behaviors described with standard formulation [START_REF] Halphen | [END_REF]Nguyen, 1975,Germain et al., 1983]. In such a formulation, the behavior is described in terms of two complementary parts (state equations and evolution laws) that are defined by means of convex and dual (in the Legendre-Fenchel sense) potentials, and use internal variables. The associated thermodynamics framework is recalled here.

The first principe of thermodynamics, which describes conservation of energy, globally reads for any continuous medium ! ⇢ Ω:

d dt (E int + E kin )=P ext + Q or dE int dt = P int + Q (5.2)
where we introduced:

• the internal energy E int , with specific density e such that E int = R ω ⇢e;

• the kinetic energy E kin = 1 2 R ω ⇢v 2 , with v the velocity;

• the received external power

P ext = R ω f v • v + R ∂ω σn • v;
• the received heat rate Q, coming from volume heat source r and heat flux q such that Q = R ω r R ∂ω q • n;

• the internal power P int = R ω σ : rr s v.

This principle thus locally reads:

⇢ ė = σ :ε + r r • q (5.
3)

The second principle of thermodynamics is an evolution principle that involves entropy S of the material system; it reads:

dS dt Z ω r T Z ∂ω q • n T (5.4)
Introducing the specific entropy density s, this principle locally reads:

⇢ ṡ + r • ( q T ) r T 0 (5.5)
Combining the two previous local equations, we obtain the fundamental inequality:

σ :ε + ⇢(T ṡ ė) q • rT T 0 (5.6)
Last, introducing the Helmholtz free (stored) energy = ⇢(e Ts) as a new thermodynamic quantity, the following Clausius-Duhem inequality becomes:

σ :ε ˙ q • rT T 0 ( Ṫ =0here) (5.7)
It is at the heart of the consistent definition of constitutive laws, as thermodynamically admissible processes should satisfy this inequality at each time of the material evolution.

Observable variables in the proposed formalism are temperature T and total strain ε. When considering dissipative phenomena, the material state also depends on past history which is represented by internal variables. For plasticity problems of interest, elastic strain ε e and plastic strain ε p variables are introduced; they are linked with the partition ε = ε e + ε p (so that in practice only the plastic strain ε p keeps the status of internal variable). Other phenomena (such as hardening) are represented by additional internal variables X k (k =1, 2,...,K) that can be scalar, vector, or tensor quantities. They are gathered in vector X, so that the full list of state variables is (T,ε, ε p , X).W e may also use the notation X 0 = ε p . In the following, we do not include the temperature variable T in the notations anymore, due to the assumption on isothermal conditions.

State equations are derived from the Helmholtz free energy used as a thermodynamical potential (ε, ε p , X); it reads (ε ε p , X)= (ε e , X) for elasto-(visco-)plasticity (i.e. strain only appears under its elastic part). It is assumed to be convex (and concave with respect to T in a more general case) in order to get a sufficient condition to satisfy stability conditions of the second principle. The Clausium-Duhem inequality then reads:

(σ @ @ε e ):ε e + σ :ε p K X k=1 @ @X k . Ẋk q • rT T 0 (5.8)
Considering a reversible elastic transformation with uniform temperature, we directly get σ = ∂ψ ∂εe which corresponds to the elastic Hooke law. By analogy, we define variables Y k = ∂ψ ∂X k which are thermodynamic forces associated with internal variables X k in a duality pairing (k =1 , 2,...,K); they are gathered in a vector Y. The previous relations between (ε e , σ) and (X, Y) are the so-called state equations.

Remark . The free energy is usually written as the sum of an elastic contribution e and a plastic contribution p :

(ε e , X)= e (ε e )+ p (X)

(5.9) Also, it is shown in [Ladevèze, 1999] that for a large class of material behaviors, sets X and Y can be defined (using a change of variables if required) such that state equations are linear i.e. the Helmholtz free energy is quadratic (normal formulation).

Remark . It may be useful to work with a potential written in terms of stress rather than strain; this is the Gibbs free energy G (σ, X), obtained from the Helmholtz free energy by using a partial Legendre-Fenchel transform on elastic strain (with given internal variables):

G (σ, X)=sup εe [σ : ε e (ε e , X)] (5.10) State equations then read ε e = ∂ψ G ∂σ and Y k = ∂ψ G ∂X k .
In order to describe dissipative processes, in particular the evolution of internal variables, a complementary formalism is required. It is derived from the Clausius-Duhem inequality which reduces to (using state equations): (5.11) with Φ 1 the intrinsic (or mechanical) dissipation, representing plastic power amputated by the rate of energy stored in the material, and Φ 2 the thermal dissipation.

⇣ σ :ε p Y • Ẋ⌘ q • rT T = Φ 1 + Φ 2 0
The definition of evolution laws may be performed by postulating the existence of a dissipation potential '(ε p , Ẋ) expressed as a scalar convex function, non-negative and zero at the origin, such that σ = ∂ϕ ∂ εp and Y k = ∂ϕ ∂ Ẋk (normality rule). This is a sufficient condition to ensure Φ 1 '(ε p , Ẋ) 0.

To summarize, defining global flux variables e e =[ε e , X] T , e p =[ε p , X] T (with e =[ε, 0] T = e e + e p ), and global thermodynamics forces s =[σ, Y] T , the constitutive behavior is described by:

• a set of state equations: s = @ @e e or e e = @ ⇤ @s (= Λs) (5.12)

• a set of evolution laws:

s = @' @ ėp or ėp = @' ⇤ @s (= B(s)) (5.13)
where Λ is a linear, symmetric and positive operator, while B is a positive but possibly nonlinear and multivalued operator (as in plasticity). Using these notations, the intrinsic dissipation reads Φ 1 = s • ėp . Dual potentials ⇤ and ' ⇤ are defined from the Legendre-Fenchel transform:

⇤ (s)=sup ee [s • e e (e e )] ; ' ⇤ (s)=sup ėp [s • ėp '( ėp )] (5.14)
The problem of material modeling is in the determination of analytical expressions for the state potential (or ⇤ ) and the dissipation potential ' (or ' ⇤ ), and their identification from experiments (which is the topic of this chapter).

Remark . The first relation in evolution laws (5.13) drives the (visco-) plasticity flow εp while other relations drive the evolution of internal variables such as hardening.

Remark . When ' ⇤ is not differentiable at some points (usual case in elasto-plasticity with behavior independent of velocity), the evolution laws should be replaced by ėp 2 @ s ' ⇤ where @ s ' ⇤ denotes the sub-differential of ' ⇤ defined as:

@ s ' ⇤ = { ėp such that ' ⇤ (s) ' ⇤ (s) ėp • (s s) 8s} (5.15)
Introducing the convex yield function f (s)  0 associated with the indicatrix function ' ⇤ (that is ' ⇤ =0if f<0 and

' ⇤ =+1 if f =0)
, one gets: ėp = ˙ @f @s with ˙ 0 and ˙ f =0(consistency condition) (5.16) This defines associated models as surface f (defining the elasticity domain) also corresponds to the flow potential.

It is equivalent to the Hill principle of maximal work indicating that the rate ėp maximizes the intrinsic dissipation Φ 1 = s • ėp ; in this case, ˙ corresponds to a multiplier in Kuhn-Tucker conditions with constraint f  0.

Definition of the mCRE functional

A suitable CRE measure defined by Legendre-Fenchel residuals and sharing similarities with the so-called symmetrized Bregman divergence [START_REF] Chen | Metrics defined by bregman divergences[END_REF] can be defined from the two potential pairs defined in the previous section [Ladevèze, 1998, Ladevèze et al., 1999, Ladevèze, 2001, Ladevèze and Pelle, 2005]. It reads in its local in space and time version:

e 2 CRE|x,t (ê e , êp , ŝ)=⌘ ψ (ê e , ŝ)+ Z t 0 ⌘ ϕ ( ėp , ŝ) 0 8(x,t) 2 Ω ⇥ [0,T] (5.17) with ⌘ ψ (ê e , ŝ)= (ê e )+ ⇤ (ŝ) ŝ • êe 0; ⌘ ϕ ( ėp , ŝ)='( ėp )+' ⇤ (ŝ) ŝ • ėp 0 (5.18)
the residuals on state equations and evolution laws, respectively. Terms y • x correspond to the duality product between variables x and y.

Residuals ⌘ ψ and ⌘ ϕ are local in space and time quantities. They are associated with an admissible solution (ê e , êp , ŝ) such that êe + êp = ê. When they vanish, constitutive laws are locally satisfied [Moreau, 1966]. A global CRE measure E 2 CRE (ê e , êp , ŝ) is obtained by integration of e 2 CRE|x,t (ê e , êp , ŝ) over the space-time domain (see (1.24)). It vanishes when the constitutive relations are verified at any space-time point.

Remark . An alternative CRE measure, introduced in [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF], Ladevèze, 1998, Ladevèze and Moës, 1999] for visco-elastoplastic materials and referred to as dissipation error, is based on the residual associated with dissipation phenomena (i.e. evolution laws) alone. In this context, state equations are inserted in admissibility conditions. This framework was widely used for model verification purposes [START_REF] Pelle | An efficient adaptive strategy to master the global quality of viscoplastic analysis[END_REF], Chamoin and Ladevèze, 2007, Chamoin and Ladevèze, 2008, Ladevèze, 2008, Ladevèze et al., 2012]. It may be used in the context of model identification or updating when elastic properties of the material are well-known or properly identified in a first step; this will not be considered here. Also, a more general framework for the definition of dissipation error was proposed [Ladevèze, 2008] when evolution laws are not given by potentials (but provided the operator B still remains monotonic).

From the definition (5.17) of the CRE measure, a direct and natural extension of the mCRE functional to standard materials is proposed. It reads for the DIC-mCRE method: 5.19) where N t denotes the number of data assimilation time points. As it is naturally derived from the thermodynamics framework, this functional keeps advantages associated with convexity properties. It is applied to an admissible solution in (A d ), that should now also include initial conditions.

E 2 mCRE (ê e , êp , ŝ; p)=E 2 CRE (ê e , êp , ŝ; p)+ ↵ 2 . 1 N DIC .N t Nt X nt=1 (d(û nt ) U nt DIC ) T ΣΣ 1 DIC (d(û nt ) U nt DIC ) with E 2 CRE (ê e , êp , ŝ; p)= Z T 0 Z Ω ✓ ⌘ ψ (ê e , ŝ; p)+ Z t 0 ⌘ ϕ ( ėp , ŝ; p) ◆ (
As indicated in previous chapters, the solution to the inverse problem corresponds to the parameter set p sol that satisfies the following nested minimization problem:

p sol = argmin p2P " min (êe,êp,ŝ)2(A d ) E 2 mCRE (ê e , êp , ŝ; p) # (5.20)
It is in practice solved by means of an alternated minimization scheme in which at iteration k +1:

•

Step 1: an optimal admissible set (ê

(k+1) e , ê (k+1)
p (k+1) = argmin p2P E 2 mCRE (ê (k+1) e , ê (k+1)
p , ŝ(k+1) ; p) (5.22)

Some model examples

We list here a set of elastoplastic constitutive models which will be used in numerical applications.

Prandtl-Reuss model

We first consider the Prandtl-Reuss plastic model with isotropic hardening, for which e e =[ ε e ,p] T , e p =[ ε p , p] T , and s

=[ σ,R] T , with p = R t 0 kε p kdt the cumulative inelastic strain (k • k =( • : •) 1/2
) and R the associated thermodynamic force (isotropic hardening variable on additional yield stress).

The associated free energy potential reads:

(ε e ,p)= 1 2 Kε e : ε e + g(p)

(5.23) with g a function that characterizes the hardening law (g(p)= 1 2 kp 2 for linear hardening, with k a strictly positive material parameter). We thus obtain the following state laws: σ = @ @ε e = Kε e ; R = @ @p = g 0 (p) (5.24) and the dual potential reads:

⇤ (σ,R)= 1 2 K 1 σ : σ + g ⇤ (R) (5.25)
with g ⇤ the Legendre-Fenchel transform of function g.

The dissipation potential ' ⇤ (σ,R) is the indicator function of the elasticity domain

C f = {(σ,R),z(σ,R)  0,R 0}
, that is:

' ⇤ (σ,R)=χ C f (σ,R)= 8 > > < > > : 0 if (σ,R) 2 C +1 if (σ,R) / 2
'(ε p , ṗ)=R 0 kε p k + Ψ Ce (ε p , ṗ) (5.27)
Remark . In the viscoplastic case (with linear hardening), the dissipation potentials should be changed to:

' ⇤ (σ,R)= k n +1 hzi n+1 + ; '(ε p , ṗ)=R 0 kε p k + k n 0 +1 ✓ kε p k k ◆ n 0 +1 + Ψ Ce (ε p , ṗ) (5.28)
with n 0 =1/n. The notation h•i + indicates the positive part, so that the elasticity domain is defined by hzi + =0.

When now considering (linear) kinematic hardening, the variable sets are complemented as e e =[ ε e , ↵↵ ,p] T , e p =[ ε p ↵↵ , p] T , and s =[ σ, X,R] T with X the kinematic hardening variable (backstress tensor) and ↵↵ the associated internal variable. The free energy is modified in: 5.29) and the von Mises yield function is rewritten as:

(ε e ,p)= 1 2 Kε e : ε e + 1 2 C↵↵ : ↵↵ + g(p) (
z(σ, X,R)=kσ D Xk + 1 2 AX : X (R + R 0 ) (5.30)
with A and C some constants. The normality rule then yields the following constitutive relations (yielding and hardening laws):

εp = ˙ σ D X kσ D Xk ; ↵↵ = ˙ ✓ X σ D kσ D Xk + AX ◆ ;ṗ = ˙ ; Ẋ = C(ε p AX ṗ); Ṙ = g 00 (p)ṗ (5.31) with ˙ = ⇣ σ D X kσ D Xk :σ ⌘ . ⇣ C( σ D X kσ D Xk AX):( σ D X kσ D Xk AX)+g 00 (p) ⌘ 1 .

Marquis-Chaboche model

This model, used in the case of several loading cycles [Ladevèze, 1999], is slightly different from the previous one.

It involves variable sets e e =[ ε e , ↵↵ , p] T , e p =[ ε p ↵↵ , p] T , and s =[ σ, X, r] T with r 0. In its viscoplastic version, the dissipation potential is defined as:

' ⇤ (σ, X, r)= k n +1 hzi n+1 + with z = kσ D Xk + 1 2 AkXk 2 (R(r)+R 0 )
k, n and A are positive material parameters, and R(r) is a concave increasing function. One thus gets the following evolution laws:

εp = ˙ σ D X kσ D Xk ; ↵↵ = ˙ ✓ X σ D kσ D Xk + AX ◆ ; ṗ = ˙ @R @ r with ˙ = khzi n + .
A priori considering Tr[ε p ]=0, the dual dissipation potential ' is given by:

'(ε p , ↵↵ , ṗ)=k n n +1 ✓ kε p k k ◆ 1+1/n + k ↵↵ εp k 2 2Akε p k + ( ṗ, kε p k) (5.32) with ( ṗ, kε p k)=kε p k(R 0 + R( ṗ/kε p k) ṗ[ ∂R ∂ r ] 1 ( ṗ/kε p k).
Elasto-plasticity is a limit case where n tends to +1; one then has:

' ⇤ (σ, X, r)=χ C f (σ, X, r) with C f = {(σ, X, r), hzi + =0} (elasticity domain) (5.33) and previous evolution laws hold with ˙ such that ˙ 0, z  0, and ˙ z =0.

Development of a specific solver

We now introduce a dedicated numerical strategy, taking into account the mathematical structure of the constructed nonlinear mCRE functional, to drive Step 1 of minimization process, that is to solve (5.21). This strategy, which shares similarities with the LATIN method [Ladevèze, 1989, Ladevèze, 1999], addresses material nonlinearities in the inverse problem. It enables to perform inversions in an effective way and with a reasonable CPU cost.

The LATIN method

The LATIN method was developed [Ladevèze, 1985,Ladevèze, 1989,Boisse et al., 1990,Ladevèze, 1999] in order to address the numerical complexity associated with nonlinear computational mechanics. As an alternative to classical step-by-step methods with decomposition of the loading path into small increments, the iterative LATIN solver takes into account the entire loading path in a single large time increment. It allows the effective solution to be performed for a large variety of material behaviors; it is used here for elastoplasticity.

The

(n) 2 (A d ) is known, a new solution s (n+1) 2 (A d
) is computed from the two following stages (see Figure 5.1):

• Local stage: from s (n) 2 (A d ), build s (n+1/2) 2 (Γ) such that s (n+1/2) s (n) 2 H + (i.e. (σ (n+1/2) σ (n) + H + (ε (n+1/2) ε (n) )=0
). H + is a positive define symmetric operator that can be interpreted as an ascent search direction; it may be chosen with infinite slope so that ε (n+1/2) = ε (n) (vertical direction). In practice, this local stage boils down to solving a nonlinear problem on [0,T] at any integration point in space, which is suited to parallel computing.

• Linear stage: from s (n+1/2) 2 (Γ), build s

(n+1) 2 (A d ) such that s (n+1) s (n+1/2) 2 H (i.e. σ (n+1) σ (n+1/2) H (ε (n+1) ε (n+1/2) )=0
). H is a positive define symmetric operator that can be interpreted as a descent search direction; it may be chosen as the tangent operator K, or such that H = K (elastic direction). This linear stage boils down to a global but linear problem, parametrized by time t.

The algorithm is stopped when a convergence criterion is reached; it may be based on a CRE measure (i.e. E CRE (s (n+1) )  " 0 ) or on a norm of s (n+1) s (n+1/2) . Convergence properties of the iterative algorithm can be found in [Ladevèze, 1999].

Remark . Another important feature of the LATIN method is the use of radial approximation (separated-variable representation, similar to POD and PGD) when solving the global linear problem parametrized by time t in the global stage. This enables for the use of the LATIN method as a suitable tool for building reduced models in nonlinear mechanics [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF], Relun et al., 2013, Néron et al., 2015, Ladevèze, 2016, Vitse et al., 2019]. This feature will not be investigated here.

Illustrative example

A beam in tension serves as an example illustrating how the LATIN method works. The uniaxial problem is given in We assume that the loading is quasi-static and cyclic, and that the material is viscoplastic and described by a variant (VP2) of the Marquis-Chaboche model. The strain is ✏ = u d /L with u d the prescribed displacement of the end of the beam and L its length. Besides the inelastic strain ✏ p , the internal variables are (↵,p) and their conjugates (X, R); they all are scalar variables.

The state equations read: .34) and we choose E = 137, 600 MPa, C = 24, 800 MPa, = 10, and a = 80 MPa.

✏ e = E ; X = C↵ ; R = a[1 exp( p)] ( 5
The evolution laws are described by the potential ' ⇤ such that:

' ⇤ ( ,X,R)= k n +1 hzi n+1 + ,z = | X| +(A/C)X 2 (R + R 0 ) (5.35) On the other hand, dz = X | X| d( X)+ 2A C XdX d  a R 1/2 2a (2 R 1/2 2a ) = X | X| d( X)+ 2A C XdX 1/2 (1 R 1/2 2a )dR d 2 z = 2A C dXdX + 2a dRdR (5.41) so that d 2 ' ⇤ = khzi n +  2A C (dX) 2 + 2a (dR) 2 + knhzi n 1 +  X | X| d( X)+ 2A C XdX 1/2 (1 R 1/2 2a )dR 2 (5.42)
Consequently, the descent direction is defined by

H = khzi n + 2 6 6 6 6 4 00 0 2A C 0 SY M γ 2a 3 7 7 7 7 5 +knhzi n 1 + 2 6 6 6 6 4
1 σ X |σ X| ( σ X |σ X| + 2A C X) σ X |σ X| (1 Rγ 1/2 2a ) 1/2 ( σ X |σ X| + 2A C X) 2 ( σ X |σ X| + 2A C X)(1 Rγ 1/2 2a ) 1/2 SY M (1 Rγ 1/2 2a ) 2 3 7 7 7 7 5
(5.43)

The problem is to find (✏

p , ↵ (1) ,p (1) , (1) ,X (1) ,R (1) ) defined on [0,T] satisfying:

✏ = E + ✏ p = U d (t) L ; X = C↵ ; R = ap ; 2 6 6 6 6 4 ✏p ✏(1/2) p (↵ ↵(1/2) ) (ṗ ṗ(1/2) ) 3 7 7 7 7 5 = H 2 6 6 6 6 4 (1/2) X X (1/2) R R (1/2) 3 7 7 7 7 5
(5.44)

with ✏ p = ↵ = p =0at t =0. (1/2)

Writing ✏ p = U d (t) L σ E ,
X (1/2) R (1/2) 3 7 7 7 7 5 2 6 6 6 6 4
✏(1/2) p Ud /L ↵(1/2) ṗ(1/2) 3 7 7 7 7 5
(5.45)

together with = X = R =0at t =0.
By integrating the previous differential equation, we construct (1) , X (1) and R (1) on [0,T]. Then, using the state equations, we obtain on [0,T]:

✏ (1) p = U d L σ (1) E , ↵ (1) = X (1) C ,p (1) = R (1) a (5.46)
The evolution of (✏

p , ↵ (1) ,p (1) , (1) ,X (1) ,R (1) ) on [0,T] is shown in Figure 5.5.

These two sets enable to reformulate the initial constrained minimization problem as two coupled minimizations in reduced spaces:

min Σψ 2(A ,ψ d ) E 2 ψ ( Σψ ; p); m i n Σϕ2(A ,ϕ d ) E 2 ϕ ( Σϕ ; p) (5.50)
The only variable which is computed in both minimizations is the stress tensor σ. invites to use a solution scheme which is similar to the LATIN method, with alternating solutions based on local and linear stages, in order to find the optimal set Σopt of admissible fields which is Σopt = S ψ \ S ϕ . Convergence of the scheme is ensured with similar arguments as in the LATIN method [Ladevèze, 1999, Section 4.5.3, Theorem 1].

Consequently, after introducing up and down directions H + and H , respectively, the following stages are iteratively performed:

• Local stage in which internal variables are searched at the integration point level (it can thus be parallelized):

given

Σ(n) ψ 2 S ψ , find Σ(n+1/2) ϕ 2 S ϕ such that: ( Σ(n+1/2) ϕ Σ(n) ψ ) 2 H + (5.51)
Among the various choices for H + , we consider here the one that prescribes the same strain tensor, that is ε(n+1/2) =ε (n) (so-called infinity direction in the literature, with vertical slope). This way, εp = εe .

• Linear stage associated with a global in space problem: given

Σ(n+1/2) ϕ 2 S ϕ , find Σ(n+1) ψ 2 S ψ such that: ( Σ(n+1) ψ Σ(n+1/2) ϕ ) 2 H (5.52)
The optimal choice for H is the local tangent direction to S ϕ (similarly to a classical Newton method) but this requires to compute a tangent operator at each iteration. For the sake of simplicity, we rather consider an elastic down direction that relates variations in stress and elastic strain tensors.

An illustration of the proposed new version of the LATIN algorithm within the mCRE framework is shown in Remark . A practical explanation of the proposed LATIN-based strategy is the description of the material integration step. In the described numerical approach, the material integration step is performed in two stages:

1. local integration: it corresponds to the integration of the evolution laws, i.e. the computation of the nonlinear evolutions. This stage only implies computations at the integration point level, so that this stage can be strongly parallelized. More precisely, at a given iteration, the entire evolution of the internal variables characterizing the nonlinearities is computed from the whole history coming from the linear space S ψ ;

2. linear integration: it corresponds to the evaluation of the state equations. This stage implies computations at the integration point level and an assembly stage to compute a global linear system . This step is performed knowing the entire evolution of the internal variables coming from the local integration step. The time aspect is taken into account by solving an incremental linear system explicitly.

Regarding

Step 2 of the mCRE strategy, which is the minimization (5.22) of the mCRE functional with respect to parameters (for a fixed set Σ of admissible fields), it is performed in a classical manner, with the adjoint-state method and a steepest-descent algorithm.

We mention that over the whole iterative mCRE minimization process, we can fully make benefit of the multiresolution aspect of the LATIN algorithm, with restart procedure [START_REF] Néron | Time-space pgd for the rapid solution of 3d nonlinear parametrized problems in the many-query context[END_REF]. Once material parameters p are updated, the previously computed space-time LATIN solution (optimal admissible fields obtained with p (k) ) are in practice reused as the initialization of another iterative LATIN procedure with the new set p (k+1) of material parameters (see Figure 5.9).

Figure 5.9: Graphical sketch of the employed multi-resolution procedure.

Implementation details

We give here some details on the implementation of linear and local stages in the LATIN algorithm associated with mCRE, considering a simple Prandlt-Reuss model with linear isotropic hardening.

Linear stage

In the linear stage, the following minimization is performed:

min Σψ 2(A ,ψ d ) ( Σψ Σ(n+1/2) ϕ )2H E 2 ψ ( Σψ ; p) (5.53)
Assuming that boundary conditions are reliable, this constrained minimization is performed with the following Lagrangian functional:

L( Σψ , ; p)=E 2 ψ ( Σψ ; p) Z T 0 Z Ω σ : ε( ) Z Ω f v d • Z ∂2Ω f s d • (5.54) with E 2 ψ ( Σψ ; p)= Z T 0 Z Ω ⌘ ψ ( Σψ ; p)+ ↵ 2 . 1 N DIC .N t Nt X nt=1 (d(û nt ) U nt DIC ) T ΣΣ 1 DIC (d(û nt ) U nt DIC ) (5.55) and ⌘ ψ ( Σψ ; p)= ✓ 1 2 K(p)(ε εp ):(ε εp )+ k 2 p2 ◆ + ✓ 1 2 K 1 (p)σ :σ + 1 2k R2 ◆ h σ :(ε εp )+ Rp i = 1 2 (σ K(p)(ε εp )) : K 1 (p)(σ K(p)(ε εp )) + 1 2k ( R k p) 2 0
(5.56)

Finding the saddle-point, with frozen internal variables, leads to the system:

Z T 0 Z Ω K(p)ε( û):(ε(û) εp )+↵. 1 N DIC .N t Nt X nt=1 d( ûnt ) T ΣΣ 1 DIC (d(û nt ) U nt DIC )=0 8 û Z T 0 Z Ω σ : K 1 (p)σ (ε εp ) ε( ) =0 8 σ Z T 0 Z Ω σ : ε( ) Z Ω f v d • Z ∂2Ω f s d • =0 8
(5.57)

After substituting σ = K(p)(ε + ε( ) εp ), and discretizing with FEM, we obtain a linear system of the form (2.9) which has to be solved at each time t 2 [0,T].

Local stage

In the local stage, the following minimization is performed:

min Σϕ2(A ,ϕ d ) ( Σϕ Σ(n) ψ )2H + E 2 ϕ ( Σϕ ; p) (5.58) with E 2 ϕ ( Σϕ ; p)= Z T 0 Z Ω Z t 0 ⌘ ϕ ( Σϕ ; p) (5.59) and ⌘ ϕ ( Σϕ ; p)=R 0 k εp k + χ Ce ( εp , ṗ)+χ C f (σ, R) σ : εp + R ṗ 0 (5.60) We recall that C e = {(ε p , ṗ), kε p k ṗ  0, Tr[ε p ]=0} and C f = {(σ,R), ||σ D k (R + R 0 )  0,R 0}
. This leads to a set of ODEs that drive the time evolution of internal variables at each integration point, in particular εp = ṗ σD kσ D k . In practice, they are solved with an implicit Euler scheme.

Model selection using error indicators

The strategy to get indicators on modeling and discretization errors is similar as that employed in Chapter 3. After tuning the weighting scalar ↵ in order to satisfy the Morozov principle, and after defining the reference energy E ref , an indicator on the mathematical model error ⌘ mod is computed from the CRE term of the mCRE functional, with the discretized admissible quantities: GPa, which are very close to the reference values (with relative errors of 1.5%, 2.5%, 6.3%, and 6.2%

⌘ 2 mod = E h 2 CRE (ê h e , êh p , ŝh ; p)= Z T 0 Z Ω ✓ ⌘ ψ (ê h e , ŝh ; p)+ Z t 0 ⌘ ϕ ( ėh p , ŝh ; p) ◆ ( 5
In order to confirm the interest of the modeling error indicator used in the adaptive model algorithm, we indicate in Table 5.1 the value of the relative indicator obtained at the end of the identification process from the 3 models considered. As expected, we observe that the indicator value is much larger than 1 when the model is not compatible with observations, while it is of the order of 1 when the reference model is selected.

Chapter 6

Extension to multiscale analysis with

MsFEM

In this last chapter, we introduce multiscale tools for identification from full-field measurements, using the Multiscale Finite Element Method (MsFEM). The proposed multiscale imaging-through-analysis approach enables to effectively capture fine scale details of the microstructure, in a computationally tractable procedure. The performance of the approach is shown on preliminary numerical experiments.

Motivation and basic concepts

Multiscale analysis

Materials with complex heterogeneous microstructures are commonly encountered in mechanical engineering.

Model identification and simulation in this context may be tricky, or even unfeasible with conventional approaches.

The development of multiscale numerical methods for such materials is therefore an active research field, and these methods are becoming a standard approach in material sciences and computational mechanics. As an alternative to solving a full fine-scale problem, with a usually prohibitive computational cost, multiscale modeling aims at linking the different scales for the accurate description of physical phenomena and/or the prediction of macroscopic (effective) properties. A main goal is to capture the impact of the smaller scales on the larger scales, in order to observe the influence of the microscopic structure on the macroscopic behavior.

In the last 20 years, much effort has been put in the design of multiscale approaches. They are inspired by the framework of the FEM but take into account scale separation between macroscopic and microscopic features. A pioneering work using multiscale finite element basis functions is reported in [Babuska et al., 1994a]. It is closely related to another development, the Generalized Finite Element Method (GFEM) [START_REF] Strouboulis | The generalized finite element method[END_REF]. These contributions share the idea of adapting the FE space to the particular fine-scale features of the problem by means of handbook functions and partition of unity. Other classical multiscale methods include the Multiscale Finite Element Method (MsFEM) [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF]] also based on modified basis functions obtained from fine-scale equations (see [START_REF] Efendiev | Multiscale Finite Element Methods: Theory and Applications[END_REF] for a review), the Variational Multiscale Method (VMS) [START_REF] Hughes | The variational multiscale method -a paradigm for computational mechanics[END_REF]], the sparse FEM [START_REF] Hoang | High-dimensional finite elements for elliptic problems with multiple scales[END_REF], multigrid methods [Hackbusch, 1985], and multilevel finite element methods such as the Heterogeneous Multiscale Method (HMM) [E et al., 2003] or the FE 2 method [Feyel, 2003]. We also mention the large literature on numerical homogenization from RVEs developed in the structural mechanics and engineering communities to obtain macroscopic constitutive laws [Suquet, 1987].

In the context of model identification from full-field measurements, multiscale approaches have been proposed [START_REF] Nguyen | Multiscale identification of the random elasticity field at mesoscale of a heterogeneous microstructure using multiscale experimental observations[END_REF], Passieux et al., 2015a, Zhang et al., 2020]. They are mainly based on a multiscale vision of images (multiscale field measurements), coupling macroscale and mesoscale analyses. Here, we propose an alternative approach that focuses on a multiscale numerical post-processing of image information. This is based on MsFEM which is a powerful tool to capture information at various scales without resorting to a full fine-scale computation [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF], Hou et al., 1999, Efendiev and Hou, 2009, Efendiev et al., 2013]. It defines an approximate solution in a finite dimensional space, related to a macroscopic mesh and generated by basis functions which encode details of the fine-scale heterogeneities. MsFEM then performs the computations in a two-stage procedure:

(i) an offline stage in which multiscale basis functions are computed solving local fine-scale problems; (ii) an online stage in which an inexpensive Galerkin approximation problem is solved from these basis functions.

The objective here is to use such a small set of multiscale basis functions that can effectively upscale fine-mesh behavior at the pixel size. It is in the spirit of approaches, mainly developed for crack analysis [START_REF] Roux | Stress intensity factor measurements from digital image correlation : post-processing and integrated approaches[END_REF], Réthoré et al., 2008, Roux et al., 2009, Rannou et al., 2010], in which specific basis functions are used to take into account the relevant mechanical information which is available. MsFEM is particularly relevant as it may be used with a quite coarse macro mesh, much decreasing the online computational cost and circumventing the issue on measurement uncertainty when using a fine mesh in DIC.

Basics on MsFEM

For simplicity reasons, we detail the MsFEM method on a fine-scale scalar elliptic problem (typically associated with heat transfer in composite materials or Darcy flow in porous media); the extension to elasticity, which will be considered in numerical results, is direct. The problem is characterized by strong material heterogeneities, with a characteristic scale denoted " (" ⌧ 1). The associated direct boundary value problem reads:

r • [A ε ru ε ]=f v d in Ω,u ε = u d on @ 1 Ω, (A ε ru ε ) • n = f s d on @ 2 Ω (6.1)
We assume that the body loading f v d 2 L 2 (Ω) and the traction loading f s d 2 L 2 (@ 2 Ω) are given and are slowly varying (they do not depend on "). The second-order diffusion tensor A ε 2 (L 1 (Ω)) d⇥d is a rapidly oscillating function. We assume that A ε is a symmetric matrix and that it is uniformly elliptic and bounded, in the sense that there exist ↵ > 0 and > 0 such that 8" 0, 8⇠ 2 R d , ↵|⇠| 2  (A ε (x)⇠) • ⇠  |⇠| 2 a.e. in Ω (6.2)

Introducing a mesh T H with characteristic size H, a conventional FE solution to the problem may be computed under the form u H = P i u i 0 i , where 0 i denote classical Lagrange FE basis functions. It is then well-known that obtaining an accurate approximation u H requires to choose H ⌧ " (as a priori error analysis yields ku ε u H k H 1 (Ω)  C(H/"), which leads in practice to a prohibitive computational cost.

From an engineering perspective, it may be sufficient to only predict the macroscopic (effective) properties of the solution to (6.1). Powerful homogenization methods can be used for this purpose, based on a scale separation assumption, i.e. considering the regime " ⌧ L where L is the characteristic size of the macroscopic variations of u ε inside Ω (see [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF], Sanchez-Palencia, 1980, Kanouté et al., 2009]. The idea is to approximate the behavior of the highly heterogeneous medium with an averaged behavior at a macroscopic scale. From a mathematical point of view, homogenization consists in identifying the limit of the operator r • [A ε r•] 1 when " tends to zero.

We consider in the following an alternative multiscale numerical approach, the MsFEM, which has the capability to approximate the full fluctuating solution u ε (with fine-scale details) in a convenient manner. Introduced in [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF], the main idea in the MsFEM approach is to construct a set { ε i } of local multiscale FE basis functions that encode small scale information within each element of a coarse mesh T H (H "). The basis functions ε i , associated with each node i of the coarse mesh T H (except those with prescribed Dirichlet boundary conditions), are adapted to the local properties of the operator. They are pre-computed in an offline stage, over each element K of the coarse mesh (see Figure 6.1), as solutions to local elliptic equations of the form

r • ⇥ A ε r ε i ⇤ =0 in K (6.3)
complemented with various boundary conditions discussed below. We note that the basis functions are independent of the loads f v d and f s d . In addition, problems (6.3) are decoupled one from each other. The MsFEM approach is thus well-suited to massively parallel computers; this is an essential feature of the approach.

Once the multiscale basis functions are computed, the MsFEM strategy is the same as that of a classical FE approach. It consists in performing, in an online stage, a Galerkin approximation with the finite dimensional space

V ε H = Span{ ε i }.
The small scale information incorporated in the basis functions is thus brought to the large scales through couplings in the global stiffness matrix K ε . The assembly of this matrix is inexpensive since it reuses local K K Figure 6.1: Left: a coarse mesh of the domain Ω is introduced, with elements of diameter H much larger than the small characteristic size " of the heterogeneities (here the diameter of the inclusions). Right: on each coarse element, we solve a local problem, using in practice a discretization at the size h adapted to the heterogeneities. matrices computed and stored in the offline stage to solve (6.3).

Remark . The philosophy of MsFEM, in which the fine-scale equations are solved inside coarser elements and are thus totally decoupled one from each other, is closely related to that of the VMS method [START_REF] Hughes | The variational multiscale method -a paradigm for computational mechanics[END_REF]] and the GFEM method [START_REF] Strouboulis | The generalized finite element method[END_REF]. However, these methods differ on several points, in particular in the definition of the local fine-scale problems (computational domain, boundary conditions).

The boundary conditions associated with (6.3) in the construction of basis functions are of critical importance, since they determine the behavior of the numerical approximation u ε H on element edges, and therefore significantly affect the accuracy of the MsFEM approach. In a conforming version of MsFEM, the choice is made to impose a linear evolution of ε i along @K as for classical first-order FE basis functions. Local problems then read: (6.4) Note that the support of ε i is identical to that of 0 i . In the regime of interest H ", and assuming that the oscillations in the material behavior are periodic (and with regularity assumptions), the following a priori error estimate is shown [START_REF] Hou | Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients[END_REF]:

r • ⇥ A ε r ε i ⇤ =0 in K, ε i = 0 i on @K
ku ε u ε H k H 1 (Ω)  C ⇣ H + p "/H + p " ⌘ (6.5)
where C is a positive constant independent of " and H. This result shows that the MsFEM approach converges to the correct solution in the homogenization limit (namely when first " ! 0 and next H ! 0, which is the regime of interest). The ratio "/H reflects the so-called cell resonance error, which occurs when H ⇡ ", and which is caused by the mismatch between the mesh size and the perfect sample size. In practice, this is often the leading error term.

In order to decrease the mismatch between u ε (which oscillates everywhere in Ω) and the MsFEM approximation u ε H (which does not oscillate in the vicinity of @K as a consequence of the non-oscillatory boundary conditions ε i = 0 i on @K), a possible strategy is to presribe (affine) Dirichlet boundary conditions not on @K, but farther away, on the boundary of a domain S K which is slightly larger than K, and to only use the interior information on K to construct the MsFEM basis functions (see Figure 6.2). This is the well-known oversampling technique that leads to a non-conforming MsFEM approximation since V ε H 6 ⇢ H 1 (Ω). It permits the boundary layer on which the small scale information is not transmitted to be reduced.

The approach thus consists in first solving the local problems

r • ⇥ A ε r ε,K i ⇤ =0 in S K , ε,K i
is affine on @S K , ε,K i (s j )= ij (6.6) where S K is a domain which is homothetic to the mesh element K, and s j the coordinates of the vertices of S K .

Second, the MsFEM basis functions are defined as ε i = ε,K i|K in K. In [START_REF] Efendiev | Convergence of a nonconforming multiscale finite element method[END_REF], an a priori error estimate is derived for the oversampling variant of the MsFEM approach. When H " and in the periodic setting, assuming that the distance between K and @S K is of order H, the estimate reads

ku ε u ε H k H 1 (Ω)  C 1 H + C 2 " H + C 3 p " (6.7)
where C 1 , C 2 and C 3 are independent of " and H, and C 2 depends on the oversampling ratio. Since the approach is non-conforming, the error is measured in the broken

H 1 -norm k • k H 1 (Ω) = ⇣ P K2T H k • k 2 H 1 (K) ⌘ 1/2 . S K S K Figure 6
.2: Illustration of the oversampling technique: a domain S K , which is slightly larger than K, is introduced to solve local fine scale problems.

Remark . A higher-order MsFEM approximation can also be constructed using local Kozlov harmonic coordinates [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF]]. In the case of a periodically oscillating coefficient A ε and with some regularity assumptions, an associated a priori error estimate was derived [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF]; it reads:

ku ε u ε H k H 1 (Ω)  C ⇣ H k + p "/H + p " ⌘ (6.8)
where k is the MsFEM approximation order, and C a constant independent of " and H.

Illustrative example

As an illustration of the MsFEM philosophy, we consider a one-dimensional problem defined in Ω =( 0 , 1), with homogeneous Dirichlet boundary conditions at x =0and x =1, and a non-uniform distributed load f (x)=x 2 .W e consider a non-periodic diffusion coefficient of the form A ε (x) = 5 + 50 sin 2 (⇡x 2 /") (6.9) and we set " =0.025. The evolution of A ε (x) is shown in Figure 6.3. We compute an approximate MsFEM solution using a coarse mesh T H composed of 5 elements with equal size (H K =0 .2). Oversampling is not needed in the one-dimensional context. We show in Figure 6.3 the computed MsFEM shape functions considering a fine mesh size h K = "/20 (identical over all coarse elements K).

The MsFEM solution, in terms of primal field u ε H but also gradient du ε H /dx and flux A ε du ε H /dx, is reported in Figure 6. The objective is to identify the value of E 0 , by means of artificial measurements derived from a reference solution computed with E 0 =1 . This reference solution, numerically obtained from an overkill mesh (with 500 ⇥ 500 Q4 elements), is shown in Figures 6.7, 6.8, and 6.9. When starting the identification process by using MsFEM, an initial coarse mesh T H made of 5 ⇥ 5 macro elements is considered. Choosing h K = H K for any K, and with no oversampling, the optimal admissible MsFEM solution is represented in Figure 6.10.

At the end of the identification process with the initial coarse MsFEM model, the error indicator ⌘ dis on discretization error indicates that this error source is too large, so that adaptive refinement is performed. In Figure 6.11, we represent a MsFEM shape function used at the beginning of the identification process, and how it is changed after applying the adaptive strategy. In Figure 6.12, we show the distribution of the initial error indicator over the macro mesh, as well as the final configuration obtained for discretization parameters h K (related to the macro mesh), d K (related to oversampling) and h K (related to the micro mesh) at the end of the adaptive process. This configuration is obtained after 8 iterations, and it results from an optimized compromise between model accuracy and measurement uncertainties.

It is found that the discretization is mainly refined in regions located close to the boundary of Ω, while inexpensive computations with coarse discretization are sufficient in zones where the fine-scale features of the solution are not activated.

Eventually, we show in Figure 6.13 the evolution of the ratio E/E 0 between the identified value of E and the reference value E 0 , after each adaptive iteration. We observe that the E is correctly identified (less than 5% error) with the optimized model, after starting from a large error (higher than 30%). computed with E 0 = 64 9 p 17 (so that E =1). This reference solution, numerically obtained from an overkill mesh (with 500 ⇥ 500 Q4 elements), is shown in Figures 6.15,6.16,and 6.17. When starting the identification process by using MsFEM, an initial coarse mesh T H made of 8 ⇥ 8 macro elements is considered. At the end of the identification process with the initial coarse MsFEM model, the error indicator ⌘ dis on discretization indicates that this source is too large, so that adaptive refinement is performed. In Figure 6.18, we show the distribution of the error indicator over the macro mesh, as well as the final configuration obtained for discretization parameters h K (related to the macro mesh), d K (related to oversampling) and h K (related to the micro mesh) at the end of the adaptive process. This configuration is obtained after 6 iterations, and it results from an optimized compromise between model accuracy and measurement uncertainties. We observe the adaptive refinement of the macro mesh T H is performed close to the crack tip, as traditional scale separation assumptions no longer hold in the vicinity of the crack tip (high variation of the gradients), while a coarser discretization is used in the remainder of the physical domain.

Eventually, we show in Figure 6.19 the evolution of the ratio E/E 0 between the identified value of E and the reference value E 0 , after each adaptive iteration. We observe that the E is correctly identified (less than 5% error) with the optimized model, after starting from a large error (higher than 20%). Remark . For this example, it would of course be interesting to use a MsFEM version of the XFEM method, which would consists in defining the MsFEM basis functions by means of local problems (posed on each coarse element

Conclusions & prospects

We developed a general methodology around the duality-based CRE concept in order to produce robust and effective tools for model updating and selection from full-field measurements. We first proposed a consistent definition of the mCRE functional, in association with DIC or in an integrated version, with a measurement term that is properly scaled to take measurement noise into account. We thus defined a relevant metric in which both measurement and modeling errors are mixed. We then focused on modeling error by developing indicators on various error sources (model class, discretization), and we designed an associated adaptive algorithm in order to compute right, and at the right cost, with regards to experimental information. In this context, we implemented ROM tools in order to effectively address multi-query procedures and enhance the overall performance. We also proposed an extension to nonlinear models by referring to thermodynamical bases and Legendre-Fenchel duality, as well as by designing a specific solver. Last, a multiscale approach was introduced in order to address complex material descriptions at the micro scale.

We think this PhD work brings several scientific advances in the context of model identification from full-field measurements, and that it strengthens the link between material modeling, numerical methods, and experimental techniques. It has also some potential for industrial applications related to structural health monitoring, as it offers a relevant strategy to manage modeling and simulation in association with experimental data.

As it provides for optimized numerical strategies in the identification procedure, that may be useful in many contexts (e.g. to address complex material descriptions, large experimental data, or real-time applications), the work naturally invites for further studies and research prospects.

A very first one is the validation of the interest of the approach on more practical situations with experimental data. It is also the use of the proposed identification approach for complex multiscale behaviors, in particular anisotropic damage in composites. In this context, an adaptive representation of the damage parameter field could be investigated using tools on modeling error estimation and adaptive mesh to represent this parameter field (as used in [START_REF] Bangerth | Adaptive finite element methods for the solution of inverse problems in optical tomography[END_REF], Bangerth, 2008, Puel and Aubry, 2011, Buhan and Darbas, 2017] for other contexts of biomechanical imaging or electromagnetics). This would thus lead to a multiscale updating strategy in which a global-scale coarse computation (with global low-resolution sensor information) would be first performed, before locally zooming and considering additional higher-resolution data and models for further analysis in regions Eventually, a fundamental aspect which needs to be further analyzed is in the practical choice of a mathematical model. From the associated error indicator, we here picked a richer mathematical empirically, but an objective could be to add more insight into the selection of this new model. This may be performed with a proper filtering and deeper analysis of the CRE term (e.g. using IA and machine learning algorithms [START_REF] Russel | Artificial Intelligence: a modern approach[END_REF]), in order to better understand the origins of modeling error. This may also be coupled with information theory, which indicates that for similar modeling error levels, the best model is the one associated with lowest a priori input. In any case, this refers to a smart learning of model ignorance from data, using the physics-based mCRE concept in complement to other advanced techniques such as neural networks [START_REF] Fang | Structural damage detection using neural network with learning rate improvement[END_REF], Raissi et al., 2019]. In this context, the recent attempts for full or highly data-based modeling and simulation from data mining techniques [START_REF] Le | Computational homogenization of nonlinear elastic materials using neural networks[END_REF], Brunton et al., 2016, Kirchdoerfer and Ortiz, 2016, Bessa et al., 2017, Kirchdoerfer and Ortiz, 2017, Leygue et al., 2018, Lopez et al., 2018, Gonzalez et al., 2019, Ladevèze et al., 2019] deserve to be analyzed in association with mCRE, in order to reconsider the issue of the constitutive law representation and introduce the least a priori constraints and empirism in the identification procedure. Titre: Nouvelles stratégies numériques pour l'identification robuste, cohérente et efficace de modèles à partir de mesures de champs Mots clés: Mesures de champs, Problèmes inverses, Erreur de modèle, Techniques adaptatives, Réduction de modèle, Analyse multiéchelle

Résumé:

Le travail de recherche a pour ambition l'identification robuste, cohérente, et efficace de paramètres matériaux à partir de mesures de champs obtenues par les techniques expérimentales avancées basées sur la corrélation d'image numérique (CIN). Dans ce contexte, la simulation numérique joue un rôle majeur car des calculs éléments finis sont effectués en parallèle de l'acquisition des données expérimentales afin de mener la procédure d'identification. L'objectif de la thèse est d'étudier et d'optimiser cette procédure en définissant une architecture numérique dédiée en terme de choix de modèle mathématique et du maillage, en fonction de l'information et du bruit contenus dans les données expérimentales. Ce choix vise à calculer juste au juste coût. L'idée est de relier les diverses sources d'erreur présentes dans la procédure d'identification: 1) bruit de mesure lié à la résolution des images; 2) erreur de modèle venant de la sélection d'un modèle basé sur la physique qui reste une représentation imparfaite de la réalité; 3) erreur de discrétisation due à l'utilisation de maillages générant des approximations numériques. On cherche donc à définir des modèles et des maillages dont le coût numérique et la précision sont gouvernés par l'information expéri-mentale disponible, avec des niveaux d'erreur similaires. Pour cela, on propose une méthodologie basée sur l'analyse duale, issue de développements du laboratoire depuis de nombreuses années, et valide pour une large gamme de modèles linéaires ou nonlinéaires. L'outil d'identification obtenu, basé sur la fiabilité de l'information, est construit à partir d'une formulation mathématique hybride avec une fonction coût composée de termes liés au modèle ou aux observations. On montre qu'il permet de définir des indicateurs d'erreur spécifiques et de mettre en place une procédure adaptative afin de corriger les biais de la modélisation numérique et la rendre cohérente avec les caractéristiques des mesures de champs. Un aspect majeur de la procédure est le calcul de champs admissibles vérifiant l'information pertinente. L'approche proposée est dans un premier temps validée sur des modèles d'élasticité linéaire avant d'être appliquée avec des comportements matériau non-linéaires afin de mener le processus de sélection et d'adaptation de modèle. De plus, l'approche est rendue plus efficace par l'utilisation complémentaire de méthodes numériques avancées telles que la réduction de modèle, des solveurs non-linéaires spécifiques, ou l'analyse multi-échelle.

I

  will never forget all the chats, laughs and pinches along the way with my lab companions, fellow PhD candidates, professors, administrative coordinators of the Laboratory of Mechanics and Technology (LMT) at Ecole Normale Supérieure Paris-Saclay (Gif-sur-Yvette, France). Notably, I could not start this research journey without full financial support from the Vietnam International Education Development (VIED), under the Ministry of Education and Training (MoET), and the Consultancy Firm for Civil Engineering Design and Quality Verification (CFCE) under the Water Resources University (Hanoi, Vietnam).

1

  General overview of V&V activities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Configuration of an inverse problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Configuration of the reference problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Geometrical representations of the CRE concept. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Geometrical representation of the CRE measure for nonlinear material behaviors. . . . . . . . . . . . 1.5 Geometrical interpretation of the CRE concept in inverse problems, with (A + d ) \ (Γ p )=;........ 1.6 Geometrical interpretation of the mCRE strategy, with (A d ) \ (Γ +obs p )=;.. . . . . . . . . . . . . . . . 1.7 Ill-posed problem on an elastic bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Schematic view of FEMU identification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Different steps of weighted FEMU-U and IDIC-U. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2. 27

 27 Used boundary conditions for FEMU and mCRE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.28 Identification results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.29 Speckle on the cross-shaped specimen (left), and FEM discretization for DIC with coarser (center) and finer (right) meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.30 Horizontal (left) and vertical (right) components of the measured displacement field when using the coarser DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.31 Morozov curve obtained from mCRE and the coarser DIC mesh. . . . . . . . . . . . . . . . . . . . . . 2.32 Evolutions of the mCRE functional and its two components with =0.. . . . . . . . . . . . . . . . . 2.33 Spatial representation of admissible fields û (left), v (center), and their difference (right), for =0 and with the coarser DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.34 Spatial representation the measurement field (left), the admissible field û (center), and their difference (right), for =0and with the coarser DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.35 Horizontal (left) and vertical (right) components of the measured displacement field when using the finer DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.36 Morozov curve obtained from mCRE and the finer DIC mesh. . . . . . . . . . . . . . . . . . . . . . . 2.37 Evolutions of the mCRE functional and its two components with = 0.5.. . . . . . . . . . . . . . . 2.38 Spatial representation of admissible fields û (left), v (center), and their difference (right), for = 0.5 and with the finer DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.39 Spatial representation the measurement field (left), the admissible field û (center), and their difference (right), for = 0.5 and with the coarser DIC mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 2 2 CRE

 22 Geometrical illustration of the proposed model selection strategy from mCRE. . . . . . . . . . . . . . . 3.3 Optimal configuration which is searched at the end of the adaptive process, with normalized components on modeling E (in blue) and measurement 1 N DIC (d(û) U DIC ) T ΣΣ 1 DIC (d(û) U DIC ) (in red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Scheme of the adaptive process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Reference model with heterogeneous material. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Normalized error estimator and indicators at each iteration of the adaptive algorithm. . . . . . . . . . . 3.7 Configurations obtained along the six iterations of the adaptive process (from left to right): when the model is adapted, highest contributions of the estimate ⌘ 2 mod are indicated in red. . . . . . . . . . . . . 3.8 Distribution of the indicator on discretization error ⌘ 2 dis (left), and mesh adaptation process with marking of elements (center and right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Illustration of the oversampling technique: a domain S K , which is slightly larger than K, is introduced to solve local fine scale problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6. 7

 7 Components of the reference displacement field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Components of the reference strain field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Components of the reference stress field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Components of the initial admissible displacement field. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Nodal MsFEM shape function at iteration 0 (left) and after adaptation of the local micro meshsize. . . 6.12 Distribution of the discretization error indicator at first iteration (top left) and final configuration after applying the adaptive process in terms of macro mesh (top right), oversampling size (bottom left), and micro mesh (bottom right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Value of E/E 0 along the adaptive modeling iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 Configuration of the crack domain (left), and evolution of E(x) in this domain (right). . . . . . . . . . . 6.15 Components of the reference displacement field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 Components of the reference strain field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 Components of the reference stress field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.18 Distribution of the discretization error indicator at first iteration (top left) and final configuration after applying the adaptive process in terms of macro mesh (top right), oversampling size (bottom left), and micro mesh (bottom right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.19 Value of E/E 0 along the adaptive modeling iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: General overview of V&V activities.

Figure 1 . 1 :

 11 Figure 1.1: Configuration of an inverse problem.

Figure 1 . 2 :

 12 Figure 1.2: Configuration of the reference problem.

Figure 1 . 3 :

 13 Figure 1.3: Geometrical representations of the CRE concept.

Figure 1 . 4 :

 14 Figure 1.4: Geometrical representation of the CRE measure for nonlinear material behaviors.

Figure 1 . 5 :

 15 Figure 1.5: Geometrical interpretation of the CRE concept in inverse problems, with (A + d ) \ (Γ p )=;.

  can again be geometrically interpreted as shown in Figure 1.6; the minimal distance is searched between an admissible solution (û, σ) 2 (A d ) and a manifold (Γ +obs p ) generated by the parametrized constitutive model and noisy observations. Here again, (A d ) \ (Γ +obs p )=; in most situations.

Figure 1 . 6 :

 16 Figure 1.6: Geometrical interpretation of the mCRE strategy, with (A d ) \ (Γ +obs p )=;.

Figure 1

 1 Figure 1.7: Ill-posed problem on an elastic bar.

  DIC (x)) is the updated deformed image at iteration k, and u (k) = u (k+1) DIC u (k) DIC is the correction at iteration k +1.

  (k) DIC (1.56) with M DIC = P x2ROI N T rf r T f N the symmetric semi-positive DIC matrix (computed once for all), and b (k) DIC = P x2ROI (f g)N T rf the residual vector at iteration k.

Figure 1 .

 1 Figure 1.10: Schematic view of FEMU identification.

Figure 1 .

 1 Figure1.12: Tensile test on dog-bone sample with a hole (left), and associated FE mesh with measured and extruded Dirichlet boundary conditions (right), from[START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF].

Figure 1 .

 1 Figure1.13: Global DIC results for the analyzed test: displacement fields expressed in pixels (top and center left) and normalized gray level residual with respect to the dynamic range of the picture in the reference configuration (bottom left); corresponding in-plane strain fields (right), from[START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF].

Figure 1 .

 1 Figure1.15: I-DIC residuals for the analyzed test: displacement residuals, normalized gray level residuals, and load residuals for the first (left) and converged (right) solutions, from[START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF].

Figure 2 . 1 :

 21 Figure 2.1: Reference model related to the studied identification problem.

2 0

 2 allows one to determine whether the model is valid or not. If one can not reduce the mCRE error below the required quality level by acting on model parameters, this indicates that errors are present in the model itself, not only in its parameters (see Chapter 3). In such a case, a stagnation criterion should be preferred to stop the iterative algorithm of the inversion process.

Figure 2 . 2 :

 22 Figure 2.2: Schematic view of the DIC-mCRE identification process.

Figure 2 .

 2 Figure 2.3: Schematic view of mI-DIC identification.

Figure 2 . 4 :

 24 Figure 2.4: Model problem with tensile test on a square plate.

Figure 2 . 5 :

 25 Figure 2.5: mCRE cost function E h 2mCRE and its two components as a function of E/E 0 , for various noise levels (from 0% to 5%).

Figure 2 . 6 :

 26 Figure 2.6: mCRE cost function E h 2mCRE and its two components as a function of E/E 0 , for a 20% noise level.

Figure 2 .

 2 Figure 2.7: mCRE cost function E h 2mCRE and its two components as a function of E/E 0 , for a 50% noise level.

Figure 2 . 8 :

 28 Figure 2.8: Synthetic displacement data on a pixel-like fine grid, with horizontal (left) and vertical (center) components, and sub-pixel discretization to compute the deformed images (right).

Figure 2 . 9 :

 29 Figure 2.9: Mean value and standard deviation of E/E 0 with respect to the noise level.

Figure 2 .

 2 Figure 2.10: Model problem with tensile test on a plate with hole (left), and associated mesh (right).

Figure 2 .

 2 Figure 2.11: Reference solution used to obtain synthetic data for identification.

Figure 2 .

 2 Figure 2.12: Noise distribution (left) and resulting noisy displacement field used as measurements (right), for a 10% noise level.

ure 2 . 1 DIC = 1 σ 2 R

 212 23), to which Gaussian noise is added. In the present case, the noise is correlated with Σ Σ ROI N T N (mean field assumption). A measurement grid of size 30 ⇥ 70 is then created, over which all integral computations are performed.

Figure 2 .

 2 Figure 2.13: Evolution of the term on discrepancy with measurements as a function of , for a 10% noise level.

Figure 2 .

 2 Figure 2.14: Evolutions of the mCRE functional and its two components with =0, for a 10% noise level.

Figure 2 .

 2 Figure 2.15: Spatial representation of admissible fields û (left), v (center), and their difference (right), for a 10% noise level.

Figure 2 .

 2 Figure 2.16: Spatial representation the admissible field û (left), measurement field (center), and their difference (right), for a 10% noise level.

Figure 2 .

 2 Figure 2.17: Noise distribution (left) and resulting noisy displacement field used as measurements (right), for a 30% noise level.

Figure 2 .

 2 Figure 2.18: Evolution of the term on discrepancy with measurements as a function of , for a 30% noise level.

Figure 2 .

 2 Figure 2.19: Evolutions of the mCRE functional and its two components with =0, for a 30% noise level.

Figure 2 .

 2 Figure 2.20: Spatial representation of admissible fields û (left), v (center), and their difference (right), for a 30% noise level.

Figure 2 .

 2 Figure 2.21: Spatial representation the admissible field û (left), measurement field (center), and their difference (right), for a 30% noise level.

Figure 2 .

 2 Figure 2.22: Model problem with three-point bending test on a beam.

Figure 2 .

 2 Figure 2.23: Reference calculation.

Figure 2 .

 2 Figure 2.24: Evolution of the mCRE cost function for various noise levels: 5%, 10%, 20% (from left to right).

Figure 2 .

 2 Figure 2.25: Analyzed test cases with restricted measurement zone.

Figure 2 .

 2 Figure 2.26: Identification with a restricted measurement zone.

Figure 2 .

 2 Figure 2.27: Used boundary conditions for FEMU and mCRE.

Figure 2 .

 2 Figure 2.28: Identification results.

Figure 2 .

 2 Figure 2.29: Speckle on the cross-shaped specimen (left), and FEM discretization for DIC with coarser (center) and finer (right) meshes.

Figure 2 .

 2 Figure 2.30: Horizontal (left) and vertical (right) components of the measured displacement field when using the coarser DIC mesh.

Figure 2 .

 2 Figure 2.31: Morozov curve obtained from mCRE and the coarser DIC mesh.

Figure 2 .

 2 Figure 2.32: Evolutions of the mCRE functional and its two components with =0.

Figure 2 .

 2 Figure 2.33: Spatial representation of admissible fields û (left), v (center), and their difference (right), for =0and with the coarser DIC mesh.

Figure 2 .

 2 Figure 2.34: Spatial representation the measurement field (left), the admissible field û (center), and their difference (right), for =0and with the coarser DIC mesh.

Figure 2 .

 2 Figure 2.35: Horizontal (left) and vertical (right) components of the measured displacement field when using the finer DIC mesh.

Figure 2 .

 2 Figure 2.36: Morozov curve obtained from mCRE and the finer DIC mesh.

Figure 2 .

 2 Figure 2.37: Evolutions of the mCRE functional and its two components with = 0.5.

Figure 2 .

 2 Figure 2.38: Spatial representation of admissible fields û (left), v (center), and their difference (right), for = 0.5 and with the finer DIC mesh.

Figure 2 .

 2 Figure 2.39: Spatial representation the measurement field (left), the admissible field û (center), and their difference (right), for = 0.5 and with the coarser DIC mesh.

Figure 3 . 2 :

 32 Figure 3.2: Geometrical illustration of the proposed model selection strategy from mCRE.

Figure 3

 3 Figure 3.3: Optimal configuration which is searched at the end of the adaptive process, with normalized components on modeling E 2 CRE (in blue) and measurement 1 N DIC (d(û) U DIC ) T ΣΣ 1 DIC (d(û) U DIC ) (in red).

Figure 3 . 4 :

 34 Figure 3.4: Scheme of the adaptive process.

Figure 3 . 5 :

 35 Figure 3.5: Reference model with heterogeneous material.

Figure 3 . 6 :

 36 Figure 3.6: Normalized error estimator and indicators at each iteration of the adaptive algorithm.

Figure 3 . 7 :

 37 Figure 3.7: Configurations obtained along the six iterations of the adaptive process (from left to right): when the model is adapted, highest contributions of the estimate ⌘ 2 mod are indicated in red.

Figure 3

 3 Figure 3.8: Distribution of the indicator on discretization error ⌘ 2 dis (left), and mesh adaptation process with marking of elements (center and right).

Figure 4 .

 4 Figure 4.1: PGD modal decomposition of the solution for a transient thermal problem: space functions (top), time functions (bottom).

Figure 4 .

 4 Figure 4.3: Relative energy norm of each PGD mode.

Figure 4 . 4 :

 44 Figure 4.4: Configuration at the level of a 3-node triangle element, with linear tractions fK|Γ defined on element edges.

Figure 4 . 5 :

 45 Figure 4.5: Mapping between reference and parametrized elements for a 3-node triangle element.

Figure 4 . 6 :

 46 Figure 4.6: Flow chart of the hybrid-flux approach without (left) or with (right) the use of PGD in Step 2.

Figure 4 .

 4 Figure 4.7: First PGD modes of ⇢ 1xΓ12,m : space functions i (x) (top), and geometry parameter functions  x,i (x 3 ) (center) and  y,i (y 3 ) (bottom), for i =1, 2, 3 (from left to right).

Figure 4 . 8 :

 48 Figure 4.8: Representation of the approximate PGD solution ⇢ 1x Γ12,m for various geometry configurations.

p

  , ŝ(k+1) ) is searched by solving the first partial minimization over (A d ), for given p = p (k) :(ê (k+1) e , ê(k+1) p , ŝ(k+1) ) = argmin (êe,êp,ŝ)2(A d ) E 2 mCRE (ê e , êp , ŝ; p (k) ) (5.21) • Step 2: a new parameter set p (k+1) is computed by solving the second partial minimization over P, for the previous admissible fields:

  kσ D k (R+R 0 ) for linear hardening and z = kσ D k ( 1 α+1 R α+1 +R 0 ) for power hardening, σ D the deviatoric part of the stress tensor, and R 0 0 the yield stress. Introducing the convex set C e = {(ε p , ṗ), kε p k ṗ  0, Tr[ε p ]=0 } with associated indicator function Ψ Ce , the dual dissipation potential reads (for linear hardening):

Figure 5 . 1 :

 51 Figure 5.1: General frame of the LATIN algorithm.

Figure 5

 5 Figure 5.2.

Figure 5 . 2 :

 52 Figure 5.2: A uniaxial example.

Figure 5

 5 Figure 5.8.

Figure 5 .

 5 Figure 5.8: New illustration of the LATIN algorithm in the mCRE framework.

Figure 5 .

 5 Figure 5.11: Relative evolution of stress-strain components ( xx , ✏ xx ) during the loading history.

Figure 5 .

 5 Figure 5.12: Maps of the cumulative plastic strain field at three loading increments.

Figure 5 .

 5 Figure 5.13: Admissible displacement field obtained from an elasticity model, with horizontal (left) and vertical (right) components.

Figure 6 . 3 :

 63 Figure 6.3: Evolution of A ε for the considered 1D non-periodic case (left), and MsFEM basis functions (right).

4 .Figure 6 . 4 :

 464 Figure 6.4: Comparison between the exact solution and the approximate MsFEM solution: primal fields u ε H and u ε (left), gradients du ε H /dx and du ε /dx (center), and fluxes A ε du ε H /dx and A ε du ε /dx (right).

Figure 6 . 6 :

 66 Figure 6.6: Artificial speckle considered (left), and evolution of E(x) in the domain (right).

Figure 6 . 7 :

 67 Figure 6.7: Components of the reference displacement field.

Figure 6 . 8 :

 68 Figure 6.8: Components of the reference strain field.

Figure 6 . 9 :

 69 Figure 6.9: Components of the reference stress field.

Figure 6 .

 6 Figure 6.12: Distribution of the discretization error indicator at first iteration (top left) and final configuration after applying the adaptive process in terms of macro mesh (top right), oversampling size (bottom left), and micro mesh (bottom right).

Figure 6 .Figure 6 .

 66 Figure 6.13: Value of E/E 0 along the adaptive modeling iterations.

Figure 6 .

 6 Figure 6.15: Components of the reference displacement field.

Figure 6 .

 6 Figure 6.16: Components of the reference strain field.

Figure 6 .

 6 Figure 6.17: Components of the reference stress field.

Figure 6 .

 6 Figure 6.18: Distribution of the discretization error indicator at first iteration (top left) and final configuration after applying the adaptive process in terms of macro mesh (top right), oversampling size (bottom left), and micro mesh (bottom right).

Figure 6 .

 6 Figure 6.19: Value of E/E 0 along the adaptive modeling iterations.

  1 N DIC (Π T ΣΣ 1 DIC Π) aa . This value of the CRE term statistically tends to a scalar value denoted E 2 ref (↵). For a given ↵, E 2 ref (↵) merely depends on measurement noise statistics and on the considered model; it can thus be computed in an offline phase (using for instanceMonte-Carlo samplings). The value E 2 ref (↵) corresponds to the reference value of the CRE term due to measurement uncertainties, taking into account the impact of measurement noise on admissible fields through a model assumed to be perfect. It is a value of choice for properly scaling the mCRE functional (see below).

  m (x ref , p geo , p)+ f jy K|Γ ⇢ jy Γ,m (x ref , p geo , p) T (p geo )rr ] s ⇢ m (x ref , p geo , p, {

		i	(4.32)
	It provides a parametrized equilibrated stress field at the element level:
	σm|K (x ref , p geo , p, {	f jx/y K|Γ })=K(p)[J f jx/y K|Γ })	(4.33)

  LATIN method makes use of properties of equations of a direct nonlinear mechanical problem, which can be split into two groups: linear (potentially global in space) equations, such as equilibrium or boundary conditions, defining a linear manifold (A d ) (admissibility space); local in space (potentially nonlinear) equations, such as constitutive equations, defining a manifold (Γ). The definitions of (A d ) and (Γ) do not exactly correspond to those introduced in the CRE context (as linear state equations are here placed in (A

d

)). The initial problem is then recast as finding s =( ε, σ) 2 (Γ) \ (A d ) as the converged solution of an iterative strategy that avoids solving a global time-space nonlinear problem all at once. An iteration is divided into a local stage and a linear stage, each computing an approximate solution over the whole space-time interval. Starting from an initial elastic calculation and assuming that a solution s

  (using Lagrange multiplier fields) is thus associated with a problem which is global over the whole spacetime domain but which is linear. The solution manifold corresponds to the linear space S ψ = {argmin Σψ 2(A ,ψ The constraint Σϕ 2 (A ,ϕ d ) involves initial conditions alone. The constrained minimization of E 2 ϕ is thus associated with a problem which is local in space (and global in time), but which is nonlinear. The solution manifold corresponds to the nonlinear space S ϕ = {argmin Σϕ2(A ,ϕThe previous split of the mCRE functional thus leads to a formulation with the solution to: (i) a linear and global in space problem; (ii) a nonlinear and local in space problem. Both problems are global in time. This formulation

	The constraint Σψ 2 (A ,ψ d	) involves equilibrium equations and boundary conditions. The constrained minimiza-
	tion of E 2 ψ d	) E 2 ψ ( Σψ ; p)}.
		d	) E 2 ϕ ( Σϕ ; p)}.

Table 5 .

 5 1: Relative indicator on modeling error obtained with the three considered models. k =0.015 MPa. When using the LATIN method in order to recover optimal admissible fields, about 15 sub-iterations are necessary for first iterations of the mCRE minimization, while between 2 to 5 sub-iterations are used for other minimization steps; this is due to the use of previously computed admissible fields to initialize the process (restart procedure). The identification procedure is performed with 8 iterations.

Last, we change the model and consider the Prandtl-Reuss plastic model with kinematic hardening (i.e. the one used to get measurements). The identified values are E = 203 GPa, R 0 = 164 MPa, A =0 .59 MPa 1 , and C = 45

Acknowledgements

"If we knew what it was we were doing, it would not be called research, would it?"

List of Tables

Enhanced numerical framework using reduced order modeling

The numerical approach introduced in previous chapters is associated with the computations of admissible fields for both identification (within the mCRE functional) and mesh adaptation (within the CRE functional). These computations are multi-query in the sense that they require many solutions of the same problem, but with various parameter values (with iterative procedure, or change in the computational domain geometry). The multi-query aspects imply computationally intensive procedures with a large number of evaluations of a deterministic finite element model. In order to circumvent this issue, we investigate in this chapter the use of reduced order modeling (ROM) techniques, and in particular the PGD technique, associated with an offline-online strategy. We show that these can be advantageously coupled to the proposed mCRE inversion approach in order to facilitate implementation and decrease computational cost, thus leading to a faster and cheaper V&V procedure.

Reduced order modeling with PGD

Basis on ROM

Applied mathematics, computational mechanics, and computer sciences contributed in the last two decades to new modeling and simulation procedures in which reduced-order modeling (ROM) techniques are one of the major achievements. These advanced techniques address complex high-dimensional engineering problems, with a large set of parameters, which are out of reach or remain very costly despite constant enhancements in computing resources. The issue comes from the exponential growth of complexity when using grid-based discretization strategies; this is the so-called curse of dimensionality. Conversely, ROM provides for surrogate simulation models that are computationally manageable while retaining high levels of accuracy on predicted outputs.

Numerical results

We illustrate the approach by considering again the test-case of Section 3. When implementing the mCRE identification process, an overall speed-up of 56 is obtained in the online phase.

Construction of equilibrated fields using PGD

Presentation

The hybrid flux (or EET) equilibration technique [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF], Ladevèze and Maunder, 1996, Pled et al., 2011] is used in Section 3.2.2 in order to recover the indicator ⌘ dis on the discretization error in the mCRE identification process. It consists in post-processing the stress field σ( V)=Kε(N V) at hand (equilibrated in a In order to capture the discretization error, a fully equilibrated stress field is recovered using the hybrid-flux technique with direct post-processing of the FE stress field at hand (at each loading increment point). This enables an estimate ⌘ tot on the whole modeling error to be obtained, then an indicator ⌘ 2 dis = ⌘ 2 tot ⌘ 2 mod on the discretization error.

The quantities (⌘ tot , ⌘ mod , ⌘ dis ) are further used to drive model selection and adaptation with respect to the observed data, and among a hierarchy of possible nonlinear models with increasing complexity (multi-fidelity approach). This is performed by means of an adaptive greedy algorithm similar to the one detailed in Section 3.3.

Starting from a coarse model, higher-fidelity models are progressively inserted when dictated by data.

Remark . The choice of a new higher-fidelity mathematical model is here performed quite empirically among a list of possible models. It may be possible to design a more elaborate procedure in which admissible fields (which are independent, a priori, of any constitutive equation) are tested with several CRE functionals, in order to select the best one. Also, results obtained from lower-fidelity models may be reused in order to regularize the identification of higher-fidelity models.

Numerical results

We consider a tensile test on a specimen with notches, as represented in Figure 5.10. A uniaxial loading with magnitude F is applied on left and right edges; other edges are free. The considered loading is pseudo-cyclic, with successive loading-unloading-reloading steps such that F goes higher than the threshold loading F ⇤ where plasticity is effectively observed in the FEM simulation (theoretically it appears immediately due to the singularity). The associated evolution curve of longitudinal stress-strain components is shown in Figure 5.11, while the progressive apparition of cumulative plastic strain is shown in Figure 5.12.

In the following, we perform identification with the DIC-mCRE method using concurrent mathematical models.

The only a priori information which is introduced in the model class is that material parameters are homogeneous.

We also assume that the Poisson ratio is known.

Multiscale mCRE strategy 6.2.1 Implementation of MsFEM into mCRE

We assume the objective is to identify constitutive parameters of a material described at the microscale, with known microstructure. The idea is to integrate this knowledge in the basis functions which are used through MsFEM. The previously developed mCRE-based identification framework is employed with multiscale admissible fields approximated from the MsFEM basis. They read ûε (x)=N ε (x) Û (kinematically admissible field), and vε (x)=N ε (x) V

(statically admissible field), with N ε defined from multiscale MsFEM basis functions. The resulting multiscale mCRE then reads (in its discretized form):

when DIC-mCRE is used, and:

when mI-DIC is used. The stiffness matrix K ε is computed from MsFEM basis functions.

The classical algorithm is then used to minimize the mCRE functional and recover the identified value p sol of material parameters. It is based on the two following steps, at iteration n +1(and considering DIC-mCRE):

•

Step 1: Compute optimal admissible fields ( Û(n+1) , V(n+1) ) for given

•

Step 2: Update model parameters for fixed admissible fields

Error indicators and adaptive process

Here again, the value of the CRE term E ε 2 CRE is an indicator on the quality of the numerical model. It can be postprocessed in an effective manner in order to conduct adaptive modeling. Out of error on the mathematical model, discretization error comes from three sources in the MsFEM context (macro mesh T H , boundary layers related to oversampling, and local micro mesh T h ). Three parameters are then tuned in an adaptive procedure in order to improve the accuracy:

• the local size H K of the macro mesh T H ; if the set {H K } is not chosen correctly, the coarse-scale discretization error may dominate the overall approximation error.

• the size d K (minimum patch radius) of the computational domains S K used to solve the microscale problems; if the set {d K } is not chosen correctly, the artificial boundary conditions set on @S K may strongly dominate the overall approximation error.

• the local size h K of the fine mesh used to solve the microscale problems (in each S K ); if the set {h K } is not chosen correctly, a fine-scale discretization error may dominate the overall approximation error.

In order to derive an error estimator on the overall modeling error (with construction of fully equilibrated stress field from the approximate MsFEM solution σ( V)=K ε ε(N ε V)) as well as error indicators on each of the MsFEM error sources, we follow the strategy developed in [START_REF] Chamoin | A posteriori error estimation and adaptive strategy for the control of msfem computations[END_REF]. It is compatible with the offline/online paradigm of MsFEM, with fine-scale computations at the element level only, and it enables an adaptive algorithm to be driven in which the various MsFEM discretization parameters are automatically chosen.

In this strategy, full equilibration with the hybrid-flux method is first perform to recover and estimate ⌘ dis on the overall MsFEM discretization error. A specificity is that the equilibrated tractions fK constructed on the edges of macro elements are oscillating (Figure 6 We indicate that most of the computations which are required to get the indicators are again performed in the offline stage of MsFEM. Each iteration in the algorithm consists of three groups of computations: (i) microscale (possibly costly) computations that can be performed in the offline phase of the iteration; (ii) macroscale computations that can be performed in this offline phase; (iii) macroscale computations that need to be performed in the online phase of the iteration, but which are relatively cheap. The adaptive procedure is in practice associated with the following considerations:

• the algorithm is initialized with a regular mesh T (0) H , no oversampling (S K = K) and a rough fine mesh size (namely h K = H K );

• when modifying the parameters h K , two values are considered in the adaptive process: h K = "/5 and h K = "/20 (finest mesh size at the microscale);

• when modifying the parameters d K from their initial value d K = H K i.e. S K = K, the oversampling size is determined by adding layers of progressive thickness ", 2", 3", . . . around the element K in the adaptive process. Such an adaptation is not performed when h K = H K ;

• refining the mesh T H , i.e. modifying the parameters H K , is performed using a quadtree method with nested elements. This step requires to handle hanging nodes.

The higher-order multiscale technique [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF] based on composition rules is beneficially used throughout the adaptive process. It enables for independent computations without coming back to the fine scale offline computations (no additional costly computations). The technique is used in its p-refinement version to compute a fully equilibrated field at the element level, while it is used in its H-refinement version when adapting the macro mesh. It thus allows us to refine the coarse mesh without the need to solve new fine-scale problems. However, when modifying parameters h K or d K in the adaptive process, new fine-scale problems need to be solved in order to define the updated multiscale basis functions.

Numerical results

We investigate the performance of the proposed methodology on two synthetic examples. For each of them, the mI-DIC version of the identification approach is used. Artificial images are used by considering a random distribution of gray levels over a 200 ⇥ 200 regular grid (Figure 6.6), and applying (from the imwarp Matllab command) a simulated displacement obtained from an overkill solution. A 2% measurement noise is added to images.

Example 1: periodic distribution of the Young modulus

We first consider a problem defined in the unit square domain Ω =( 0 , 1) 2 . The loading consists in a uniform body loading f v d =( 1, 0) T , while homogeneous Dirichlet boundary conditions are applied on @Ω. The material is assumed to be isotropic linear elastic, with a heterogeneous distribution of the Young modulus of the form:

1 (2 + P cos(2⇡(x 0.5)/")) (2 + P cos(2⇡(y 0.5)/")) (6.14) where damage occurs. We mention that the use of mCRE with parameter fields was investigated [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], Bonnet and Aquino, 2015, Babaniyia et al., 2017], and that a multiscale approach with mCRE was designed for composites [Huang, 2016].

Also, sparse regularization techniques (with L1-norm) coming from compressed sensing could be beneficially used with mCRE, in order to give a priori on the sparsity of damage phenomena being identified [START_REF] Mascarenas | Compressed sensing techniques for detecting damage in structures[END_REF], Zhang and Xu, 2015, Nagarajaiah and Yang, 2017, Zhang et al., 2017, Guo et al., 2018].

A second prospect is the extension to 3D full-field measurements with a large amount of noisy experimental data coming from DVC technologies (several Gigabytes or Terabytes). In this context, and for a given mathematical model, the proposed mCRE framework could be employed for data pruning, going from big data to smart data by collecting and selecting right data through the numerical model. This would be in complement to other tools (e.g. [Passieux et al., 2015b]) already introduced to facilitate the numerical processing of large data.

Another prospect is the integration of the proposed strategy in the general framework of DDDAS [Darema, 2015],

with sequential data assimilation and dynamic updating on evolutionary systems (with evolving properties) within a global feedback loop. On potential application is the design of numerically-assisted experimental tests on materials, with controlled loading, or the tracking of structural damage for structural health monitoring (e.g. see [START_REF] Ding | Dynamic data driven fault diagnosis of wind turbine systems[END_REF], Allaire et al., 2012, Prudencio et al., 2015, Kapteyn et al., 2020]). For this purpose, Kalman filtering with prediction-correction algorithm could be coupled with mCRE as envisioned in [Alarcon, 2012,Marchand et al., 2016] in order to estimate recursively and in real-time the material state (from optimal admissible fields of mCRE) as well as parameters (from the Kalman filter algorithm) by processing data sequentially. In this procedure, the stateobservation formulation of the dynamical system is reconsidered by changing the metric of the observation space compared to the classical metric used with Kalman filtering, and information on modeling error is directly given by mCRE. This information could be further used for optimal and dynamical management of computing resources, searching a trade-off between model accuracy and complexity according to experimental information, in order to be compatible with the real-time constraint.

We mention that real-time sequential data assimilation on mechanical structures with DIC measurements was recently addressed [Rubio et al., 2019a] for structural integrity on damageable concrete structures (with propagation of a single crack). It used a full stochastic strategy (Bayesian inference coupled with transport map sampling) and reduced order modeling to update model parameters, make fast and accurate predictions on outputs of interest, and perform online control. In order to extend this framework to more complex systems (e.g. with damage at the structural scale), the coupling between mCRE and Kalman filtering seems to be a relevant option.

In addition, a goal-oriented version of mCRE (as investigated in [START_REF] Djatouti | Coupling a goal-oriented inverse method and proper generalized decomposition for fast and robust prediction of quantities of interest in building thermal problems[END_REF], Djatouti et al., 2021] for thermal building) could be used in order to focus the identification effort on outputs of interest which are related to control aspects of DDDAS, and therefore further reduce the computational cost.

Appendices

Appendix A

List of personal publications

National & international conferences