Bayesian Identification of Elastic Constants in Multi-Directional Laminate from Moiré Interferometry Displacement Fields

Abstract

The ply elastic constants needed for classical lamination theory analysis of multi-directional laminates may differ from those obtained from unidirectional laminates because of three dimensional effects. In addition, the unidirectional laminates may not be available for testing. In such cases, full-field displacement measurements offer the potential of identifying several material properties simultaneously. For that, it is desirable to create complex displacement fields that are strongly influenced by all the elastic constants. In this work, we explore the potential of using a laminated plate with an open-hole under traction loading to achieve that and identify all four ply elastic constants (E ₁ , E ₂ , ν ₁₂ , G ₁₂ ) at once. However, the accuracy of the identified properties may not be as good as properties measured from individual tests due to the complexity of the experiment, the relative insensitivity of the measured quantities to some of the properties and the various possible sources of uncertainty. It is thus important to quantify the uncertainty (or confidence) with which these properties are identified. Here, Bayesian identification is used for this purpose, because it can readily model all the uncertainties in the analysis and measurements, and because it provides the full coupled probability distribution of the identified material properties. In addition, it offers the potential to combine properties identified based on substantially different experiments. The full-field measurement is obtained by moiré interferometry. For computational efficiency the Bayesian approach was applied to a proper orthogonal decomposition (POD) of the displacement fields. The analysis showed that the four orthotropic elastic constants are determined with quite different confidence levels as well as with significant correlation. Comparison with manufacturing specifications showed substantial difference in one constant, and this conclusion agreed with earlier measurement of that constant by a traditional four-point bending test. It is possible that the POD approach did not take full advantage of the copious data provided by the full field measurements, and for that reason that data is provided for others to use (as on line material attached to the article).

Appendices

Appendix 1: Proper Orthogonal Decomposition

The objective of this appendix is to provide an overview of the proper orthogonal decomposition method that is used for the dimensionality reduction of the full-fields. First the theoretical foundations of the method are presented followed by some results of its application to the plate with a hole problem.

Let us consider U ⁱ ∈ ℝⁿ, which is the vector representation of a field (e.g. displacement field). Note that n is usually several thousands. We seek, based on N sample vectors {U ⁱ}_{i = 1..N}, a reduced dimensional representation of the fields’ variations with some input parameters.

The aim of the proper orthogonal decomposition (POD) method is to construct an optimal, reduced dimensional basis for the representation of the sample vectors (also called snapshots). In the POD approach the snapshots need to have zero mean, if this is not the case the mean value needs to be subtracted.

We denote ${\left\{ {{\Phi_k}} \right\}_{{_{{k = 1..K}}}}}$ the vectors of the orthogonal basis of the reduced dimensional representation of the snapshots. The POD method seeks to find the basis vectors Φ _k that minimize the representation error:

$$ \min \frac{1}{2}\sum\limits_{{i = 1}}^N {\left\| {{U^i} - \sum\limits_{{k = 1}}^K {{\alpha_{{i,k}}}{\Phi_k}} } \right\|}_{{L2}}^2 $$

(6)

Because ${\left\{ {{\Phi_k}} \right\}_{{_{{k = 1..K}}}}}$ is an orthogonal basis, the coefficients α _i,k are given by the orthogonal projection of the snapshots onto the basis vectors. As a result we have the following reduced dimensional representation ${\widetilde{U}^i}$ of the vectors of the snapshot set:

$$ {\widetilde{U}^i} = \sum\limits_{{k = 1}}^K {{\alpha_{{i,k}}}{\Phi_k} = \sum\limits_{{k = 1}}^K {\left\langle {{U^i},{\Phi_k}} \right\rangle {\Phi_k}} } $$

(7)

The reduction in dimension is from N to K. The truncation order K needs to be selected such as to maintain a reasonably small error in the approximate representations ${\widetilde{U}^i}$ of U ⁱ. Selecting such a K is always problem specific and an error criterion is given further down.

The main advantage of the POD method is that it provides a simple procedure for constructing the basis from the samples {U ⁱ}_{i = 1..N}. The procedure guarantees that for a given truncation order we cannot find any other basis that better approximates the snapshots subspace.

The basis ${\left\{ {{\Phi_k}} \right\}_{{_{{k = 1..K}}}}}$ is constructed using the following matrix:

