Toward 4D Mechanical Correlation

Background: The goal of the present study is to illustrate the full integration of sensor and imaging data into numerical procedures for the purpose of identiﬁcation of constitutive laws and their validation. The feasibility of such approaches is proven in the context of in situ tests monitored by tomography. The bridging tool consists of spatiotemporal (i.e., 4D) analyses with dedicated (integrated) correlation algorithms. Methods: A tensile test on nodular graphite cast iron sample is performed within a lab tomograph. The reconstructed volumes are registered via integrated digital volume correlation (DVC) that incorporates a ﬁnite element modeling of the test, thereby performing a mechanical integration in 4D registration of a series of 3D images. In the present case a non-intrusive procedure is developed in which the 4D sensitivity ﬁelds are obtained with a commercial ﬁnite element code, allowing for a large versatility in meshing and incorporation of complex constitutive laws. Convergence studies can thus be performed in which the quality of the discretization is controlled both for the simulation and the registration. Results: Incremental DVC analyses are carried out with the scans acquired during the in situ mechanical test. For DVC, the mesh size results from a compromise between measurement uncertainties and its spatial resolution. Conversely, a numerically good mesh may reveal too ﬁne for the considered material microstructure. With the integrated framework proposed herein, 4D registrations can be performed and missing boundary conditions of the reference state as well as mechanical parameters of an elastoplastic constitutive law are determined in fair condition both for DVC and simulation. regularization (cid:7) m voxels. The corresponding standard displacement uncertainty is again equal to voxel. When the T4-DVC code is run, the standard displacement uncertainty is equal to 0.07 voxel when (cid:7) = 37 voxels. A gain of about a factor 3 is observed in the present case, which is close to the a priori estimate (i.e., 2 √ 3 ≈ 3 . obtained from the uncertainty analysis (see “Mesh of the ROI” section).


Background
The emergence of simulation-based engineering sciences calls, among many challenges [1], for robust validation procedures and uncertainty quantifications to achieve reliable predictions. Full-field measurements are one way of bridging experimental and computational mechanics. Their advantage lies in the fact that the comparison is now achieved by using huge amounts of data, including 3D imaging, to probe the predictive capacity of material models [2] and numerical frameworks [3]. The aim of the paper is to show that a seamless procedure, hereafter called "integrated 4D registration," can be formulated to analyze an in-situ test performed in a lab tomograph for the purpose of identifying a nonlinear constitutive law and unknown boundary conditions. Computed microtomography allows 3D images of materials to be obtained, which reveal the microstructure in the bulk in a non-destructive way [4][5][6][7][8]. Very early on, mechanical tests were performed in-situ [9][10][11][12]. For instance, the damage development in particulate composites could be analyzed [11,13,14]. Creep has also been studied with such techniques [15][16][17][18]. One additional feature is to quantify kinematic fields via digital volume correlation (or DVC [19,20]). This additional piece of information can be used to analyze fatigue crack propagation [21,22] or to validate finite element simulations of cracked samples [3].
Very few studies deal with the identification of material parameters based on volumetric analyses. A first reason is related to the measurement uncertainties and biases that are usually higher than those typically encountered in 2D analyses [23,24]. Second, the computational environment needed to solve these inverse problems is generally very intrusive (e.g., constitutive equation gap method [25][26][27], virtual fields method for nonlinear constitutive laws [28], equilibrium gap method [29]). Third, nonlinear constitutive equations require spatiotemporal (i.e., 4D) analyses to be considered, which are both experimentally and computationally very demanding. The aim of the present work is to show the feasibility of such a framework to study the elastoplastic behavior of spheroidal graphite (SG) cast iron.
One of the goals of the present study is to achieve full integration of sensor and measurement data in the developed numerical procedures [1]. This type of analysis corresponds to so-called integrated approaches [30,31], which were developed up to now mostly with 2D images of sample surfaces [32][33][34][35][36] in which material parameters are directly measured from image registrations. A first extension to 3D images was recently proposed to analyze an indentation experiment on plasterboard monitored via tomography and DVC [37]. When material parameters are sought, the corresponding kinematic and static sensitivities [38] are needed. Finite element simulations are used in particular to generate a set of kinematic fields that are further used for DVC purposes. It was shown that non-intrusive spatiotemporal schemes can be used for 2D images [39], whereby the incorporation of time leads to 3D integration. In a similar spirit but extended to 3D (spatial) images, the following analysis corresponds to a first step toward 4D integration.
In "Methods" section, the experimental configuration is presented. DVC analyses based upon discretizations with 4-noded tetrahedra (i.e., T4-DVC) are used to measure kinematic fields without any mechanical regularization. From the acquired scan in the reference configuration, a mesh is adapted to the actual shape of the sample. From the kinematic measurements, the measured nodal displacements of the top and bottom faces of the region of interest (ROI) become the boundary conditions of the integrated scheme in which the material parameters are sought. A non-intrusive setting is developed. The previous framework is applied to analyze the described tensile test on spheroidal graphite cast iron to determine unknown boundary conditions, and to identify elastoplastic parameters ("Results and discussion" section). Uncertainty quantifications are performed in addition to convergence analyses in which the mesh quality is assessed as well.

In situ mechanical test
A dog-bone sample is obtained via electrodischarge machining (EDM) from a 1.6-mm thick plate made of SG cast iron. To ensure that the specimen will break in the ligament area and not in the grips, the central part is thinned with a radius of 20 mm (Fig. 1a). This zone is 9-mm high, while the smallest cross-section of the sample is 1.6 × 1 mm 2 . The tensile test reported herein has been carried out in the X50+ tomograph of LMT. The specimen is loaded in nine incremental loading and unloading steps (Fig. 1b) in a custom made testing machine, which is placed on the turntable of the tomograph [40]. It is worth noting that the reference scan is acquired at a load level equal to 165 N. The first two scans of the reference configuration allow displacement uncertainties to be evaluated.
A complete scan of the sample corresponds to a 360 • rotation along its vertical axis, during which 1000 radiographs are acquired with a definition of 1280×1860 pixels. Each scan lasted less than 40 min. The physical voxel size is 6.4 µm, and the reconstructed volume is encoded with 8-bit deep gray levels. This test was analyzed via global approaches to DVC with regular meshes made of 8-noded cubes to reveal the damage micromechanism (i.e., nodule debonding) with the correlation residuals [41]. In the following, the underlying macroscopic behavior will be studied.

Digital volume correlation
DVC is used to measure 3D displacement fields in the bulk of samples [42,43]. In its incremental version, DVC consists of registering two volumes by evaluating displacement fields that yield the best possible match. In the early developments, small zones of interest or ZOIs (i.e., small interrogation volumes [19,20,44,45]) in the considered region of interest (ROI) are registered. This type of approach is nowadays referred to as local (i.e., the only information that is kept is the mean displacement assigned to each analyzed ZOI center) and incremental since only two volumes are considered. Each registration is performed over the ZOI with various kinematic hypotheses, possibly accounting for its warping. This type of analysis is not considered herein since it does not offer an interpretation of the kinematic measurement that is "congruent" with that of FE analyses.
Global and incremental approaches, which appeared more recently [46], consist of performing a registration over the whole ROI. The kinematic bases are then defined over the whole ROI (i.e., Rayleigh-Ritz formulations) or over a discretized ROI (i.e., Galerkin formulations). The fact that continuity is enforced allows the measurement uncertainty to be reduced in comparison with the same discretization in which the nodes are not shared [47]. Another important advantage of global DVC in comparison with local DVC is its direct link with numerical simulations of mechanical problems [3,21]. When dealing with digital volumes it is natural to consider regular meshes made of 8-noded cubes (i.e., C8-DVC [46]) even though unstructured meshes may also be considered [48]. In the following, unstructured meshes made of 4-noded tetrahedra are chosen to allow for a more faithful description of external surfaces ("Mesh of the ROI" section) in comparison with structured C8 meshes.
The incremental registration problem (i.e., digital volume correlation) consists in minimizing the sum of squared differences over the considered ROI for the reconstructed volumes in the reference configuration I 0 (x) and deformed configuration I t (x). γ 2 I denotes the variance associated with acquisition and reconstruction noise of the gray level volumes. For the sake of simplicity and in line with the chosen formulation (1), the latter is assumed to be white and Gaussian. This is a rather crude assumption for computed tomography, for which as its very name indicates, the raw data (i.e., radiographs) are processed to provide the reconstructed volume [49], and thereby the original noise in the radiograph displays in the reconstructed volume a nonhomogeneous character and exhibits strong and anisotropic spatial correlations. However, even if this assumption is not optimal, from previous experience, it revealed sufficient for handling DVC for standard tomographic images. Correlation residuals however do reveal systematically more pronounced levels close to the rotation axis and artifacts such as rings [23,50].
Equation (1) expresses the underlying assumption of gray level conservation for each considered voxel x of the ROI. This least squares problem is nonlinear. It is solved with an iterative Gauss-Newton scheme [51] in which the stationarity condition is recast in the following variational formulation [52].
and where ρ(x) = I 0 (x) − I t (x + u(x)) is the current gray level residual, u the current estimate of the sought displacement field, δu the correction to the current displacement field, and δv an increment to the trial displacement field.
The subspace to which the measured and trial displacements belong is chosen to be the vector space generated by a set of fields ϕ n (x) so that where υ n are the sought amplitudes at the considered time t (υ n (0) = 0), which are gathered in the column vector {υ}. The minimization procedure consists of successive linearizations and corrections {δυ} to the measured degrees of freedom [53] [M]{δυ} = { b} (4) with where { ρ} gathers all voxel-scale dimensionless residuals ρ(x)/ √ 2γ I , and † denotes the transpose of a matrix or column vector. The DVC matrix [M] and the corresponding residual vector {b} are constructed with the matrix [ ] that gathers all scalar products of the basis fields by the volume contrast (i.e., (x, n) = ϕ n (x) · ∇I 0 (x)/ √ 2γ I ). The covariance matrix [Cov υ ] of the measured degrees of freedom is given by [54] [ In the following analyses, 4-noded tetrahedra will be considered. Consequently the kinematic basis made of ϕ n (x) fields corresponds to the shape functions of T4 elements. Such registration approach will be referred to as incremental T4-DVC.

Mesh of the ROI
To perform T4-DVC analyses, the ROI of the sample needs to be meshed based on its true geometry. This type of procedure usually requires the acquired scan to be processed with mathematical morphology and ad hoc procedures [55][56][57]. When multiphase microstructures are studied, this operation becomes more complex to make sure that the mesh is compatible with the morphology of the phases. For instance, adaptive level-set method can be used [58]. In the present case, the simulations will not be performed at the microstructural scale but using a macroscopic description of the material behavior. Consequently, the mesh will not be adapted to the underlying microstructure. However, it needs to comply with the external surface of the sample.
A digital image correlation procedure is implemented to adjust a regular mesh made of 4noded quadrilaterals to the actual surface of the sample [59]. Different transverse sections are considered. From these registrations, a 3D discretization made of 8-noded hexahedra (H8) is obtained. It is subsequently converted into a T4 discretization by subdividing each H8 element into six T4 elements of equal volume. The two meshes are exactly nodeequivalent, i.e., the number of nodes and the node positions are unchanged.
With a given discretization, one critical issue in DVC analyses is the measurement uncertainty. Because the registration is an inverse problem, it is expected that finer meshes will lead to poorer measurement uncertainties [47]. This effect can be understood by estimating the uncertainty associated with acquisition noise. Equation (6) shows that the covariance of measured degrees of freedom is related to the inverse of the DVC matrix. The latter is assembled by using elementary matrices [M e ] of the T4 elements used in the mesh. The latter contains scalar products of the volume contrast (i.e., ∇I 0 (x)) and the vector shape functions (i.e., N i (x)e p , with i = 1, 4 and p = 1, 3 where VOI e corresponds to the integration domain of the considered T4 element, N i the (linear) shape functions of the T4 element, and δ pq Kronecker's delta. If the correlation length of the gray level volume is smaller than a characteristic length scale of the element a mean field approximation can be made [47] (e p · ∇I 0 )(e q · ∇I 0 ) ≈ 2 where μ is the physical size of one voxel. The standard uncertainty for each degree of freedom reads Equation (11) shows that the measurement uncertainty is related to the noise to signal ratio γ I / I within the considered ROI, and is inversely proportional to the square root of the number of voxels within the considered T4 element. The standard displacement uncertainty γ υ e is two times lower than that of a regular hexahedron with the same volume [47]. The elementary matrices [M e ] are then assembled to evaluate the global DVC matrix [M] and the corresponding covariance matrix [Cov υ ]. In the assembly process, the nodal displacements will now be evaluated over more than one element (i.e., all the elements sharing the considered node). Consequently, the number of voxels 3 g considered for the evaluation of the nodal displacements will become larger than that of each individual element. For a regular mesh made of C8 elements, the corresponding volume is equal to (2 ) 3 [47]. With the chosen procedure to transform H8 meshes into T4 meshes, inner T4 nodes are shared by 24 elements. This discussion shows that when using a T4 mesh with all its connectivities enforced, a reduction factor of the order of 2 √ 6 is to be expected on the standard displacement uncertainty in comparison with the same discretization in which each T4 element is considered independently. An additional gain is to be expected thanks to the overall continuity requirement of the displacement field [47,54]. Last, if the same element sizes as defined herein are considered for T4 and C8 meshes, the measurement uncertainties of the former are expected to be 2 √ 3 times lower than the latter. The present results show that the finer the mesh (i.e., the smaller ), the higher the measurement uncertainties. From a purely metrological point of view, finer meshes are to be discarded since they will lead to very high measurement uncertainties or even to singular DVC matrices that will not lead to trustworthy results, if any is obtained. However, on the mechanical side of the problem, coarse meshes may degrade the quality of the simulations.
In the context of the Finite Element Method, verification procedures have been introduced since the late 1970s [60][61][62]. They enable a posteriori discretization error estimates to be computed and mesh adaptivity to be driven. In this work, a verification procedure based on the concept of constitutive relation error (CRE) [61] is implemented. It leads to robust error estimates for both linear and nonlinear time-dependent problems [63][64][65][66]. The CRE concept, which is particularly suited to Computational Mechanics models, is based on duality and rests on a simple idea, namely, after constructing so-called admissible fields satisfying all equations of the model except (part of) the constitutive law, the residual associated with the constitutive relationship is evaluated.
The constitutive equation investigated herein is Ludwik's law [67], which is written within the continuum thermodynamics framework [68]. The corresponding state potential reads [69] where is the total strain tensor, p the plastic strain tensor, p the cumulated plastic strain, C Hooke's tensor, and (K , n) the hardening parameters. The state laws are derived from the state potential The intrinsic dissipated power density D becomes The pseudo potential of dissipation ϕ * (σ, R) is the indicator I f of the elastic domain. The yield function is written as where J 2 is the second invariant of the deviatoric stress tensor, and σ y the yield stress. The growth laws becomė whereλ is the plastic multiplier that is obtained from the consistency conditionḟ = 0. Let us introduce the dual pseudo potential of dissipation ϕ(˙ p , −ṗ) defined in the Legendre-Fenchel sense as Within this standard material formulation with internal variables, the dissipation error is the relevant tool in the CRE concept [65]. It consists of dividing constitutive equations into two groups: • constitutive equations related to the free energy (i.e., balance equations, kinematic constraints, initial conditions, state laws) that define the concept of admissibility; • constitutive equations related to the dissipation (i.e., growth laws).
For a given admissible solution (˙ p ,ṗ,σ,R) satisfying the first group of equations, the local dissipation error functional is constructed as and the global dissipation error, θ CRE , which is used as a posteriori error estimate, is obtained from Performing the appropriate change of variables (p, R) → (p,R) so that a normal formulation with linear state lawR = γp is obtained [63], it can be shown that the error estimate (19) is directly related to a gap between exact and FE solutions [65]. This key result is analogous to the Prager-Synge theorem for elasticity problems [70]. Moreover, a relative error is defined as with The construction of an admissible solution (˙ p ,ṗ,σ,R) is performed by post-processing the computed FE solution [63,65]. The main technical part is the construction of an admissible stress field that satisfies balance equations in the strong sense. The technique used herein is based on a domain decomposition approach [61,65,71].
The uncertainty analysis and the verification steps are not necessarily compatible (i.e., the discretization used for DVC purposes requires not too fine meshes to be considered). Therefore a compromise needs to be found if the two approaches are performed sequentially, otherwise the two meshes will be different and the interpolation quality from one mesh to another will not be checked against experimental data. In the following, it will be shown that integrated approaches will no longer require this trade-off to be implemented since both steps will be performed in a unique calculation in which the mesh is no longer a limitation since the number of unknowns will be drastically reduced.

4D mechanical correlation
Up to now, the implicit regularization of the correlation procedure is related to the definition of ZOIs (in a local approach) or elements (in a global approach) and the displacement interpolations therein since the inversion problem cannot be solved at the voxel scale. Incremental FE-based DVC can be seen as a strong regularization of local incremental DVC.
An additional way of regularizing correlation procedures is to consider mechanical admissibility. This is, for instance, possible by requiring that the measured displacement field minimize a weighted sum of the DVC objective functional and the equilibrium gap by assuming, say, an elastic behavior [72]. These so-called mechanics-aided (i.e., regularized) DVC procedures enable the high frequency kinematic components that are not mechanically admissible to be filtered out [47]. This penalization acts as mechanical filter in registration procedures. Because mechanical admissibility controls the high spatial frequencies, fine meshes (that cannot be considered in standard DVC calculations) become accessible [41,72]. They also allow registrations to be performed at the voxel level [73]. In this extreme situation, implementations on Graphics Processing Units (GPUs) enable very large computations to be run [47]. One alternative route consists in considering PGD-based DVC [74].
Another regularization route consists in performing spatiotemporal registrations [36,75,76]. This type of approach will be referred to as 4D correlation. The spatiotemporal registration problem aims to minimize the sum of squared differences over space and time The displacement field is expressed in a spatiotemporal way as where φ m (t) are the temporal discretization functions, and υ mn the spatiotemporal kinematic unknowns. A final step consists in requiring the measured displacement fields to be fully mechanically admissible (i.e., they satisfy equilibrium for a chosen constitutive law). The first 3D developments were based on elastic simulations for which a non-intrusive framework was developed [37]. The finite element code is then used to generate kinematic fields that are needed for DVC purposes. When material parameters are sought, the corresponding spatiotemporal sensitivities [38] are also needed for nonlinear constitutive postulates. The 4D kinematics is thus parameterized with the sought corrections {δp} to the current material parameters { p} where [∂u/∂{p}] is the matrix gathering all the spatiotemporal sensitivity fields to the sought material parameters. In the present case, the number of acquired scans is small (i.e., 9). Consequently, the time regularization described by the temporal functions φ m (t) will not be considered since loading / unloading sequences are performed and very few scans are available for each loading / unloading step (see Fig. 1b). Conversely, the fact that the sensitivity fields are available still regularizes very strongly the 4D (mechanical) correlation procedure. Thus the discretized displacement field becomes in which the kinematic degrees of freedom υ n are not independent but all linked via their sensitivities to the sought parameters The kinematic sensitivities associated with a chosen spatial discretization are gathered in matrix [S υ (t)] (i.e., (S υ (t)) nm = (∂υ n /∂p m )(t)] that is evaluated for each time step t. The 4D mechanical correlation then consists of minimizing the sum of squared differences directly with respect to {p} As in incremental DVC, a Gauss-Newton scheme is implemented for which the corrections {δp} to the sought material parameters satisfy at a current iteration with where is the instantaneous weighted kinematic Hessian [39]. The covariance matrix [Cov p ] of the measured parameters reads with The second term of Eq. (30) is due to all the correlations associated with the volume in the reference configuration. If the cross-correlations are neglected, Eq. (30) reduces to [39] [ In the present case, load measurements gathered in vector {F meas } are also available. They can therefore be compared with the FE predictions in which the reaction force vector {F FE } is due to the fact that measured displacements are prescribed on the top and bottom boundaries of the considered ROI. These Dirichlet boundary conditions are measured with incremental T4-DVC and applied to the FE model. The global equilibrium gap over the whole loading history reads where γ 2 F denotes the variance of the measurement uncertainty of load cells. The latter is assumed to be white and Gaussian.
The minimization of ρ 2 F alone corresponds to updating the FE model by considering only global equilibrium to estimate the sought parameters {p}. It is referred to as loadbased finite element model updating (or FEMU-F [77,78]). As for displacement-based approaches, the load sensitivity matrix [S F (t)] to the sought parameters is computed to update {δp} from the current estimate {F FE ({ p})} of the reaction forces. A Gauss-Newton scheme can also be implemented so that the corrections {δp} to the sought material parameters satisfy at a current iteration /γ 2 F the instantaneous static Hessian [39], and The covariance matrix [Cov p ] of the measured parameters is given by In the following, both sets of data, namely, reconstructed volumes and measured loads will be considered in a single approach. Since the physical nature of the data (i.e., gray levels and force) are different, special care has to be exercised. A Bayesian framework is considered hereafter [39,59] in which each considered information is weighted by its merit (i.e., variance and covariance with all the other data). 4D mechanical correlation finally consists of minimizing χ 2 with respect to the sought parameters with and where n F denotes the number of load measurements per scan (i.e., n F = 1 in the present case). The global system to solve for each Gauss-Newton iteration reads and The covariance matrix [Cov p ] of the measured parameters is given by [39] [ whose approximation when neglecting correlations associated with the reference volume reduces to The integrated DVC code, Correli 3.0 [79], is a Matlab implementation that uses in a non-intrusive way a finite element code to compute the spatiotemporal sensitivity fields (Fig. 2a) in addition to the current estimates of the displacement fields and reaction forces. The inputs to such simulations are the measured displacements of some boundaries of the ROI. In the present case, the commercial finite element package Abaqus/Standard is used with its built-in constitutive laws. C++ kernels compute all the other data needed to perform DVC analyses, namely, the DVC matrix [M], the instantaneous voxel residuals ρ(x, t), and the instantaneous nodal residual vector {b(t)}. Binary MEX files are generated and called in the Matlab environment to compute the Hessians and residual vectors to be utilized in the Gauss-Newton schemes introduced above.
Since measured boundary conditions drive the finite element simulations, a first T4-DVC analysis is required. It uses the same MEX files as I-DVC (Fig. 2b). However, there is no need for any mechanical finite element code since no mechanical regularization is performed herein via T4-DVC. If the mesh used in I-DVC is finer than that used in T4-DVC, the measured displacement fields are interpolated by using the T4 shape functions.

Results and discussion
In situ mechanical test Figure 3 shows the reconstructed volume in the reference configuration (i.e., F = 165 N). The coarse microstructure of cast iron can be observed in which the graphite nodules appear dark and the ferritic-pearlitic matrix in light gray levels. This is due to the fact that the X-ray absorption of iron is significantly higher than that of carbon. The mean volume fraction of nodules is equal to 15 %.

Mesh of the ROI
In the reconstructed 3D images, an ROI is defined in which the displacement field is to be measured. The size of the ROI is adjusted to the sample geometry within a volume size of 500 × 420 × 1215 voxels. Figure 4a shows the H8 mesh, which is then transformed into a T4 mesh (Fig. 4b) that will be considered in the results reported herein. The mesh density was chosen to allow T4-DVC to converge. It is the finest discretization that could be considered because of the rather coarse microstructure of cast iron (Fig. 3). This density will be labeled as 1 (Table 1).  Finer meshes will also be considered in integrated analyses (Table 1). For each mesh, two different equivalent lengths are given. The first one corresponds to the cubic root of the mean volume of the T4 elements. This quantity was introduced for the evaluation of the measurement uncertainty when each T4 element is considered independently. However, to define the spatial resolution associated with displacement measurements of global DVC a second length is also evaluated. It corresponds to the cubic root of the mean number of voxels considered for the nodal displacement measurement. This length g defines the spatial resolution of the measurement technique. Table 1 shows that g is less than three times the element size .

Measured displacement fields
Before studying the displacement field between two scans, it is important to evaluate the level of uncertainty attached to the natural texture of the 3D images. Since the studied SG cast iron has a coarse microstructure it is necessary to make a compromise between the spatial resolution and displacement uncertainties [41]. For C8-DVC it is found that the smallest element size is 32 voxels to reach convergence. The standard displacement uncertainty is equal to 0.2 voxel. With regularized C8-DVC [47,72], 16-voxel elements can be considered with a regularization length m = 16 voxels. The corresponding standard displacement uncertainty is again equal to 0.2 voxel. When the T4-DVC code is run, the standard displacement uncertainty is equal to 0.07 voxel when = 37 voxels. A gain of about a factor 3 is observed in the present case, which is close to the a priori estimate (i.e., 2 √ 3 ≈ 3.5) obtained from the uncertainty analysis (see "Mesh of the ROI" section).
The root mean square (RMS) residual is equal to 6.5 gray levels, which is an estimate of acquisition and reconstruction noise (i.e., √ 2γ I ). The T4-DVC code is then run incrementally to estimate the total displacements for the 9 analyzed steps. Figure 5 shows the change of the dimensionless RMS gray level residual ρ 2 ROI /2γ 2 I for all analyzed scans. These levels are virtually constant for the whole analysis except the first one where no motion has occurred between the two acquisitions. Furthermore the dimensionless residuals remain very close to 1 (i.e., an increase of 30 % at most is obtained). This trend was also observed in C8-DVC and regularized C8-DVC (i.e., RC8-DVC) analyses [41].
Longitudinal displacement fields are reported in Fig. 6 for three different load levels. As the applied load level increases so do the displacement amplitudes. It is worth noting that even for the last load level no localization is observed. The material is presumably still in the hardening regime of plasticity.

Integrated DVC
In the sequel all reported results will consider displacement fields that are mechanically admissible in a finite element sense. Consequently, the first question to address is the  quality of the finite element simulations. The next issues will be related to the analysis of the experiment per se. It is worth remembering that when finite element simulations are carried out, the material model and the parameters need to be known for cast iron. Two sets of parameters (see Table 2) will be considered for Ludwik's law ("Mesh of the ROI" section). The first one was obtained via weighted FEMU-UF [41] to analyze a cyclic tensile test on the same material at the macroscale. The second one was determined from the analysis of a standard monotonic tensile test [40].

Verification
Let us first analyze the discretization error level with respect to the FE mesh. Five uniform meshes with different densities are considered. The relative errors θ are reported in Table 1. First, measured (Dirichlet) boundary conditions are directly applied to the FE models. As expected, the discretization error decreases with the mesh size [65,80]. However, the mean error remains very high. To understand this result, the associated local errors in space θ 2 E = T 0 VOI e η dxdt are shown for each element e of the mesh in Fig. 7. Irrespective of the mesh density, most of the discretization errors concentrate in the vicinity of the two boundaries where the (measured) boundary conditions are applied.
To evaluate the effect of boundary conditions, the measured fields are fitted with low order polynomials on both surfaces where they are prescribed. Several interpolations, which are referred to as Q#, are defined by the number of monomials picked in the list {1, x, y, xy, x 2 , y 2 , x 2 y, xy 2 }, where x and y denote the two spatial coordinates over which the interpolation is performed. Therefore, Q1 corresponds to constant interpolation, Q3 to linear interpolation, and Q6 to quadratic interpolation. For a linear interpolation, there is a clear gain in terms of discretization errors when compared to the initial results (Fig. 7). Table 3 shows that the RMS interpolation error between the measured displacements and their interpolations is very high for the uniform case (i.e., more than 36 times the standard displacement uncertainty). This is due to the fact that rigid body rotations are not accounted for. The interpolation errors for the other descriptions are close and significantly lower than in the previous case (i.e., of the order of 3 times the standard displacement uncertainty). For a given mesh density (i.e., 1), the associated quality θ is also given in Table 3. There is a clear influence of the interpolation on the mesh quality since a fivefold gain is observed when the measured boundary conditions are linearly interpolated. This is also confirmed when the local contributions to the CRE error are compared (see Fig. 8). The effect on the CRE error is even more pronounced when finer meshes are considered, namely, up to an eightfold gain is observed for mesh densities ranging from 3 to 5 (Table 1). This result proves that most of the CRE errors are due to Dirichlet boundary conditions constructed with the measured displacement fields whose high frequency fluctuations are due to noise.

Boundary conditions
Since the reference scan was acquired for a load level equal to 165 N (Fig. 1), the zero load configuration is not known. This lack of information is of no consequence when standard DVC analyses are run (i.e., C8-DVC and T4-DVC) or even regularized DVC. Conversely, when integrated approaches are performed in which the constitutive law is nonlinear, the zero load configuration has to be either known or estimated. In the present case, the second route has to be followed. Since the measured displacement fields u meas (x, t − t 0 ) have been evaluated with respect to a nonzero load configuration (at time t 0 for which the load level is equal to 165 N), an additional field u(x, t 0 ) has to be added to get the experimental displacement from which the Dirichlet boundary conditions are extracted to drive FE simulations This additional field being unknown, it will be evaluated via FEMU-F. Since the first increment is elastic, it is assumed that where t 1 corresponds to the 300-N load level, and α an unknown amplitude to be determined. For a mesh density of 1 and the first set of material parameters, Table 4 shows the effect of different values of the amplitude α on the static residual χ F . There is a clear influence of α on χ F , and more importantly, the standard uncertainty γ α is very low in comparison with the reported level of α. The fact that χ F is significantly larger than unity is an indication of model errors. This first analysis allows the load residuals to be decreased from 26 (when α = 0) to 18 when α is optimized via the load residuals, which corresponds to a 30 % decrease.
The next question to address is whether a more complex correction can be proposed. Instead of having the same amplitude on the longitudinal and transverse components, it is proposed to look for three different amplitudes α x , α y and α z for each direction. Before performing the identification, a sensitivity analysis is carried out [81]. It consists in computing the dimensionless static Hessian for a given variation of the sought parameters (here chosen equal to 10 %) where [δF] gathers all reaction force increments associated with the parameter changes. Any linear combination of parameters leading to an eigen value of [H F ] greater than or equal to 1 will more sensitive than the noise level. The load cell of the testing machine used herein has a standard uncertainty equal to 3.5 N. In the present case, the dimensionless static Hessian reads and has only one eigen value greater than 1. The corresponding eigen vector is essentially aligned with the third parameter direction (i.e., α z ). The other two eigen values are at least three orders of magnitude lower. The new parameterization is therefore not sensitive enough with the data at hand. Similarly, when α x and α y are assumed to be identical, the largest eigen value is greater than 1 and its eigen direction is again aligned with the direction of α z . The second eigen value is three orders of magnitude lower. This last parameterization is not compatible with the available information. Consequently only the first one is considered hereafter.

Effect of mesh density
When using integrated approaches, the limitation associated with the uncertainty / spatial resolution no longer applies [47]. Five different mesh densities are therefore investigated (Table 1). I-DVC is run for a given set of material parameters and with the same α value.
Since two different types of data are combined, two residuals are first analyzed, namely, the static residual χ F and the registration residual χ I . Figure 9 shows the change of both residuals as functions of the mesh density. There is a slight increase of the load residual and a small decrease of the correlation residual with the mesh density. Consequently, the effect of the mesh density is minimal in the present case. Further, since measured boundary conditions are prescribed and the sample geometry remains rather simple explains that the levels of χ I are significantly lower than those of χ F . Last, if the combined residual χ is computed, it is virtually identical to χ I since the number of voxels in the ROI (more than 45 × 10 6 ) is very large in comparison with the load measurement per scan (i.e., 1). If Bayesian weighting is followed, it is concluded that mesh density 5 leads to the best results. Since T4-DVC results are also available for the coarsest mesh, the correlation residuals are compared to those of I-DVC with different mesh densities. In the present case, the history of instantaneous dimensionless registration residuals is again analyzed (Fig. 5). Contrary to T4-DVC it is observed that the registration quality degrades as the scan number increases, i.e., more yielding has occurred. The fact that this trend is identical for all five mesh densities is an indication of a material model error and not a discretization error. To obtain the results discussed in this section the computation time ranges from a few minutes for mesh density 1 to two hours for the finest mesh on a workstation (with 2.6-GHz Intel Xeon E5-2650 processor). Let us emphasize that this computation time remains small as compared to the duration of the experiment, and even smaller if preparation time is considered.

Comparison of two sets of material parameters
Having based 4D registrations on a single set of material parameters, the next question to address is their influence on the results. A second set of parameters (see Table 2) is used in the same analysis as in the previous sections. Figure 10 shows the values of the unknown loading parameter α as functions of the mesh densities for the two sets of material parameters when using FEMU-F. This convergence analysis shows that for the two sets of material parameters, the mesh densities 4 and 5 lead to approximately the same value of α and χ F . There is a small influence (i.e., less than 3 %) of the set of material parameters on the value for the finest meshes. This is to be expected since the first load level is assumed to be elastic and the corresponding Young's moduli are very close. Further, the standard uncertainty γ α decreases by about 15 % from density 1 to 5. Table 5 shows the load and registration residuals for the two sets of material parameters. For the load residual, the second set of material parameters is better than the first one. There is a clear sensitivity of χ F to the two sets. It is worth noting that the finer the meshes the sightly higher the load residuals in both cases. For the registration residuals, the first set of material parameters is better than the second one. This is also true for the global residual χ. These results indicate that even though the number of scans is very small (i.e., less than ten), the two residuals are sensitive to the material parameters.
In both cases, it is worth noting that the load, registration and total residuals are virtually unchanged for the last two densities. All of them have converged so that density 4 can be considered as the finest mesh needed for identification purposes.

Identification of material parameters
In the following, only the plastic parameters (i.e., σ y , K and n) are identified. By performing the same sensitivity analysis as in "Boundary conditions" section, it is concluded that two eigen values of [H F ] are greater than 1 for a parameter variation of 10 %. Instead of setting the least sensitive parameter to its initial guess, a Tikhonov-type of regularization [82] is preferred. It consists in adding a penalization term from the initial guess {p 0 } of the sought parameters, say in the FEMU-F procedure, where λ F is set in such a way that at the end of the relaxation procedure is has the same order of magnitude as the second eigen value of [H F ]. The relaxation procedure [59] consists in iterating the present scheme by starting with λ F equal to one thousandth of the maximum eigen value of [H F ]. Once convergence is reached (i.e., each relative parameter variation is less than 10 −5 ), the regularization length λ F is divided by ten and the procedure is repeated until λ F has the same order of magnitude as the second eigen value. Table 6 shows the results of the identification via FEMU-F. The values of the optimal parameters are virtually identical for the four mesh densities. In the present case, the quality of the mesh has only a minor influence on the identification results. The simple geometry of the studied sample may explain such results. As expected from the sensitivity analysis, the maximum uncertainties are of the order of 10 % at most. Two out of three parameters have varied significantly from the initial guess (i.e., from the first set of parameters) whereas the hardening modulus has remained virtually unchanged. The static residual has been lowered to 14.2, which is still a very high value in comparison with  Effect of the interpolation of the Dirichlet boundary conditions on the identified material parameters (mesh density: 1). The reference set of parameters is that when no interpolation is applied to the measured boundary conditions acquisition noise. However, it is lower than what was observed for the two sets of material parameters (Table 5). Figure 11 shows an opposite dependence of the two identification residuals with the mesh density. However, the relative variations of χ F are smaller than those observed for the two sets of parameters (Table 5). Conversely, χ I is slightly more sensitive to the mesh density.
In the following, 4D registration is coupled with load measurements to identify material parameters. If a Bayesian framework is followed, the weight associated with load data is very small and the results are identical to those in which no load data is accounted for. As a consequence, the load residual reaches very high values (i.e., more than 100), which is not acceptable. It was therefore chosen to give equal weight to all the static and kinematic data by minimizing χ 2 = 1/2(χ 2 I + χ 2 F ) instead [39]. Table 3 shows that the residuals at convergence for mesh density 1 are close to those observed in FEMU-F for  Table 6.

Sensitivity to boundary conditions
The verification procedure has shown that the errors are located close to the boundaries where the measured displacements are prescribed. To evaluate the effect of the prescribed displacement interpolation on the identification results, a first analysis is performed for the coarsest mesh density. Since the boundary conditions have been interpolated, they may induce changes on the identified material parameters. Figure 12 shows the relative changes of the three identified material parameters for the different interpolations. The reference set of parameters is that when no interpolation is performed. Except for the uniform interpolation, the total residual varies less than 1 %. The registration and load residuals are also similar for all analyzed cases but the uniform interpolation ( Table 3). The latter is clearly not acceptable in terms of registration residuals. From this analysis, the linear interpolation (Q3) is the best compromise between mesh quality and identification results when mesh density 1 is studied.
In the following, only the linear interpolation is studied for the different mesh densities. When the load and registration residuals are computed for all these cases with the same material parameters (identified when the measured boundary conditions are considered), the registration residuals are identical and the load residuals slightly increase with the mesh density as already observed (see Fig. 11). These results indicate that the identification procedure does not need such fine meshes. One of the reasons of this result is likely to be related to the fact that the load residuals are still very high (i.e., there is a significant model error). Figure 13 summarizes the identification results for the finest mesh (i.e., leading to the smallest CRE error). The gray level residuals are very low in comparison with the reference microstructure (see Fig. 3b). This is an additional proof that the registration is successful not only globally but also locally. When the measured load history is compared to the corresponding predictions, higher differences are observed. This result is to be expected from the global residuals, which are significantly higher than the measurement uncertainties.

Conclusions
It has been shown herein that full integration of sensor and measurement data can be achieved in numerical procedures via so-called 4D mechanical correlation. In the present case, load measurements are combined with reconstructed volumes thanks to computed tomography. Local and global residuals have been constructed in which each measured or predicted result is compared to the corresponding noise level. In 4D registrations, the mechanical finite element code was utilized in a non-intrusive way to generate elastoplastic fields. This type of approach is very generic and can be deployed under very different conditions.
Within the proposed framework, uncertainty quantifications have been carried out for the measured kinematic degrees of freedom and the sought parameters. For incremental correlation procedures, it is shown that the finer the mesh, the more uncertain the measured quantities. This trend is opposite to that of verification procedures where refined meshes better capture mechanical fields. This limitation can be overcome by considering kinematic fields that are mechanically admissible, i.e., to integrated DVC (i.e., 4D mechanical correlation). Consequently, convergence studies can be conducted when dealing with experimental data.
For the analyzed tensile test on nodular graphite cast iron, unknown boundary conditions and plastic parameters have been identified. By performing sensitivity analyses, it is shown that only some of the parameters can be identified for a given uncertainty. This is due to the fact that only a very limited number of scans is available for the reported analyses. However, even under these very difficult conditions, it is possible to carry out identifications with rather low uncertainties. Further, the verification analyses also show that most of the errors are due to the fact that measured (i.e., noisy) boundary conditions are prescribed on the numerical model. Yet, the identification results are virtually insensitive to the quality of the mesh because the geometry of the tested sample is very simple. In the present analyses the boundary conditions were measured and directly applied to the finite element model to compute kinematic bases. However, the boundary conditions may themselves become part of the set of unknowns. The present framework is suitable to such extensions.
Alternative 4D routes may be considered as well. For instance, by resorting to projectionbased registration procedures [48], it is possible to reduce the number of projections needed to evaluate 3D kinematic fields. When extended to (many) more time steps, it would allow the experimentalist not to interrupt the test and acquire the radiographs on the fly. A truly 4D formulation would be required in which reconstruction and registration procedures would be coupled with mechanically admissible fields.
Another direction of progress can be envisioned thanks to the scale at which microtomography is performed, namely, the study of the mechanical behavior at the microscale. For the studied material, this would allow the behavior of the ferrite/pearlite matrix to be analyzed in conjunction with the debonding of the interface between the brittle nodules and their surrounding matrix. The same framework as that introduced herein can be used. However, the meshes will needed to be significantly finer and adapted to the microstructure [58].