Physics-Enhanced Multi-fidelity Learning for Optical Surface Imprint

As shown in Figure 1A, one formal-shaped indenter (here is four-fold) is pushed into the material surface to create one minute crater (also called pile-up) (Figure 1C). There are two types of indentation techniques. The first technique is to apply the constant load ($P$) and only use the division of load ($P$) and projection area ($A$) of the crater to describe material’s property. That is called hardness ($H$).

While this method is convenient, the final measured result $H$ is not directly related to material parameters, i.e.,, the stress-strain relation. To acquire more information, the second method records the load evolution($P$) v.s.,indenting depth($h$) (Figure 2A right) and tries to map from the load-depth curve into material stress-strain relation ((Figure 2A left). However, many former research reveal that the load-depth curve may not be informative enough (Chen et al., 2007; Campbell et al., 2018). Here we show that the first technique (hardness measuring) plus optical profile can already make the inverse problem well-conditioned, even without the need for load-depth curve.

Our goal is to utilize the optical profilometer (Fainman et al., 1982) to measure the surface height map of the residual imprint. Specifically, as shown in Figure 1B, a beam of light from a single source is split by the interferometer into two separate beams. Each of these beams travel separate paths, one onto a reference surface and the other onto the surface to be measured. The beams are then recombined resulting in an interference pattern. An imaging device, usually a CCD array, is used to collect this information. By moving the interferometer vertically away from the measurement surface, the point at which this interference occurs can be found for each pixel of the CCD. By tracking the position of the interferometer during this process a 3D map (Figure 1D) of the surface can be formed.

Figure 2 schematically illustrates the problem set of our study. Figure 2A (left side) is a schematic diagram showing a typical stress-strain response of a power-law strain-hardening material which can be used for many engineering metallic materials. The elastic behavior follows Hook’s law, whereas the plastic response is approximated by different constitutive models (Hertelé et al., 2011). One assumption is the three-parameter Hollomon model (fitting parameters: $E$, $\sigma_{y}$, $n$), in which true stress $\sigma$ and true strain $\varepsilon$ are related as:

, while another assuming model is the four-parameter Ludwik model (fitting parameters: $E$, $\sigma_{y}$, $n$, $K$), displayed as:

, where E is the elastic modulus, $\sigma_{y}$ is the yield stress, $K$ is the work hardening coefficient, n is the work hardening exponent, and $\varepsilon_{p}$ is the equivalent plastic strain determined as $\varepsilon_{p}=\varepsilon-\dfrac{\sigma}{E}$. We find that the pre-assumed model sometimes fails to accurately fit the stress-strain curves measured in the real experiments, as shown in Appendix A. Later we will also turn into the study without pre-assuming the constitutive model and represent stress-strain curves via point-to-point linear connection.

Knowing the stress-strain relation of one material, theoretically we can uniquely determine the residual imprint and load-depth curve. This is referred as the forward problem. However, how to inversely determine the stress-strain curve from the residual imprint remains challenging.

To simulate the elasto-plastic behaviors of the system, the ABAQUS (Dassault Systèmes Simulia Corp.) software package is employed to conduct 2D and 3D FEM analysis (Reddy, 2019). For single 16-core CPU, each 2D job costs around 10 minutes, while each 3D job costs 6 hours.

Specific NN structure is displayed in Figure 2B,comprising several independent NNs connected by extra parameters. Each NN owns 6 hidden layers capable of learning the variations of pile-up features and hardness with respect to different stress-strain properties. Each hidden layer owns 32 neurons. ReLU is used as the activation function. The Adam optimizer (Kingma & Ba, 2014) is applied in training by setting the learning rate as 0.0001. Appendix B shows the dataset statistics, in which the parameters of FEM materials are wide enough to nearly cover all the metallic materials. We mainly use the dataset of Ludwik model to train the NNs in the following sections.

One fundamental question is whether the features we choose can guarantee a unique problem. Here we artificially define the feature extraction method based on curve characteristics, as illustrated in Appendix D. Specifically, we can extract three features from force curves based on the the loading curvature, initial unloading slope, and the ratio of residual unloading depth to maximum loading depth. We extract nine features from pile-up curves by finding the maximum height and calculating the volumes and weighted centres in varied parts. We also have tried using encoder-decoder structure to automatically output features, while the testing results show that both types of features perform similarly.

As shown in Figure 4A, we use 2D FEM data (Figure 3A) to build up an accurate forward prediction model, i.e.,, predicting the force and pile-up features based on input constitutive model parameters. The data number is actively augmented to ensure that the Mean Absolute Percentage Error ($MAPE$), defined as follows,

, for predicting each feature is below 2%. $T_{i}$ and $P_{i}$ are the true and prediction values of the $i^{th}$ data point, respectively. Each time 50 data points are added into the parameter range with large errors, and the total data number for forward prediction is 2450. Then we use this trained NN as a surrogate model to explore the informativeness of the indented features. To check whether one material owns siblings with quite similar features, we fix its parameters and iterate the parameters of the candidate material to continually decrease their feature differences. The iteration process uses typical BFGS optimization algorithm. The iteration stops when the $MAPE$ of all features are below the set limit or the iteration number exceeds the upper limit. Finally, we carry out FEM simulations to verify whether the hypothesized material siblings own similar features. Figure 4B1 and Figure 4B2 display two typical material siblings acquired from the above-mentioned workflow, in which the load-displacement curves are almost the same and the pile-ups reveal some difference at the highest parts.

Theoretically, we can distinguish two materials if the maximum difference among their features exceeds the possible variance, e.g.,, the experimental errors. Hence, we define the concept of distinguishing ratio as,

, where $F_{i}$ and $\bar{F}_{i}$ are the features from two material siblings, respectively. $N$ depends on how many features considered. We then take a uniform grid (grid number = 5) on the parameter space of Ludiwk model, and test the uniqueness of 625 materials. Among these 625 materials, $N_{A}$ of them own material siblings with all the relative feature differences lower than the distinguishing ratio. Then we define the non-unique ratio at this specific distinguishing ratio as $N_{A}/625$, displaying the possibility to encounter unique problems. Figure 4C shows the evolution of the non-unique ratio with the distinguishing ratio employing different features. It reveals that the unique problem is quite severe if only three force features and hardness are input as features. However, the non-uniqueness is largely mitigated if we also include the information of pile-up. Even if we choose only three features from pile-up combined with hardness, the performance is still much better than the pure force case when the distinguishing ratio is lower than 6%. Meanwhile, we find that the non-unique ratios with nine pile-up features and three force features are close to the case with nine pile-up features and hardness, both observing slight increases when the distinguishing ratio is higher than 8%. This hints that when including the information of pile-up, maybe we can represent the information of load-displacement relation with only hardness. We then verify this hypothesis by doing the inverse training with nine pile-up features and hardness. The $MAPE$ of each predicted parameter ($E$, $\sigma_{y}$, $n$, $K$) for the testing set decreases to lower than 5% when the total data number increases to 4000, as shown in Figure 4D.

We then discuss the transfer learning framework of our models. To decrease the burden of expensive data like 3D FEM simulations and experiments, we first build up a baseline model with many pure 2D FEM data, and then ameliorate the gap between 2D and 3D FEM models with some 3D FEM data. Finally, we calibrate the gap between 3D FEM simulations and experiments with several experimental data. The features we use in all the following conditions are based on nine pile-up features and hardness.

The numerical indenter simulated in the 2D axisymmetric FEM modeling is in the conical shape, in which the pile-up morphology can be represented by a one-dimensional (1D) curve, as shown in Figure 4B2. However, the indenters used in most hardness testing experiments are fourfold Vickers indenters without the axisymmetric property. Figure 5A and Figure 5B show two typical 2D images of pile-up morphologies in the experiment and 3D FEM simulation, respectively. Due to surface roughness and grain variations in real metals, the raw pile-up morphologies acquired from experiments are not perfectly smooth in heights (Figure 5A). To mitigate the height variations and transform the 2D image into the 1D curve, we divide the 2D image into many slim square strips and calculate the average height in each divided region. Then we acquire a plot of the averaged height v.s.,the distance to the imprint center and calculate related features, as shown in Figure 5C.

Since different shaped indenters will result to different pile-up morphologies and corresponding features, using the ML model trained by the 2D FEM data to predict the 3D FEM result will contain large errors. We first use 4000 2D FEM data to build up a 2D ML model and directly employ it for the prediction of 3D FEM results, denoted as $Y1$ in Figure 2. The green points in Fig. 4A display the prediction of $n$ in this case, and the errors are quite large. However, these predicted Y1 still incorporate the information of 2D FEM data, and they own a rough trend with the true values. To correct the wrong correspondence and transfer the model from 2D to 3D, here we set the predicted $Y1$ as the extra feature and put it together with the original pile-up features and hardness extracted from 3D FEM model to predict the target parameters ($E$, $\sigma_{y}$, $n$, $K$). Figure 6B shows the evolution of MAPE v.s.,the 3D data number for the case with and without the extra feature ($Y1$). We find the MAPE of the case incorporating 2D prediction will decrease to lower than 5% when the data number increases to 50, much lower than the case using pure 3D data. Figure 6A shows the correspondence of the predicted and targeted values based on this transferred model. The transfer learning framework decreases the requirement of expensive data.

After acquiring a 3D ML model, the next step is to transfer it into the real experimental ML model. However, as shown in Figure 5C and Figure 5F, the pile-up curves of 3D simulations and experiments are not perfectly consistent, even if the input stress-strain curves for simulations are acquired from the real materials. This inconsistency may be caused by many systematic errors, e.g., the indenter tip is not perfectly sharp, the over-simplified assumption of Poisson’s ratio $\nu$ to be 0.3 and friction coefficient $\mu$ to be 0.15, and the grain effects. The feature value differences between experiments and simulations are larger than 20% in some features, while the evolution trends with different materials are consistent. The destination target in this section is to use some types of materials for calibrating sim-to-real gap, and then employ this calibrated model to predict the stress-strain relation of other types of materials.

Under the constraints of the limited number of metal types in the real experiments or industries, one intuitive way to mitigate the sim-to-real gap is by incorporating some experimental data points with the 3D FEM data points together, and apply these merged data for the transfer learning process from 2D FEM model to 3D model, the same as the process mentioned in the section 2D to 3D simulation model transfer. This time the 3D model is not based on pure 2D/3D FEM data, since some experimental data are also used. Here in our experiments, each type of material has been measured 8 times repetitively. The number of data points chosen out of these 8 points for the transfer learning should be determined. We then plot the evolution of MAPE with the experiment data number in Figure 6C. The number of materials used for the training ranges from 1 to 3 (one-shot to three-shot), respectively. The MAPE are calculated from the other 20 testing materials. The specific choice of material types and experimental data points are randomly repeated for 20 times. With the experiment data number increasing, the MAPE will first decrease and then increase when the experiment data number exceeds 4. This phenomenon is reasonable since an overlarge number of data points of the same material tends to induce overfitting.

According to the above discussion, the transfer learning will be inefficient if the experimental data is rare and the sim-to-real gap is too large. In the above FEM simulations, the Poisson’s ratio $\nu$ and the friction coefficient $\mu$ are always assumed to be fixed values ($\nu$ = 0.3, $\mu$ = 0.15), which are common settings in most indentation models (Chen et al., 2007; Goto et al., 2019; Haj-Ali et al., 2008). However, our further study reveals that these two physical parameters will greatly impact the acquired features, as illustrated in Appendix E. We vary the $\nu$ to be (0.2, 0.3, 0.4) and the $\mu$ to be (0.05, 0.15, 0.25), and then calculate the evolution of pile-ups in these $3^{2}$ cases. For each type of features, the maximum difference among these 9 cases will range from 20% to 40%. The changing ratios will also vary among different materials. We then calculate the corresponding features of the materials (Al7075, SS430) and compare them to the experimental features. Figure 6E displays the total feature differences between simulations and experiments under varied $\nu$ and $\mu$. Among the nine combinations, the sim-to-real gap will be lower than others if the two physical parameters ($\nu$, $\mu$) = [(0.3, 0.15), (0.3, 0.05), (0.2, 0.15)]. Next, we build up three independent DNNs in which the input 3D FEM data are calculated by setting two physical parameters to be the above three values, respectively. Specifically, the parameters ($\nu$, $\mu$) are set as fixed values (0.3, 0.15) for 4000 2D FEM simulations, while changed to [(0.3, 0.15), (0.3, 0.05), (0.2, 0.15)] for each subset of 50 3D FEM data. The total number of 3D FEM data used increases to 150.

Then for each DNN, using the merged data of 3D FEM and experiments, one candidate value will be predicted, denoted as $Y3$, $Y3^{\prime}$, $Y3^{\prime\prime}$ in three independent subsets. These three DNNs form into a committee and finally determine what the optimal prediction value is. The specific weights for each DNN are tuned as,

Here $Y4$ is the final predicted value in our model and ($\alpha 1$, $\alpha 2$) are two parameters to be determined by experimental data. To stabilize the training, we use $Sigmoid$ function to constrain the values of ($\alpha 1$, $\alpha 2$) between 0 and 1. We use the experimental data of SS430 and Al7075 as the training data and find that the training error reaches minimum when $\alpha 1=0.19$, $\alpha 2=0.33$. Through forming a committee under varied physical parameters and fine-tuning the relative weights, the $\nu$ and $\mu$ represented in simulations can be closer to the experimental condition and other systematic biases such as tip radius effects of nominally sharp indenters can also be actively mitigated. More specifically, the reason why the settings [(0.3, 0.15), (0.3, 0.05), (0.2, 0.15)] are closer to experiments is not necessarily owing to the closer $\nu$ and $\mu$ values as that in experiments, but that these settings possibly offset the systematic biases caused by other factors. This whole process is named as ‘physical-boosting’ since the intuition is based on tuning two physical parameters in the simulations. In summary, the experimental training data has been used three times throughout the ML model, i.e., the choice of physical parameters with three least sim-to-real gaps, combining with 3D simulation data to train DNNs, fine-tuning the relation among three independent DNNs. We use this final model to predict the stress-strain relation of the left 20 other types of materials, serving as the testing set (SS304 (Figure 6F1) and Al6061 (Figure 6F2) as the examples shown). The predicted stress-strain curves are satisfyingly close to the real curves acquired by tensile testing (the mean relative errors of predicted stresses are 3.4% across all testing materials). The full accuracy table is in Appendix C.

During forming the committee of three ($\nu$, $\mu$) pairs, the initial grid search for choosing three closest ones is indispensable, since the member with quite large sim-to-real gap will even undermine the model performance.

To make the prediction more tenable and represent a much wider material space, here we add another kind of material model into the training and testing dataset, called ‘point-stress model’. In this model, the linear elastic part is determined by elastic modulus $E$ and yield stress $\sigma_{y}$, while the nonlinear plastic part is represented by nine strain-stress points, i.e., [$(\varepsilon_{i},\sigma_{i}),i=2,3,…10$], as shown in Figure 7A. The whole tensile curve is estimated between the points by linear interpolation. The yield strain $(\varepsilon_{y}=\dfrac{\sigma_{y}}{E})$ and the maximum strain $(\varepsilon_{10}=0.3)$ are then limits of the total interval of the strain values $\varepsilon_{i}$ to be estimated on the tensile curve. From this interval, the intermediate strain values, representing the positions of the $2^{nd}$ to $9^{th}$ points to estimate, are calculated by adopting a geometric progression. This choice of progression serves to obtain a higher density of points at lower strains, aiming to capture more features near the yield stress. The specific equation is defined as,

Here the common ratio $q$ ranges from 1.1 to 1.5 in different materials. A larger $q$ corresponds to denser aggregation of points near yield stress. After determining strain points $\varepsilon_{i}$, the corresponding stress values $\sigma_{i}$ are randomly generated under the constraint of softened-hardening behavior, displayed as,

In summary, the variables to be determined in the point-stress model are [$q,E,\sigma_{y},\sigma_{2},\sigma_{3},…,\sigma_{10}$].

Then the dataset for the pointwise stress-strain prediction of 2D FEM simulations comprises 2000 trials of Ludwik model, 1000 trials of Hollomon model, and 1000 trials of the point-stress model. As for the inverse prediction, we also try to predict the discrete strain-stress points and linearly interpolate them into the whole tensile curve, in which the objective parameters are [$E,\sigma_{y},(\varepsilon_{1}^{*},\sigma_{1}^{*}),(\varepsilon_{2}^{*},\sigma_{2}^{*}),…,(\varepsilon_{m}^{*},\sigma_{m}^{*})$]. Here the total number $m$ of strain-stress points are artificially determined by humans, even can be numerous if the computing resource is enough. Meanwhile, the specific positions of strains can also be flexible, different from the points in the input point-stress model. We set 6 strain positions uniformly gridded on [$\dfrac{\sigma_{y}}{E},\dfrac{\sigma_{y}}{E}+0.01$], and 10 strain positions uniformly gridded on [$\dfrac{\sigma_{y}}{E}+0.02,0.2$]. Figure 7B shows the evolution of the average mean absolute errors of all the objective parameters with the iterations (epochs). Both the training error (blue) and the testing error (red) will decrease to less than 5% when epochs arrive 2000. The inner plot shows the corresponding comparison between the true stress and the predicted stress. The predicting accuracy is magically good, hinting that each part of the stress-strain tensile curve may have its unique influences on the final pile-up morphologies.

The remaining process of transferring the 2D model into the 3D model, and then into the experimental model is the same as that stated in the section: Sim-to-real model transfer (Physical-boosting). For the 3D simulation data, apart from the original 150 trials of Ludwik model and 60 trials of Hollomon model, 60 trials of the point-stress model are added. The final predicted strain-stress points for the two training materials and two testing materials in experiments are displayed in Figure 7C and Figure 7D, respectively. The overall prediction accuracy is satisfying.