$$ X = \left( {\matrix{ {U_1^1} \hfill & \ldots \hfill &{U_1^N} \hfill \\ \vdots \hfill & \ddots \hfill & \vdots \hfill \\ {U_n^1} \hfill & \ldots \hfill &{U_n^N} \hfill \\ }<!end array> } \right) $$

(8)

The vectors ${\left\{ {{\Phi_k}} \right\}_{{_{{k = 1..K}}}}}$ are then obtained by the singular values decomposition of X, or equivalently by calculating the eigenvectors of the matrix XX ^T. The singular values decomposition allows writing that:

$$ X = \Phi \sum {{\Lambda^T}} $$

(9)

where Φ is the matrix of the column vectors Φ _k . The svd() function in Matlab was used here for the singular value decomposition.

A truncation error criterion ε is then defined by the sum of the error norms as shown in equation 10.

$$ \sum\limits_{{i = 1}}^N {\left\| {{U^i} - \sum\limits_{{k = 1}}^K {{\alpha_{{i,k}}}{\Phi_k}} } \right\|}_{{L2}}^2 \leqslant \varepsilon \sum\limits_{{i = 1}}^N {\left\| {{U^i}} \right\|}_{{L2}}^2 $$

(10)

where $\varepsilon = 1 - \left( {{{{\sum\limits_{{j = 1}}^K {\sigma_j^2} }} \left/ {{\sum\limits_{{j = 1}}^N {\sigma_j^2} }} \right.}} \right)$, and σ _j are the diagonal terms of the diagonal matrix Σ. For a derivation of this criterion and further details on POD the reader can refer to [21]. The variation of this criterion with the truncation order for the plate with a hole problem is provided in Table 3 of the main section.

The error norm truncation criterion ε, while being a global error criterion, is relatively hard to interpret physically. Furthermore the criterion is based only on the convergence of the snapshots that served for the POD basis construction. However in most cases we will want to decompose a field that is not among the snapshots, so we also want to know the convergence of the truncation error in such cases.

Accordingly we chose to construct a different error measure based on cross validation. The basic idea of cross validation is the following: if we have N snapshots, instead of using them all for the POD basis construction we can use only N-1 snapshots and compute the error between the actual fields of the snapshot that was left out and its truncated POD decomposition. By successively changing the snapshot that is left out we can thus obtain N errors. The root mean square of these N errors, which we denote by CV_RMS, is then a global error criterion that can be used to assess the truncation inaccuracy.

In order to use the cross validation technique we need to define how to measure the error between two strain fields (the actual strain field and its truncated POD decomposition). We chose the maximum absolute difference between two fields. This maximum error is computed at each of the N (N = 200 here) cross validation steps and the root mean square leads to the global error criterion CV_RMS. Table 8 provides these values for different truncation orders. The relative CV_RMS error with respect to the value of the field where the maximum error occurs is also given in Table 8.

Table 8 Cross validation CV_RMS truncation error criterion

Next we provide a graphical representation of POD results in order to have a more intuitive understanding of the modal decomposition. First an illustration of the fields obtained for a particular snapshot (snapshot 1) is shown in Fig. 6, which provides an idea of the spatial variations and order of magnitude of the displacement fields. These fields were obtained with the following parameters: E ₁  = 202.2 GPa, E ₂  = 10.84 GPa, ν ₁₂  = 0.2142, G ₁₂  = 4.989 GPa, t = 0.1312 mm.

Fig. 6

U and V displacement fields for snapshot 1

The first four POD modes that we obtained are represented graphically in Figs. 7 and 8. We note that the first modes are relatively close (but not identical even though the differences cannot be seen by naked eye) to the typical U and V displacement fields (see Fig. 6). Furthermore we see that the modes have a more complicated shape with increasing mode number, as expected for a modal decomposition basis.

Fig. 7

First 4 POD modes for the U field

Fig. 8

First 4 POD modes for the V field

Appendix 2: Error Measures of the Response Surface Approximations

The error measures used to assess the quality of the response surface approximations (RSA) constructed in “Identification Results and Discussion” are given in Table 9 for the first four POD coefficients of the U fields and in Table 10 for those of the V fields. The second row gives the mean value of the POD coefficient across the design of experiments (DoE). The third row provides the standard deviation of the coefficients across the DoE, which gives an idea of the magnitude of variation in the coefficients. Row four provides R², the correlation coefficient obtained for the fit, while row five gives the root mean square error among the DoE points. The final column gives the cross validation PRESS error [22].

Table 9 Error measures for RSA of the U-field POD

Table 10 Error measures for RSA of the V-field POD

Comparing the error measures for each coefficient to their range of variation (i.e. standard deviations) we considered that the RSA are accurate enough to be used in the identification process, with the approximation error being negligible compared to the other sources of uncertainty.

Appendix 3: Bayesian Numerical Implementation

A Bayesian identification numerical procedure able to account for measurement uncertainty, modeling uncertainty as well as uncertainty in other model parameters was previously developed by the authors in [23] and will also be used here.

The numerical procedure uses Monte Carlo simulations for uncertainty propagation. While the use of Monte Carlo simulation has the advantage of propagating uncertainties represented by arbitrary probability distribution functions it has the drawback of being computationally expensive. The POD method and response surface methodology are thus used for reducing the cost associated with the construction of the likelihood function. A flowchart overview of the utilized procedures is presented in Fig. 9.

Fig. 9

Flow chart of the procedure used to calculate the likelihood function. Cost reduction is shown in green and dimensionality reduction in red

The likelihood function is computed point by point within the prior distribution’s support (truncation bounds) at a grid in the four-dimensional space of the material properties E = {E ₁ , E ₂ , ν ₁₂ , G ₁₂ }. We chose a 17⁴ grid, which is a compromise between convergence and computational cost considerations.

At each of the grid points, E is fixed and we need to evaluate the probability density function (PDF) of the POD coefficients, ${f_{{{{\alpha } \left/ {{E = {E^{{fixed}}}}} \right.}}}}(\alpha )$, at the point α = α ^measure. The PDF of the POD coefficients is determined by propagating through Monte Carlo with 4,000 simulations the uncertainties in the other model parameters (plate thickness here) and adding a sampled value of the normally distributed uncertainty in the POD coefficients resulting from measurement and modeling uncertainty, as described in the previous subsection.

Physical considerations showed that the resulting samples must be close to Gaussian so the samples were replaced by the normal distribution, having the sample mean and variance-covariance matrix. This Gaussian nature is due to the fact that the uncertainty resulting from the measurement noise is Gaussian and the uncertainty due to thickness is proportional to 1/h, which can in this case be well approximated by a normal distribution. The distribution ${f_{{{{\alpha } \left/ {{E = {E^{{fixed}}}}} \right.}}}}(\alpha )$ was then evaluated at the point α = α ^measure, leading to ${f_{{{{\alpha } \left/ {{E = {E^{{fixed}}}}} \right.}}}}({\alpha^{{measure}}})$. In this way we obtain a discretized likelihood function, which multiplied by the prior distribution gives us the posterior distribution of the elastic constants that we seek to identify.

At this point we want to make the following note. We found that the overall uncertainty on the POD coefficients is close to normal, which means that the Bayesian identification could have been treated within a purely analytical framework, thus avoiding the need for expensive Monte Carlo simulations. The analytical treatment would however have no longer been possible if uncertainties on more complex input parameters would have been considered leading to a clearly non-Gaussian distribution on the POD coefficients. In such a case the Monte Carlo simulations based approach would still work and this is the purpose why we developed this more general approach in [23]. For convenience we reused our already developed approach here, even though with several days of computational cost this approach is clearly not numerically the most efficient. More efficient numerical implementation are certainly possible and we provide the data files of the experimental results together with this paper such as to allow interested persons (including possibly ourselves at a future point) to test other implementations on the same experimental data.

Appendix 4: Data Structure of the Provided Measurements

The authors also provide attached to the online version of the paper the experimental displacement fields obtained. These are the displacement fields obtained based on the moiré interferometry images and the automated phase extraction procedure described in “Experiment”. The data is provided in matrix form as .txt files which can be opened with any text editor. The X.txt file provides the X coordinate (in millimeters) of each point where a displacement value is provided. Similarly the Y.txt file provides the Y coordinate. The U.txt file provides the displacements in micrometers at each point in the loading direction (Y direction). The V.txt file provides the displacements in the orthogonal direction (X direction). Note that the X and Y direction are related to the actual experiment axes (Fig. 2) not to the 1 and 2-driections provided in Fig. 1.